時間項の離散化-AB法とAM法 Adoms-Bashforth法とAdoms-Molton法で離散化

非圧縮性流体

前回、かなり雑にやってしまった、時間項の離散化をしっかりとTaylor展開を使ってやっていこうと思います。

前回、こんな感じで近似する3通りが考えられるといったと思います。それぞれ、名前としてはこんな感じです。

$$
\begin{align}
&f(t+\Delta t) -f(t) = \int_t^{t+\Delta t} F dt \\\\
&f^{n+1} – f^n = \Delta t \cdot F^n \qquad Euler 陽解法 \tag{1} \label{eq:scale-factor-1} \\\\
&f^{n+1} – f^n = \Delta t \cdot F^{n+1} \qquad Euler 陰解法 \tag{2} \label{eq:scale-factor-2} \\\\
&f^{n+1} – f^n = \Delta t \frac{F^n + F^{n+1}}{2} \qquad Cronk-Nicolson法, 2次AM法 \tag{3} \label{eq:scale-factor-3} \\\\
\end{align}
$$

ただ、空間と同じようにもちろん時間項の高次近似が存在します。それぞれ、

陽解法 → Adoms-Bashforth法 (1次~n次)
陽解法 → Runge-Kutta法 (1次~n次)
陰解法 → Adoms-Molton法 (1次~n次)

といいます。Euler 陽解法は1次のAdoms-Bashforth法であり、Euler 陰解法は1次のAdoms-Molton法であります。これらは、本当に単純な高次近似であり、空間の高次近似とあまり変わりません。Runge-Kutta法は面白くて、陽解法でありながらちょっと陰解法っぽくもあり精度が高いため、陽解法の時よく使われてます。

どれも一個ずつ導出してみます。

 

まず、最初に、関数を

$$
\begin{align}
&\frac{d u}{d t} =  f \tag{4} \label{eq:scale-factor-4}
\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{5} \label{eq:scale-factor-5}
\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 \frac{u_{i-1}-2u_i+u_{i+1}}{\Delta x^2} + f_i \tag{6} \label{eq:scale-factor-6}
\end{align}
$$

このように、空間を離散化した段階では、u_i = u_i (t) となり、xの依存はなくなり、偏微分方程式から、常微分方程式になります。このように、流体の支配方程式を空間離散化した後は、時間の常微分方程式になります。

つまりは、\(f(t,x_1,x_2,x_3)\)を、f(t)に直しただけのことで時間の偏微分と空間を離散化して時間の常微分にすることは全く同じなのですが、合成関数の微分と、偏微分と常微分がごっちゃになって分からなくなる人もいるのでこのようにします。


Adoms-Bashforth法

まず、\eqref{eq:scale-factor-4}式を、例のごとくTaylor展開します。ここで、まずuを時間tについてn時点周りに展開します。

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

ここで、uについての時間微分項が欲しいので、\(u^n\)を移行して、全体を\(\Delta t\)で割ります。

$$
\begin{align}
\frac{u^{n+1} – u^n}{\Delta t} &= \frac{1}{1!} \left.\frac{d u}{d t}\right\vert^n +\frac{\Delta t}{2!} \left.\frac{d^2 u}{d t^2}\right\vert^n +\frac{\Delta t^2}{3!} \left.\frac{d^3 u}{dt^3}\right\vert^n + \cdots  \tag{8} \label{eq:scale-factor-8}
\end{align}
$$

ここで、元の\eqref{eq:scale-factor-4}式より、\(\frac{d u}{d t}=f^n\)となるので、置き換えます。

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

また、今度はfについての後退差分をとります。これについては大丈夫だと思うので、まとめて書きます。

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

この2式で時間についての差分式を作っていきます。

ここで、2次精度とすると、作る式の形は、

$$
\begin{align}
&\frac{u^{n+1} – u^n}{\Delta t} + \beta_0 f^n + \beta_{-1} f^{n-1} =0 \tag{11} \label{eq:scale-factor-11}
\end{align}
$$

となります。各係数を表にしてみます。

 \(f^n\)\(\frac{d f}{d t}\)\(\frac{d^2 f}{d t^2}\)\(\frac{d^3 f}{d x^3}\)
\(\frac{u^{n+1}-u^n}{\Delta t}\)1\(\frac{1}{2}\)\(\frac{1}{6}\)\(\frac{1}{24}\)
\(f^n\)1000
\(f^{n-1}\)1-1\(\frac{1}{2}\)\(-\frac{1}{6}\)

