時間項の離散化-Runge-Kutta法 天才か

未分類

前回、の続きでRunge-Kutta法という天才的なアルゴリズムをみていきます。前回同様、関数を

$$
\begin{align}
&\frac{d u}{d t} =  f(t,u) \tag{1} \label{eq:scale-factor-1}
\end{align}
$$

とします。つまり、一階の常微分方程式とします。流体の支配方程式は偏微分方程式として記述されるのですが、時間項の扱い時は常微分方程式として考えます。まあ別に時間の偏微分のままで考えても全く変わらないのですが、空間が入ってない方が考えやすいのでそうして考えます。そうしていい理由を説明します。ここで、1次移流拡散方程式を考えます。

$$
\begin{align}
&\frac{\partial u}{\partial t} =  – U \frac{\partial u}{\partial x} + \nu \frac{\partial^2 u}{\partial x^2} + f \tag{2} \label{eq:scale-factor-2}
\end{align}
$$

変数u=u(t,x)は、独立変数t と独立変数xの双方に依存するので、uは偏微分方程式として記述せざるを得ません。また、より一般的にするために、移流速度\(U=U(t, x)\)、拡散係数\(\nu =\nu (t, x)\)、および外力\( f=f(t,x)\)として、これらの係数も時間と空間依存するとします。ここで、空間についてだけ中心差分で離散化します。

$$
\begin{align}
&\frac{d u_i}{d t} =  – U_i \frac{u_{i+1} – u_{i-1}}{2 \Delta x} + \nu_i \frac{u_{i-1}-2u_i+u_{i+1}}{\Delta x^2} + f_i \tag{3} \label{eq:scale-factor-3}
\end{align}
$$

このように、空間を離散化した段階では、\( u_i = u_i (t) \)となり、xの依存はなくなり、偏微分方程式から、常微分方程式になります。係数についても、空間の一点についての値なので、\(U_i(t) , \nu_i(t)\)となって時間のみです。このように、流体の支配方程式を空間離散化した後は、時間の常微分方程式になります。つまりは、\( f(t,x_1,x_2,x_3) \)を、\( f(t) \)に直しただけのことで時間の偏微分と空間を離散化して時間の常微分にすることは全くの同値なのですが、Runge-Kuttaの場合は手順的にもこうしますし、分かりやすいのでこうします。


Runge-Kutta法

Runge-Kutta法だけはまた毛色が違うアルゴリズムになっています。
一応、Runge-Kutta法は陽解法と書きましたが、陰解法のRunge-Kutta法もあります。一般にRunge-Kuttaと言えば陽的Runge-Kuttaを指しますが、同じ発送で陰的Runge-Kuttaを作ることもできます。Runge-Kutta法がすごいのは、発散しない時間ステップ幅(\Delta t)の安定幅が高次になるに比例して広がっていくこと、比較的少ない計算で高い精度が得られるからです。また、係数が簡単なこと、4次までは比例して精度が伸びていくことから、4次のRunge-Kutta法は頻繁に使われます。まずは陽解法の一般的に指すRunge-Kutta法をやっていきます。まず、発想を理解します。冒頭にあった、陽解法や陰解法の基本的な考えがあったと思いますが、Runge-Kuttaの考え方はこうです。

細かく途中で刻んでいくイメージです。本当は途中の点などの重み付け平均なので違いますが、あくまでイメージです。それをうまく使って、精度を上げます。そして、関数が\eqref{eq:scale-factor-1}式のように、

$$
\begin{align}
&\frac{d u}{d t} =  f(t,u)
\end{align}
$$

となっているとします。fは当然、uの関数でもあり、tの関数でもあります。さて、Adoms-Bashforthの時のように、uをtについて、n周りで展開します。ここで、\( (u^{t+1}, u^t) \)のように書いてきたものを一旦、\( ((t + \Delta t), u(t)) \)などと書きます。意味は同じですが、Runge-Kuttaはこう書かないとやりずらい時もあるので説明のために一旦こうします。まあ別にどっちでやっても書き方が違うだけですが、最初の導入はこっちの方が僕は好きなだけです。

