📂数値解析

ミッドポイントメソッド

メソッド

$D \subset \mathbb{R}^2$で定義された連続関数に対し、初期値問題$\begin{cases} y ' = f(x,y) \\ ( y( x_{0} ), y (x_{1}) ) = ( Y_{0} , Y_{1} ) \end{cases}$が$f$に与えられている。区間$(a,b)$を$a \le x_{0} < x_{1} < \cdots < x_{n} < \cdots x_{N} \le b$のようなノードポイントで分けたとする。特に十分に小さい$h > 0$に対し$x_{j} = x_{0} + j h$とすると、初期値$y_{0} = Y_{0}$に対して $$y_{n+1} := y_{n-1} + 2 h f ( x_{n} , y_{n} )$$

導出 ¹

$Y’(t)$を$x_{n}$に対してテイラー展開すると $$Y’(t) = y ' ( x_{n} ) + ( t - x_{n} ) Y’’ (x_{n} ) + {{( t - x_{n} )^2} \over {2}} Y^{(3)} ( \xi_{n} )$$ を満たす$\xi_{n} \in [ x_{n-1} , x_{n} ]$が存在する。両辺に$t= x_{n-1}$から$x_{n+1}$までの定積分を取ると $$\int_{x_{n-1} }^{ x_{n+1} } Y’(t) dt = \int_{x_{n-1} }^{ x_{n+1} } \left[ y ' ( x_{n} ) + ( t - x_{n} ) Y’’ (x_{n} ) + {{( t - x_{n} )^2} \over {2}} Y^{(3)} ( \xi_{n} ) \right] dt$$ $x_{n-1} = x_{n} - h$かつ$x_{n+1} = x_{n} + h$なので $$\int_{x_{n-1} }^{ x_{n+1} } Y’(t) dt = ( x_{n} + h - (x_{n} - h) ) y ' ( x_{n} ) + {{ h^3 - ( - h^3 ) } \over {6}} Y^{(3)} ( \xi_{n} )$$ 微積分学の基本定理により $$Y(x_{n+1}) - Y(x_{n-1})= 2h y ' ( x_{n} ) + {{ h^3 } \over {3}} Y^{(3)} ( \xi_{n} )$$ $Y’ = f$と仮定し、$\displaystyle {{ h^3 } \over {3}} Y^{(3)} ( \xi_{n} )$を除外すると $$y_{n+1} = y_{n-1} + 2 h f ( x_{n} , y_{n} )$$

■

オイラーメソッドとの比較

誤差分析

オイラーメソッドと比較すると、2個のデータを使う分、計算は増えるが誤差は小さくなると推測できる。オイラーメソッドの誤差は、多くの条件が与えられても$O (h)$にバウンドされているだけだ。

一方、ミッドポイントメソッドはマルチポイントメソッドとして、その誤差を計算すると $$\max_{ x_{0} \le x_{n} \le b} | Y (x_{n} ) - y_{h} (x_{n} ) | \le e^{2K(b-x_{0})} \eta (h) + \left[ {{ e^{2K (b - x_{0} ) } - 1 } \over {2K}} \right] \left[ {{1} \over {3}} h^2 | Y^{(3)} |_{\infty} \right]$$ ここで初期誤差は$\displaystyle \eta (h) = \max \left\{ |Y_{0} - y_{0}| , | Y(x_{1}) - y_{h} (x_{1}) | \right\}$で、$y_{0} =Y_{0}$だから実際には$\eta (h) = | Y(x_{1}) - y_{h} (x_{1}) |$になる。$Y$をテイラー展開すると $$Y (x_{1}) = Y(x_{0}) + h Y(x_{0}) + {{h^{2} } \over {2}} Y’’ (\zeta) = y_{0} + h f(x_{0} , y_{0}) + {{h^{2} } \over {2}} Y’’ (\zeta)$$ この時、オイラーメソッドによると$y_{0} + h f(x_{0} , y_{0}) = y_{1}$であるから、 $$Y(x_{1}) - y_{h} (x_{1}) = {{h^{2} } \over {2}} Y’’ (\zeta) = O(h^{2} )$$ 結果としてミッドポイントメソッドはオイラーメソッドよりも精度面で優れていると言える。$O (h^2)$にバウンドされているためだ。

安定性の問題

しかし、ミッドポイントメソッドには致命的な欠陥がある。以下はRでミッドポイントメソッドを実装したコードだ。

Midpoint<-function(f,Y\_0,Y\_1=NA,a,b,h=10^(-3))
{
  Y <- t(Y\_0)
  if(is.na(Y\_1))
  {
    Y<-rbind(Y,Y[1,]+h*f(a,Y[1,]))
  }
  else
  {
    Y<-rbind(Y,Y\_1)
  }
  node <- seq(a,b,by=h)
  
  for(x in node[-1])
  {
    Y<-rbind(Y,Y[length(Y[,1])-1,]+2*h*f(x,Y[length(Y[,1]),]))
  }
  
  return(Y)
}
 
f<-function(x,y) {1/(1+x^2) + - 2*(y^(2))}
out<-Midpoint(f,seq(0,0,len=1),0,2,h=0.1)
out[,1]
 
g<-function(x,y) {x-y^2}
out<-Midpoint(g,seq(0,0,len=1),0,3.25,h=0.25)
out[,1]

実行してみると、f()と違い、g()は適切に解が得られない。これはミッドポイントメソッドの安定性に問題があるという例を示している。