あとは、各係数のつじつまを合わせていきます。それぞれの式を1倍、\(\beta_0\)倍、\(\beta_{-1}\)倍したとして式を作ります。0次( \(f^n\) )と1次の微分項を消せば良いです。\(f^n\)については、\(\beta_0)\)を掛けた方だけ残したいので、縦にみた時に消す形になります。表の左側が左辺、上側が右辺と考えてください。

\( \begin{eqnarray} \left\{ \begin{array}{l} 1 + \beta_0 + \beta_{-1} = 0 \\ \frac{1}{2} + 0 \times \beta_0 – \beta_{-1}= 0  \end{array} \qquad \Leftrightarrow \qquad \right.\left\{ \begin{array}{l} \beta_0 =-\frac{3}{2} \\ \beta_{-1} = \frac{1}{2} \end{array} \right. \end{eqnarray} \)

したがって、

$$
\begin{align}
& \frac{u^{n+1} – u^n}{\Delta t} – \frac{3}{2}f^{n} + \frac{1}{2} f^{n-1} = 0 \times f^n +  O(\Delta t^2) \\\\
&\Leftrightarrow \frac{u^{n+1} – u^n}{\Delta t} = \frac{1}{2} ( 3f^{n} – f^{n-1} ) + O(\Delta t^2)) \tag{12} \label{eq:scale-factor-12}
\end{align}
$$

となります。空間差分の時と似ているけど、微妙に違いますね。

そのまま、3次精度も一応やってみます。

 \(f^n\)\(\frac{d f}{d t}\)\(\frac{d^2 f}{d t^2}\)\(\frac{d^3 f}{d x^3}\)
\(\frac{u^{n+1}-u^n}{\Delta t}\)1\(\frac{1}{2}\)\(\frac{1}{6}\)\(\frac{1}{24}\)
\(f^n\)1000
\(f^{n-1}\)1-1\(\frac{1}{2}\)\(-\frac{1}{6}\)
\(f^{n-2}\)1-22(- \frac{8}{6}\)

\( \begin{eqnarray} \left\{ \begin{array}{l} 1 + \beta_0 + \beta_{-1} + \beta_{-2}= 0 \\ \frac{1}{2}  – \beta_{-1} – 2\beta{-2}= 0 \\ \frac{1}{6}  + \frac{1}{2}\beta_{-1} + 2\beta{-2}= 0 \end{array} \qquad \Leftrightarrow \qquad \right.\left\{ \begin{array}{l} \beta_0 =-\frac{23}{12} \\ \beta_{-1} = \frac{16}{12} \\ \beta_{-2} = – \frac{5}{12} \end{array} \right. \end{eqnarray} \)

したがって、

$$
\begin{align}
& \frac{u^{n+1} – u^n}{\Delta t} – \frac{23}{12}f^{n} + \frac{16}{12} f^{n-1} – \frac{5}{12} f^{n-2} =  O(\Delta t^3) \\\\
&\Leftrightarrow \frac{u^{n+1} – u^n}{\Delta t} = \frac{1}{12} ( 23f^{n} – 16f^{n-1} + 5f^{n-2} ) + O(\Delta t^3)) \tag{13} \label{eq:scale-factor-13}
\end{align}
$$

ちなみに、誤差項の具体的な値を計算したければ、式にしたがって足していけば出ます。
4次以降は省略します。4次と5次の形だけ書いときます。

$$
\begin{align}
&\frac{u^{n+1} – u^n}{\Delta t} = \frac{1}{24} ( 55f^{n} – 59f^{n-1} + 37f^{n-2} – 9f^{n-3} ) + O(\Delta t^4)) \tag{14} \label{eq:scale-factor-14} \\\\
&\frac{u^{n+1} – u^n}{\Delta t} = \frac{1}{720} ( 1901f^{n} – 2774f^{n-1} + 2616f^{n-2} -1274f^{n-3} + 251f^{n-4}) + O(\Delta t^5)) \tag{15} \label{eq:scale-factor-15} \\\\
\end{align}
$$
時間の場合は、2ステップ先の値など知れる方法がないので、空間の時と違って後ろにどんどん項を増やしていく形になります。これでやっと、一つ終わりました。

 


Adoms-Molton法

Adoms-Molton法は、 先ほどのAdoms-Bashforth法を陰解法にしただけなので、サクサクいきます。\eqref{eq:scale-factor-9}式と、\eqref{eq:scale-factor-10}式は同じで、fについての前進差分の式が追加されるだけです。

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

ここからもbashforthと同じです。係数を調整していきます。2次と3次まとめてやってみます。

求める形としては、2次は

$$
\begin{align}
&\frac{u^{n+1} – u^n}{\Delta t} + \gamma_1 f^{n+1} + \gamma_{0} f^{n} =0 \tag{17} \label{eq:scale-factor-17}
\end{align}
$$