$$
\begin{align}
u(t+\Delta t) &= u(t) + \frac{\Delta t}{1!} \left.\frac{d u}{d t}\right\vert^t +\frac{\Delta t^2}{2!} \left.\frac{d^2 f}{d t^2}\right\vert^t +\frac{\Delta t^3}{3!} \left.\frac{d^3 u}{d t^3}\right\vert^t + \cdots  \\\\
&= u(t) + \frac{\Delta t}{1!} f(t,u) +\frac{\Delta t^2}{2!} \left.\frac{d f}{d t}\right\vert^t +\frac{\Delta t^3}{3!} \left.\frac{d^2 f}{d t^2}\right\vert^t + \cdots  \tag{4} \label{eq:scale-factor-4}
\end{align}
$$

この時、\eqref{eq:scale-factor-4}式の第3項以下を誤差項にすればオイラー陽解法になります。

$$
\begin{align}
u(t + \Delta t) &\simeq u(t) + \Delta t f(t,u) + O(\Delta t^2)  \tag{5} \label{eq:scale-factor-5}
\end{align}
$$

精度が2に見えますが、ここから一個(\Delta t\)で割るので、精度は1です。この場合の精度の上げ方として、\eqref{eq:scale-factor-4}式の第4項以下を誤差項にして、そのまま、

$$
\begin{align}
u(t + \Delta t) &\simeq u(t) + \Delta t (f(t,u) +\frac{\Delta t}{2} \left.\frac{d f}{d t}\right\vert^t ) + O(\Delta t^3)  \tag{6} \label{eq:scale-factor-6}
\end{align}
$$

とする方法が考えられます。あとはこれをうまく閉じていきます。ここで、uはtの関数であるので、f(t,u)をtで微分する際には、tによって変化するuの変化分も考えなければなりません。よって、f(t,u)の常微分である(\left.\frac{d f}{d t}\right\vert^t )は、こうなります。

$$
\begin{align}
\left.\frac{d f}{d t}\right\vert^n &= \frac{\partial f}{\partial t} \frac{dt}{dt} + \frac{\partial f}{\partial u} \frac{du}{dt}  \\\\
&= \frac{\partial f}{\partial t} + \frac{\partial f}{\partial u} f  \\\\
&= f_t + f_u f \tag{7} \label{eq:scale-factor-7}
\end{align}
$$

ここでは、各変数の偏微分を(f_t)などとして表しています。これは、そのままtの偏微分で書いても、uはtの関数なので、偏微分でも成り立つと思うのですが、ちょっと微妙なので、厳密性のために常微分にしておきました。なので、\eqref{eq:scale-factor-6}式は、

