これまでも陰解法の議論(
もっとも単純な陰解法がImplicit Euler (Backward Euler)法である.オイラー陽解法を思い出してみると,以下の積分された常微分方程式の右辺第2項の積分区間において,
Euler陰解法の場合は,積分区間内で
当然のことながら,右辺に未知数を含んでいるため,このままでは解くことができない.
もし,離散化した後の式で,
右辺は全て既知量であるため,簡単に解くことができる.
ついでにオイラー陰解法の安定条件を調べてみる.安定性関数
オイラー陽解法とは異なり,安定領域が円の外側であることがわかる.この領域に対して,
さらにオイラー陰解法の精度を確認する.いつも通り,
次にオイラー陰解法は次のようになる.
これらの比較から区分の誤差(local error)は
実際にテスト方程式を解くプログラムを実行してみる.
program implicit_euler_method implicit none real*8 :: x=0.d0, y=1.0d0, x_end=1.0, lambda=-0.5, h=4.2 ! unstable case for explicit euler R = -1.5 integer :: n, n_max n_max = 120 write(*,*) x, y, exp(lambda*x) ! print initial condition do n=1, n_max y = y/(1-h*lambda) x = x+h write(*,*) x, y, exp(lambda*x) end do stop end program implicit_euler_method
実行結果を以下に示す.オイラー陽解法のように発散をせず,テスト方程式であっても安定していることがわかる.ただし,精度は1次精度であるため,精度を高めるためには,分割数をかなり大きくしなければならない.ただし,安定性を最優先するような計算では利用価値はある.
上述と異なり,代数的に解けない場合の方法として,後の節で説明する逐次近似法を利用する方法がある.ある初期値
さらに,また後の節で説明する非線形方程式の近似解法の一つNewton法を用いる方法もある.まず,方程式を整理して次の形にする.
ニュートン法に基づいて次のようにステップを繰り返して
ただし,この場合は
Heun法は予測子・修正子法であり,陰的な計算を行わなかった.それに対して,予測子を利用せず直接
クランク・ニコルソン法(
ここで
オイラー陰解法のように逐次近似法やニュートン法でも解けるが,代数的方法でテスト方程式を解くプログラムを作成して結果を示す.テスト方程式を含むプログラムは以下のようになる.
program crank_nicolson_method implicit none real*8 :: x=0.d0, y=1.0d0, x_end=1.0, lambda=-0.5, h=3 ! unstable case for explicit euler R = -1.5 integer :: n, n_max n_max = 120 write(*,*) x, y, exp(lambda*x) ! print initial condition do n=1, n_max y = y*(1+0.5*h*lambda)/(1-0.5*h*lambda) x = x+h write(*,*) x, y, exp(lambda*x) end do stop end program crank_nicolson_method
この方法は台形公式であるので,オイラー陽解法のような発散もなく,2次精度であるためオイラー陰解法よりも精度が高いことがわかる.