陰解法

はじめに

これまでも陰解法の議論(yi+1を含む定式化)はなされてきたが,実際はHeun法,Adams-Bashforth-Moulton法も,yi+1を予測子で近似して,陰的なスキームを陽的なスキームに変換して計算を進めた.陰解法とはそのような近似を行わず直接解く方法を指す.

オイラー陰解法 (Implicit Euler or Backward-Euler method)

もっとも単純な陰解法がImplicit Euler (Backward Euler)法である.オイラー陽解法を思い出してみると,以下の積分された常微分方程式の右辺第2項の積分区間において,f(x,y)f(xi,yi)で一定とした.

(1)y(xi+1)=y(xi)+xixi+1f(x,y)dxy(xi)+f(xi,yi)h

Euler陰解法の場合は,積分区間内でf(x,y)f(xi+1,yi+1)で一定だと考える.つまり以下のようになる.

(2)y(xi+1)=y(xi)+xixi+1f(x,y)dxy(xi)+f(xi+1,yi+1)h

当然のことながら,右辺に未知数を含んでいるため,このままでは解くことができない.

代数的方法(簡単に解ける場合)

もし,離散化した後の式で,yi+1で括って代数的に解ける場合には簡単に解くことができる.例えば,安定性の節で議論したテスト方程式,dydx=λyの場合であれば以下のように変形する.

yi+1=yi+hλyi+1(3)yi+1=yi1hλ

右辺は全て既知量であるため,簡単に解くことができる.

安定性について

ついでにオイラー陰解法の安定条件を調べてみる.安定性関数|R(z)|=|11z|1を満足する必要がある.安定性関数R(z)は(両辺に|1z|をかけることで)以下のように変形でき,他のスキームと同様な安定領域を図かすることができる.

(4)|1z|={1Re(z)}2+{Im(z)}21

オイラー陽解法とは異なり,安定領域が円の外側であることがわかる.この領域に対して,R(z)=|11z|をプロットすると次のようになる.

精度について

さらにオイラー陰解法の精度を確認する.いつも通り,yi+1のテイラー展開との比較を行う.まずyi+1のテイラー展開は以下のようになる.

(5)y(xi+1)=y(xi+h)=y(xi)+hdydx+h22!d2ydx2+=y(xi)+hf(xi,yi)+h22!f(xi,yi)+

次にオイラー陰解法は次のようになる.

(6)y(xi+1)=y(xi)+f(xi+1,yi+1)h=y(xi)+hf(xi,yi)+h2f(xi,yi)+

これらの比較から区分の誤差(local error)はO(h2)であり,全区間での誤差(global error)はO(h)となり,オイラー陽解法と同様に一次精度であることがわかる.

実際にテスト方程式を解くプログラムを実行してみる.

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次精度であるため,精度を高めるためには,分割数をかなり大きくしなければならない.ただし,安定性を最優先するような計算では利用価値はある.

(代数的に解けない場合:その1)逐次近似的な方法(収束計算)

上述と異なり,代数的に解けない場合の方法として,後の節で説明する逐次近似法を利用する方法がある.ある初期値yi+1(0)から開始して,以下の方程式を繰り返し計算してyi+1の収束解を求めるという方法がある.

(7)y(n+1)(xi+1)=y(xi)+f(xi+1,yi+1(n))h

(代数的に解けない場合:その2)ニュートン法を用いた計算

さらに,また後の節で説明する非線形方程式の近似解法の一つNewton法を用いる方法もある.まず,方程式を整理して次の形にする.

(8)F(yi+1)=yi+1f(xi+1,yi+1)hyi=0

ニュートン法に基づいて次のようにステップを繰り返してyi+1を更新していく.

(9)yi+1(n+1)=yi+1(n)yi+1(k)f(xi+1,yi+1(k))hyi1hdf(xi+1,yi+1(k))dyi+1

ただし,この場合はdfdyを自分で求めなければならない(テスト方程式などであれば簡単に求められる).

クランク・ニコルソン法 (台形公式,Crank-Nicolson method)

Heun法は予測子・修正子法であり,陰的な計算を行わなかった.それに対して,予測子を利用せず直接yi+1を利用する方法があり,クランク・ニコルソン法と呼ばれている.この方法は台形公式であり,fifi+1の重みはそれぞれ12である.一般に,重みを調整する変数α(01)を導入して以下のように表すことが多い.

(10)yi+1=yi+h{(1α)f(xi,yi)+αf(xi+1,yi+1)}

α=0.5の場合はクランク・ニコルソン法になり,α=0の場合はオイラー陽解法,α=1.0の場合はオイラー陰解法となる(この方法は記号にαの代わりにθを用いて,theta-methodとも呼ばれることがある).

安定性について

クランク・ニコルソン法(α=0.5)に関して安定性を確認する.クランク・ニコルソン法で離散化されたテスト方程式は以下のように変形することができる.

(11)yn+1=1+z21z2y(n)

ここでz=hλである.上の式から,安定性条件は|R(z)|=|1+z21z2|1ということである.これを解くのは少しややこしい(例えばMIT Math, Numerics, & Programming, pp.338など.確かに代入してみると成り立つのはわかる)ので,結果だけ示すとRe(z)0(虚軸が境界)となり,以下のような図の領域が安定領域となる.

オイラー陰解法のように逐次近似法やニュートン法でも解けるが,代数的方法でテスト方程式を解くプログラムを作成して結果を示す.テスト方程式を含むプログラムは以下のようになる.

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次精度であるためオイラー陰解法よりも精度が高いことがわかる.