$$
\begin{align}
u(t + \Delta t) &\simeq u(t) + \Delta t (f(t,u) +\frac{\Delta t}{2} \left.\frac{d f}{d t}\right\vert^t  + O(\Delta t^3) \\\\
&= u(t) + \Delta t (f(t,u) +\frac{\Delta t}{2}(f_t + f(t,u) \ f_u) \tag{8} \label{eq:scale-factor-8}
\end{align}
$$

になります。さてここで、fが、(c_2 \Delta t , a_{21} \Delta u)だけ進んだ点の2変数のTaylor展開を考えていきます。係数がちょっと変なのは後に一般化するためです。2変数のTaylor展開を一応書いておくと、

$$
\begin{align}
f (x+ \Delta x, y+ \Delta y) &= f(x,y) + \frac{1}{1!} (\Delta x \frac{\partial }{\partial x} +\Delta y \frac{\partial }{\partial y} ) f(x,y) + \frac{1}{2!} (\Delta x \frac{\partial }{\partial x} +\Delta y \frac{\partial }{\partial y} )^2 f(x,y) + \frac{1}{3!} (\Delta x \frac{\partial }{\partial x} +\Delta y \frac{\partial }{\partial y} )^3 f(x,y) + \cdots  \tag{9} \label{eq:scale-factor-9}
\end{align}
$$

今回は、uはtの関数なので、(\Delta u = \Delta t f(t,u))と近似できます。Taylor展開をさらにやるので、2次以下のものはなくなると考え、気にしません。2変数のTaylor展開を使うと、こうなります。

$$
\begin{align}
f (t+c_2 \Delta t, u(t)+ a_{21} \Delta t f(t,u)) &= f(t,u) + c_2 \Delta t \frac{\partial f}{\partial t} + a_{21} \Delta t f(t,u) \frac{\partial f}{\partial u} + O(\Delta t^2) \\\\
&= f(t,u) + c_2 \Delta t f_t + a_{21} \Delta t f(t,u) f_u + O(\Delta t^2) \tag{10} \label{eq:scale-factor-10}
\end{align}
$$

なんとなく見えてきたでしょうか?このTaylor展開で算出した式と、\eqref{eq:scale-factor-8}式の括弧の中は変数を合わせれば一致しそうです。ここで単に\(c_2=1/2, a_{21}=1/2\)だと思ってやると、修正オイラー法という方法になります。やってみます。

$$
\begin{align}
u(t + \Delta t) &= u(t) + \Delta t (f(t,u) +\frac{\Delta t}{2}(f_t + f(t,u) \ f_u) \\\\
&=u(t) + \Delta t ( f (t+c_2 \Delta t, u(t)+ a_{21} \Delta t f(t,u) ) )\\\\
&= u + \Delta t f (t+ \frac{\Delta t}{2}, u+ \frac{f \Delta t}{2} )\tag{11} \label{eq:scale-factor-11}
\end{align}
$$

このままじゃ少し足りないので、もう一工夫加えます。\(f(t,u)\)も使うことを考えます。

$$
\begin{align}
& b_1 f(t,u) + b_2 f(t+c_2 \Delta t, u(t)+ a_{21} \Delta t f(t,u)) \\\\
&= b_1 f (t,u) + b_2 ( f(t,u) + c_2 \Delta t f_t + a_{21} \Delta t f(t,u) f_u + O(\Delta t^2 ) \\\\
&= (b_1+b_2) f(t,u) + \Delta t(b_2 c_2 f_t + b_2 a_{21} f(t,u) \ f_u) \tag{12} \label{eq:scale-factor-12}
\end{align}
$$

この時、この項は\(b_1+b_2=1, b_2 c_2=1/2, b_2 a_{21}=1/2\)を満たせば、\eqref{eq:scale-factor-8}式の括弧の中と一致します。これが、2次のRunge-Kutta法の条件式です。先ほどのは修正オイラー法はRunge-Kutta法の一種ではありますが、一般には修正オイラー法と呼ばれる解法です。このように、Runge-Kutta法では2次では変数4つに対して、条件式が3つと解が無数にある状況になっています。なので亜種がたくさんあります。やはり簡単なのが正義なので、簡単になる解を決めていきます。そこで、今度は\(b_1=b_2=1/2, c_2=a_{21}=1\)とします。すると、

$$
\begin{align}
u(t + \Delta t) &= u(t) + \Delta t (f(t,u) +\frac{\Delta t}{2}(f_t + f(t,u) \ f_u) \\\\
&= u(t) + \Delta t ( (b_1+b_2) f + \Delta t( b_2 c_2 f_t+b_2 a_{21}f f_u) ) \\\\
&= b_1 f(t,u) + b_2 f(t+ c_2 \Delta t, u(t) + a_{21} \Delta t f(t,u) ) )\\\\
&= u(t) + \frac{\Delta t}{2} ( f(t,u) +  f(t+\Delta t, u+  \Delta t f(t,u) ) \tag{13} \label{eq:scale-factor-13}
\end{align}
$$

となります。これが一般的に指す2次のRunge-Kuttaです。これを分けて、

$$
\begin{align}
&s_1 = f(t,u) \\\\
&s_2 = f(t + \Delta t, u+ \Delta t s_1) \\\\
&u(t + \Delta t) =u(t) + \frac{\Delta t}{2} ( s_1 + s_2 ) \tag{14} \label{eq:scale-factor-14}
\end{align}
$$

という式で使われます。これは、まず\(s_1\)を出してから、\(s_2\)を出すというアルゴリズムのため、このように記述されることが多いです。なかなかイメージしにくいかもしれません。一旦、熱伝導方程式で考えてみます。空間を離散化して、

$$
\begin{align}
& \frac{\partial T}{\partial t} =  \alpha \frac{\partial^2 T}{\partial x^2} \\\\
& \frac{\partial T_i}{\partial t} =\alpha \frac{T_{i-1} – 2T_i+T_{i+1}}{\Delta x^2}  = f(t, T_i) \tag{15} \label{eq:scale-factor-15}
\end{align}
$$

となりますが、他のi-1やi+1の時にも成り立つので、

$$
\begin{align}
& f(T_{i-1}^n, T_i^n, T_{i+1}^{n}) = \alpha \frac{T_{i-1} – 2T_i+T_{i+1}}{\Delta x^2} \tag{16} \label{eq:scale-factor-16}
\end{align}
$$

と書くべきでしょう。tはややこしいので、Tの値を\(T_i^n\)などとして、t時点であることを表現します。

$$
\begin{align}
s_{1,i} &= f(T_{i-1}^n , T_i^n, T_{i+1}^n) = \alpha \frac{T_{i+1}^n – 2T_i^n +T_{i-1}}{\Delta x^2}^n  \\\\
s_{2,i} &= f(T_{i-1}^n + \Delta t s_{1,i-1} , T_i^n + \Delta t s_{1, i}, T_{i+1}^n + \Delta t s_{1, i+1}) \\\\
& =  \alpha \frac{ (T_{i-1}^n + \Delta t s_{1,i-1} )  – 2 ( T_i^n + \Delta t s_{1, i}) + ( T_{i+1}^n + \Delta t + s_{1, i+1} ) }{\Delta x^2} \\\\
T_i^{n+1} &= T_i^n + \frac{\Delta t}{2} ( s_{1,i} + s_{2,i} )
\tag{17} \label{eq:scale-factor-17}
\end{align}
$$

という形になります。空間が離散化されているものに対して考えるので少しややこしいですね。でもイメージは掴めたでしょうか?\(u+\Delta t s_1\)の部分で、tの時点はうまく表現されているので、tのところはあまり気にしないことです。また、\(s_{1,i} \)を先に全部求めてから、\(s_{2.i}\)を算出しにいきます。 そして、これを一般化していきます。2段2次の陽的Runge-Kutta法は以下の式で与えられます。

$$
\begin{align}
&s_1 = f(t,u) \\
&s_2 = f(t + c_2 \Delta t, u+ \Delta t  a_{21} s_1) \\
&u(t + \Delta t) =u(t) + \Delta t( b_1 s_1 + b_2 s_2 ) \tag{18} \label{eq:scale-factor-18}
\end{align}
$$

ここで、\(s_1, s_2\)などのように、何個算出するかを段と呼び、その時の精度と組み合わせて、2段2次陽的Runge-Kutta法と表現します。そして、s段の陽的Runge-Kutta法の到達可能精度p(s)は以下になることが知られています。

s12345678910\(s \geq 10\)
p(s)1234456677\(p(s) \leq s-2\)

精度と計算量の兼ね合いから、4段4次精度公式がもっとも効率がいいことがわかると思います。なので最も使われます。4段4次精度のものを求めて行こうかと思います。 2段の時と発想は同じです。

$$
\begin{align}
&s_1 = f(t,u) \\
&s_2 = f(t + c_2 \Delta t , u+ a_{21} s_1 \Delta t ) \\
&s_3 = f(t + c_3 \Delta t , u+ a_{31} s_1 \Delta t + a_{32} s_2 \Delta t ) \\
&s_4 = f(t + c_4 \Delta t , u+ a_{41} s_1 \Delta t + a_{42} s_2 \Delta t + a_{43} s_3 \Delta t) \\
&u(t + \Delta t) =u(t) + \Delta t ( b_1s_1 + b_2s_2 + b_3s_3 + b_4s_4 ) \tag{19} \label{eq:scale-factor-19}
\end{align}
$$

とおきます。ここまでの流れは先ほど、2次で、\(b_1, b_2, c_2, a_{21}\)を使って置いたものと全く同じです。ここで、\(a_{22}, a_{23}, a_{33}\)などがないのは陽的だからです。ここで、こういった項が入ると、陰的または半陰的のRunge-Kutta法になります。後は、元のTaylor展開の式に4次の項まで一致させます。Taylor展開のために、\(f = f(t, u(t))\)の時間微分を求めていきます。

$$
\begin{align}
&\frac{du}{dt} = f, \quad \frac{d^2 u}{dt^2}=Df, \quad \frac{d^3 u}{dt^3} = D^2f + (Df)f_u , \\
& \frac{d^4}{dt^4} = D^3f+(D^2 f)f_u + (Df)f_u^2  + 3(Df)Df_u \tag{20} \label{eq:scale-factor-20}
\end{align}
$$

と書けます。ただし、\(f_u = \frac{\partial f}{\partial u}\), および、

$$
\begin{align}
&D = \frac{\partial }{\partial t} + f \frac{\partial }{\partial u}, \quad D^2 = \frac{\partial^2}{\partial t^2} + 2f \frac{\partial^2}{\partial t \partial u} + f^2 \frac{\partial^2}{\partial u^2}, \\
& D^3 = \frac{\partial^3}{\partial t^3} + 3f \frac{\partial^3}{\partial t^2 \partial u} + 3f \frac{\partial^3}{\partial t \partial u^2} + 3f^3 frac{\partial^3}{\partial u^3} \tag{21} \label{eq:scale-factor-21}
\end{align}
$$

とします。このまま、u(t)のTaylor展開の高次の項まで求めていきます。

$$
\begin{align}
u(t + \Delta t) &= u(t) + \Delta t (f(t,u) +\frac{\Delta t^2}{2!} \left.\frac{d f}{d t}\right\vert^n +\frac{\Delta t^3}{3!} \left.\frac{d^2 f}{d t^2}\right\vert^n + \frac{\Delta t^4}{4!} \left.\frac{d^3 f}{d t^3}\right\vert^n  + O(\Delta t^5) \\
&= u(t) + \Delta t (f(t,u) +\frac{\Delta t^2}{2!}(Df) + \frac{\Delta t^3}{3!}(D^2f + (Df)f_u) \\
\hspace{30px} &+ \frac{\Delta t^4}{4!}(D^3f + (D^2f)f_u + (Df)f_u^2 + 3(Df)Df_u  ) + \cdots \tag{22} \label{eq:scale-factor-22}
\end{align}
$$

そしてこんな感じで(\Delta t^4)までの係数を比較すると、13個の未知数に対する、11個の方程式が得られます。ゴリゴリ書いてみます。

$$
\begin{align}
&c_2 = a_{21} \\
&c_3 = a_{31} \\
&c_4 = a_{41} + a_{42} + a_{43} \\
&b_1 + b_2 + b_3 + b_4 = 1 \\
&b_2c_2 + b_3c_3 + b_4c_4 = 1/2 \\
&b_2c_2^2 + b_3c_3^2 + b_4c_4^2 = 1/3 \\
&b_2c_2^3 + b_3c_3^3 + b_4c_4^3 = 1/4 \\
&b_3a_{32}c_2 + b_4(a_{42}c_2 + a_{43}c_3) = 1/6 \\
&b_3a_{32}c_2^2 + b_4(a_{42}c_2^2 + a_{43}c_3^2) = 1/12 \\
&b_3a_{32}c_2c_3 + b_4(a_{42}c_2 + a_{43}c_3)c_4 = 1/8 \\
&b_4a_{43}a_{32}c_2 = 1/24
\tag{23} \label{eq:scale-factor-23}
\end{align}
$$

そして、一番簡単で有名解として、1/6公式というものがあり、次式で与えられます。

$$
\begin{align}
&s_1 = f(t,u) \\
&s_2 = f(t +  \Delta t / 2 , u +  s_1 \Delta t /2 ) \\
&s_3 = f(t +  \Delta t /2 , u +  s_2 \Delta t /2 ) \\
&s_4 = f(t +  \Delta t , u+  s_3 \Delta t) \\
& u(t + \Delta t) =u(t) + \frac{\Delta t}{6} ( s_1 + 2s_2 + 2s_3 + s_4 ) \tag{24} \label{eq:scale-factor-24}
\end{align}
$$

覚えやすくていいですね。導出過程をちゃんとみたい人はRunge-Kutta法についてのノート(早川尚男)もしくは、流体伝熱基礎講座第4回をご覧ください。他にも1/8公式なんてのもあります。ここで、より一般化していきます。完全に一般化すると、

$$
\begin{align}
&s_i = f \bigl( t + c_i \Delta t , u(t) + \displaystyle \sum_{j=1}^{s} a_{ij} s_j \Delta t \bigr), \quad (1 \leq i \leq s) \\
&u(t + \Delta t) =u(t) + \Delta t \displaystyle \sum_{i=1}^{s} b_i s_i \tag{25} \label{eq:scale-factor-25}
\end{align}
$$

ただし、係数に対し

$$
\begin{align}
&c_i = \displaystyle \sum_{j=1}^{s} a_{ij}, \quad (1 \leq i \leq s) \tag{26} \label{eq:scale-factor-26}
\end{align}
$$

という条件が付きます。これは今までも実は暗に満たされていたのですが、係数を比較すると必ず出てくる条件です。そして、係数の表現方法として便利なので、ブッチャー配列を導入します。つまり、

$$
\begin{align}
& \begin{array}{c|cccc} c_1 & a_{11} & a_{12} & \cdots & a_{1s}  \\ c_2 & a_{21} & a_{22} & \cdots & a_{2s} \\  \vdots & \vdots & \vdots & \ddots & \vdots  \\ c_s & a_{s1} & a_{s2} & \cdots & a_{ss} \\ \hline  & b_1 & b_2 & \cdots & b_s  \end{array}
\end{align}
$$

ブッチャー行列だと、先ほどのcの条件がわかりやすくて便利ですね。横にみて、cとaの合計が一致してないとだめという形です。今まで、陽的Runge-Kuttaとしてやってきたものは、\(i \leq j \)に対して\(a_{ij}=0\)、つまりaの上三角成分を0にしたものです。つまり、

$$
\begin{align}
& \begin{array}{c|ccccc} c_1 & 0 & \cdots & \cdots & \cdots & 0  \\ c_2 & a_{21} & 0 & \cdots & \cdots & 0 \\ c_3 & a_{31} & a_{32} & 0 & \cdots & 0 \\  \vdots & \vdots & \vdots & \ddots & \ddots & \vdots  \\ c_s & a_{s1} & a_{s2} & \cdots & a_{ss-1} & 0 \\ \hline  & b_1 & b_2 & \cdots & b_{s-1} & b_s  \end{array}
\end{align}
$$

となっているものです。これが陽的になる理由はなんとなく分かるでしょうか?陽的とは逆行列(連立方程式)を解くというステップが必要のないそのままの計算でできる方法をさします。先ほどの4段4次のRunge-Kuttaの1/6公式はブッチャー行列で表現するとこうなります。

$$
\begin{align}
& \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0  \\ 1/2 & 1/2 & 0 & 0 & 0 \\  1/2 & 0 & 1/2 & 0 & 0  \\ 1 & 0 & 0 & 1 & 0 \\ \hline  & 1/6 & 1/3 & 1/3 & 1/6  \end{array}
\end{align}
$$

そんで、1/8公式と呼ばれるものはこんな感じです。

$$
\begin{align}
& \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0  \\ 1/3 & 1/3 & 0 & 0 & 0 \\  2/3 & -1/3 & 1 & 0 & 0  \\ 1 & 1 & -1 & 1 & 0 \\ \hline  & 1/8 & 3/8 & 3/8 & 1/8  \end{array}
\end{align}
$$

最後に、\eqref{eq:scale-factor-25}式のような形で表現してきましたが、もう一つ表現方法があります。\eqref{eq:scale-factor-25}式もう一度書いておきます。

$$
\begin{align}
&s_i = f \bigl( t + c_i \Delta t , u(t) + \displaystyle \sum_{j=1}^{s} a_{ij} s_j \Delta t \bigr) \\
&u(t + \Delta t) =u(t) + \Delta t \displaystyle \sum_{i=1}^{s} b_i s_i  \tag{27} \label{eq:scale-factor-27}
\end{align}
$$

これに対し、
$$
\begin{align}
&u^{(i)} = u^n + \Delta t \displaystyle \sum_{j=1}^{s} a_{ij} f \bigl( t + c_i \Delta t , \ u(t) + u^{(j)} \bigr), \quad (1 \leq i \leq s) \\
&u^{n+1} =u^n + \Delta t \displaystyle \sum_{i=1}^{s} b_i f \bigl( t + c_i \Delta t ,  \ u^{(i)} \bigr) \tag{28} \label{eq:scale-factor-28}
\end{align}
$$

という式も構築できます。ただ、関数の部分を出しただけです。導出が簡単なのと、最初が理解しやすいので前者の方を使っていましたが、後者の方がアルゴリズム的に書きやすい人も多いので、どちらも知っておくといいと思います。


陰的Runge-Kutta法

では次に、陰的Runge-Kutta法を見ていきます。陰的Runge-Kutta法ではある場所で\(a_{ij} \neq 0\)になります。そして、陰的になった式で連立方程式を解いていくことになります。また、\(i \\lt j \)に対して、\(a_{ij}=0\)とするものに対して半陰的Runge=Kutta法と呼びます。s段の半陰的Runge-Kuttaの到達可能精度は、\(p(s) = s+1\)になることが知られています。また、陰的Runge-Kuttaの到達可能精度は\(p(s) = 2s\)になります。とりあえず、2段陰的Runge-Kuttaに対してブッシャー配列と計算式は、

$$
\begin{align}
& \begin{array}{c|cc} c_1 & a_{11} & a_{12}  \\ c_2 & a_{21} & a_{22} \\ \hline  & b_1 & b_2   \end{array}
\end{align}
$$

および、

$$
\begin{align}
&s_1 = f(t + c_1 \Delta t, u + a_{11} s_1 \Delta t + a_{12} s_2 \Delta t) \\
&s_2 = f(t + c_2 \Delta t, u+ a_{21} s_1 \Delta t + a_{22} s_2 \Delta t) \\
&u(t + \Delta t) =u(t) + \Delta t( b_1 s_1 + b_2 s_2 ) \tag{29} \label{eq:scale-factor-29}
\end{align}
$$

もしくは、

$$
\begin{align}
&u^{(1)} = u^n + \Delta t \bigl( a_{11} f(t + c_1 \Delta t , u^{(1)}) + a_{12} f(t + c_2 \Delta t , u^{(2)}) \bigr) \\
&u^{(2)} = u^n + \Delta t \bigl( a_{21} f(t + c_1 \Delta t , u^{(1)}) + a_{22} f(t + c_2 \Delta t , u^{(2)}) \bigr) \\
&u^{n+1} =u^n + \Delta t( b_1 f(t + c_1 \Delta t , u^{1})  + b_2 f(t + c_2 \Delta t , u^{2}) ) \tag{30} \label{eq:scale-factor-30}
\end{align}
$$

となり、\(s_1\)と\(s_2\)の連立方程式、もしくは\(u^{(1)}\)と\(u^{(2)}\)の連立方程式を解いた後に\(u^{n+1}\)を計算します。でも、これ計算量が半端ないの分かるでしょうか?1つの式に対しても陰的になっているので空間のいろんな点が入った行列の連立方程式をときますが、それが2つあってさらにそれ同士の連立方程式なんで、半端ないっす。やばいっす。計算量は莫大です。これに対し、半陰的Runge-Kuttaを考えます。半陰的の時は、

$$
\begin{align}
& \begin{array}{c|cc} c_1 & a_{11} & 0  \\ c_2 & a_{21} & a_{22} \\ \hline  & b_1 & b_2   \end{array}
\end{align}
$$

$$
\begin{align}
&s_1 = f(t + c_1 \Delta t, u + a_{11} s_1 \Delta t ) \\
&s_2 = f(t + c_2 \Delta t, u+ a_{21} s_1 \Delta t + a_{22} s_2 \Delta t) \\
&u(t + \Delta t) =u(t) + \Delta t( b_1 s_1 + b_2 s_2 ) \tag{31} \label{eq:scale-factor-31}
\end{align}
$$

もしくは、

$$
\begin{align}
&u^{(1)} = u^n + \Delta t a_{11} f(t + c_1 \Delta t , u^{(1)})  \\
&u^{(2)} = u^n + \Delta t \bigl( a_{21} f(t + c_1 \Delta t , u^{(1)}) + a_{22} f(t + c_2 \Delta t , u^{(2)}) \bigr) \\
&u^{n+1} =u^n + \Delta t( b_1 f(t + c_1 \Delta t , u^{1})  + b_2 f(t + c_2 \Delta t , u^{2}) ) \tag{32} \label{eq:scale-factor-32}
\end{align}
$$

となります。この場合は、\(s_1\)と\(s_2\)が単独で連立方程式、もしくは\(u^{(1)}\)と\(u^{(2)}\)が単独で連立方程式を解いた後に\(u^{n+1}\)を計算します。なので、一段あたりの負荷がAdoms-Molton法と同等です。ハンパねえ、、、陰的Runge-Kutta、バケモン。代表的なものを書いてみます。IRK = Implicit-Runge-Kutta法(陰的Runge-Kutta)、SIRK = Semi-Implicit-Runge-Kutta(半陰的Runge-Kutta)とします。

オイラー陰解法(1段1次IRK法)  陰的中点法(1段2次IRK法)  Crank-Nicolson法(2段2次 SIRK法)

$$
\begin{align}
& \begin{array}{c|c} 1 & 1  \\ \hline  & 1  \end{array} \hspace{70pt} \begin{array}{c|c} 1/2 & 1/2  \\ \hline  & 1  \end{array} \hspace{70pt} \begin{array}{c|cc} 0 & 0 & 0  \\ 1 & 1/2 & 1/2 \\ \hline  & 1/2 & 1/2  \end{array} 
\end{align}
$$

2段4次IRK法                                                                     2段3次SIRK法(ノルセット法)

$$
\begin{align}
& \begin{array}{c|cc} (3-\sqrt{3})/6 & 1/4 & (3-2 \sqrt{3})/12 \\ (3+\sqrt{3})/6 & (3+2 \sqrt{3})/12 & 1/4 \\ \hline  & 1/2 & 1/2  \end{array} \hspace{40pt} \begin{array}{c|cc} (3+\sqrt{3})/6 & (3+\sqrt{3})/6 & 0 \\ (3-\sqrt{3})/6 & -3 \sqrt{3} & (3+\sqrt{3})/6 \\ \hline  & 1/2 & 1/2  \end{array}
\end{align}$$3段6次IRK法$$\begin{align}
& \begin{array}{c|ccc} (5-\sqrt{15})/10 & 5/36 & 2/9 – \sqrt{15}/15 & 5/36 – \sqrt{15}/30 \\ 1/2 & 5/36 + \sqrt{15}/24 & 2/9 & 5/36 – \sqrt{15}/24 \\ 5+\sqrt{15})/10 & 5/36 + \sqrt{15}/30 & 2/9 + \sqrt{15}/15 & 5/36 \\ \hline  & 5/18 & 4/9 & 5/18  \end{array}
\end{align}
$$

Crank-Nicolson法とは、2次Adoms-Molton法のことです。1段や2段のものが既存の方法と一致するのは分かるでしょうか?Crank-Nicolson法だけやってみます

$$
\begin{align}
&u^{(1)} = u^n  \\\\
&u^{(2)} = u^n + \frac{\Delta t}{2} f(t^n, u^{(1)}) + \frac{\Delta t}{2}f(t^n+\Delta t, u^{(2)}) \\\\
&u^{n+1} =u^n + \frac{\Delta t}{2} f(t^n, u^{(1)})  + \frac{\Delta t}{2} f(t^n +\Delta t, u^{(2)}) ) \tag{33} \label{eq:scale-factor-33}
\end{align}
$$

と書き直すことができます。この場合、\(u^{(1)} = u^n), (u^{(2)} = u^{n+1}\)になっていて、Crank-Nikolson法(2次Adoms-Molton法)の、

$$
\begin{align}
&u^{n+1} =u^n + \frac{\Delta t}{2}( f^{n+1} + f^n) \tag{34} \label{eq:scale-factor-34}
\end{align}
$$

になることが分かると思います。まあ陰的はあまり使うことがないでしょう。計算負荷やばいです。しかし天才的な方法ですね。普通に感動します。なんか天才が考えたヤヴァイ方法を使ってる感が出て好きですね、本当に。今回の内容は、なんども繰り返しますが流体伝熱基礎講座第4回を参考にさせていただきました。非常に素晴らしい資料で昔勉強した時にも今回も非常にお世話になりました。ありがとうございます。



[
前の記事へ]  [非圧縮性流体の目次へ]  [次の記事へ]


 

未分類

Posted by mamitsu