多段階解法

多段解法No.3

予測子・修正子法（Predictor‐Corrector Method）

これまでも繰り返し述べたように，(常)微分方程式は厳密に次に示す積分された形に変形することができる．

\begin{matrix} (1) & y (x_{i + 1}) = y (x_{i}) + \int_{x_{i}}^{x_{i + 1}} f (x, y) d x \end{matrix}

右辺の積分を既知の値だけで表現する陽解法，未知数 $y_{i + 1}$ を含む形で表現する陰解法があることを説明してきた．そして，陰解法は陽解法と違い単純な代入だけで数値解を求めることができない（一般的に繰り返し計算が必要であったり，逆行列を求めなければならない）．一方で陰解法のメリットとして，Adams-Moulton法で学んだように $f_{i + 1}$ の値を含んで被積分関数を構築（近似）することで精度を向上させることができた．

この陽解法のメリットと陰解法のメリットを組み合わせる（ある意味の妥協としての）アイディアが，予測子・修正子法である．つまり，陽解法によって $y_{i + 1}$ の予測値， ${\hat{y}}_{i + 1}$ を構築し，この値を陰解法に代入して最終的な解を得る．この時利用する陽解法を「予測子」，一方の陰解法を「修正子」とよび，組み合わせた方法を予測子・修正子法と呼んでいる．

ここで，何か既にやったことがあるような・・・となると思う．この方法は既にHeun法で利用している（Heun法はもっとも低次の予測子・修正子法である）．ここでは，Adams-Bashforth系の方法とAdams-Moulton系の方法を組み合わせた予測子・修正子法，つまり，より高次のAdams-Bashforth-Moulton法を紹介する．と言っても組み合わせるだけであるので特に難しいことはない．

組み合わせは何通りもあるので，３次のAdams-Bashforth法と４次のAdams-Moulton法を組み合わせた方法を紹介する（同じステップ数なので，記憶しておく値が共通している）．

まず，３次のAdams-Bashforth法で $\hat{y} (x_{i + 1})$ を求める．

\begin{matrix} (2) & \hat{y} (x_{i + 1}) = y (x_{i}) + \frac{h}{12} (5 f_{i - 2} - 16 f_{i - 1} + 23 f_{i}) \end{matrix}

次に，この $\hat{y} (x_{i + 1})$ の値を用いて $f (x_{i + 1}, y_{i + 1})$ の予測値 $\hat{f} (x_{i + 1}, y_{i + 1}) = {\hat{f}}_{i + 1}$ を求める．最終的に，この予測値 ${\hat{f}}_{i + 1}$ を用いて４次のAdams-Moulton法で計算値を求める．

\begin{matrix} (3) & y (x_{i + 1}) = y (x_{i}) + \frac{h}{24} (f_{i - 2} - 5 f_{i - 1} + 19 f_{i} + 9 {\hat{f}}_{i + 1}) \end{matrix}

これらの手順に基づいて，分割数に対する誤差の変化を計算するFortranプログラムは以下のようになる．（多少雑であるが）３次の予測子・修正子法が利用できるまでの初期ステップにおいては，Heun法を利用している．

-->

 program error_pc
  implicit none
  real*8 :: x, y, y_p1, x_end=1.d0, h, f_i, f_m1, f_m2, f_p1
  integer :: i, n, n_max
  
  do i=1,6
    n_max = 2**i
    h     = x_end/n_max
    x=0.d0
    y=1.0d0
    do n=1, n_max
     f_i  = f(x,y)
     
     if (n < 3) then 
       y_p1 = y+f_i*h
       f_p1 = f(x+h, y_p1)
       y    = y+0.5*h*(f_i+f_p1) ! Heun method
     else
       y_p1 = y+h/12.d0*(5*f_m2-16.d0*f_m1+23.d0*f_i) ! 3rd-order AB
       f_p1 = f(x+h, y_p1)
       y    = y+h/24.d0*(f_m2-5.d0*f_m1+19.d0*f_i+9.d0*f_p1) ! 4th-order AM
     end if
     f_m2 = f_m1
     f_m1 = f_i
     x    = x+h
    end do
    write(*,*)n_max, abs(y-exp(1.d0))/exp(1.d0)
  end do
  stop
  
  contains
  
  function f(x,y)
    implicit none
    double precision, intent(in) :: x, y
    double precision :: f
    f = 3.d0*x**2*y
  end function f
end program error_pc

Gnuplotのスクリプト

 reset
set log xy
set title 'Error by Predictor-Corrector method'
set xlabel 'n'
set ylabel 'E_{total}'
set format y "%.2e"
set title  font "Arial,20"
set xlabel font "Arial,20"
set ylabel font "Arial,20"
set key    font "Arial,12"
plot 'error_pc.txt' with linespoints linewidth 3, 1/x, 1/x**2, 1/x**3
set terminal png
set output 'error_pc.png'
replot

３次の予測子・修正子法に関する数値誤差の分割数依存性を以下に示す．分割数が増大するにつれて３次精度の傾きに近くなっていることがわかる．

	program error_pc
	implicit none
	real*8 :: x, y, y_p1, x_end=1.d0, h, f_i, f_m1, f_m2, f_p1
	integer :: i, n, n_max

	do i=1,6
	n_max = 2**i
	h = x_end/n_max
	x=0.d0
	y=1.0d0
	do n=1, n_max
	f_i = f(x,y)

	if (n < 3) then
	y_p1 = y+f_i*h
	f_p1 = f(x+h, y_p1)
	y = y+0.5h(f_i+f_p1) ! Heun method
	else
	y_p1 = y+h/12.d0(5f_m2-16.d0f_m1+23.d0f_i) ! 3rd-order AB
	f_p1 = f(x+h, y_p1)
	y = y+h/24.d0(f_m2-5.d0f_m1+19.d0f_i+9.d0f_p1) ! 4th-order AM
	end if
	f_m2 = f_m1
	f_m1 = f_i
	x = x+h
	end do
	write(,)n_max, abs(y-exp(1.d0))/exp(1.d0)
	end do
	stop

	contains

	function f(x,y)
	implicit none
	double precision, intent(in) :: x, y
	double precision :: f
	f = 3.d0x2y
	end function f
	end program error_pc

	reset
	set log xy
	set title 'Error by Predictor-Corrector method'
	set xlabel 'n'
	set ylabel 'E_{total}'
	set format y "%.2e"
	set title font "Arial,20"
	set xlabel font "Arial,20"
	set ylabel font "Arial,20"
	set key font "Arial,12"
	plot 'error_pc.txt' with linespoints linewidth 3, 1/x, 1/x2, 1/x3
	set terminal png
	set output 'error_pc.png'
	replot