となります。3次は、

$$
\begin{align}
&\frac{u^{n+1} – u^n}{\Delta t} + \gamma_1 f^{n+1} + \gamma_{0} f^{n} + \gamma_{-1} f^{n-1} =0 \tag{18} \label{eq:scale-factor-18}
\end{align}
$$

例のごとく係数を調整していきます。

 \(f^n\)\(\frac{d f}{d t}\)\(\frac{d^2 f}{d t^2}\)\(\frac{d^3 f}{d x^3}\)
\(\frac{u^{n+1}-u^n}{\Delta t}\)1\(\frac{1}{2}\)\(\frac{1}{6}\)\(\frac{1}{24}\)
\(f^{n+1}\)11\(\frac{1}{2}\)\(\frac{1}{6}\)
\(f^n\)1000
\(f^{n-1}\)1-1\(\frac{1}{2}\)\(-\frac{1}{6}\)

2次の場合は、

\( \begin{eqnarray} \left\{ \begin{array}{l} 1 + \gamma_1 + \gamma_{0} = 0 \\ \frac{1}{2} + \gamma_{1}= 0  \end{array} \qquad \Leftrightarrow \qquad \right.\left\{ \begin{array}{l} \gamma_1 = -\frac{1}{2} \\ \gamma_{0} = -\frac{1}{2} \end{array} \right. \end{eqnarray} \)

したがって、

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

3次の場合は、

\( \begin{eqnarray} \left\{ \begin{array}{l} 1 + \gamma_1 + \gamma_{0} + \gamma_{-1}= 0 \\ \frac{1}{2}  + \gamma_1 – \gamma_{-1} = 0 \\ \frac{1}{6}  + \frac{1}{2}\gamma_{1} + \frac{1}{2}\gamma_{-1} = 0 \end{array} \qquad \Leftrightarrow \qquad \right.\left\{ \begin{array}{l} \gamma_1 =-\frac{5}{12} \\ \gamma_{0} = -\frac{8}{12} \\ \gamma_{-1} = \frac{1}{12} \end{array} \right. \end{eqnarray} \)

したがって、

$$
\begin{align}
&\frac{u^{n+1} – u^n}{\Delta t} = \frac{1}{12} ( 5f^{n+1} + 8f^{n} – 1f^{n-1} ) + O(\Delta t^3)) \tag{20} \label{eq:scale-factor-20}
\end{align}
$$

となります。ついでに4次と5次も形だけ書いときます。

$$
\begin{align}
&\frac{u^{n+1} – u^n}{\Delta t} = \frac{1}{24} ( 9f^{n+1} + 19f^{n} – 5f^{n-1} + 1f^{n-2} ) + O(\Delta t^4)) \\\\
&\frac{u^{n+1} – u^n}{\Delta t} = \frac{1}{720} ( 251f^{n+1} + 646f^{n} – 264f^{n-1} + 100f^{n-2} -19f^{n-3} ) + O(\Delta t^5))
\tag{21} \label{eq:scale-factor-21}
\end{align}
$$

Adoms-Moltonの方も、陰解法と言っても、2つ以上先の時間は無理なので、後ろにどんどん項を追加していく形になります。

 


まとめ

Runge-Kutta法はとても長くなるので次回に回します。しかし、Runge-Kutta法はとても素晴らしいアルゴリズムです。

Adoms-Bashforth法も、AdomsMolton法も2次までしかあまり使われることがありません。理由は、高次になると安定領域が狭くなるからです。どういうことかというと、各項に対する発散しない\(\Delta t\)の幅を出すことができますが、Adoms-Bashforth法とAdoms-Molton法は高次になるとこの領域が狭くなります。つまり、高次精度にすると時間ステップをより狭くしなければならないということです。もちろん、精度を高めるために使っているので、時間ステップの幅も狭くするのは当然なのでまあいいのですが、高次精度にしたために発散しやすくなるのは少し厄介です。

しかし、Runge-Kutta法は高次にするほど安定領域が広がり、また計算量も4次までは精度と比例するので、4次のRunge-Kuttaは頻繁に使われます。

今回の記事は、名古屋工業大学さんの流体伝熱基礎工学第4回(森西洋平)を参考にさせていただきました。

また安定領域についてはとても詳しく解説してくれています。
めちゃくちゃいい記事なので、ぜひ読んでみてください。

流体伝熱基礎工学第4回(森西洋平)スライド
流体伝熱基礎工学第4回(森西洋平)資料

流体伝熱基礎工学

陰解法の方が精度が高いので、メッシュ数とかdtとか色々いじってみると面白いです。



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


<

p style=”text-align: center;”>