📂数値解析

数値解析におけるオイラー法

メソッド ¹

$D \subset \mathbb{R}^2$ で定義された連続関数について、初期値問題が$\begin{cases} y ' = f(x,y) \\ y( x_{0} ) = Y_{0} \end{cases}$で与えられている。区間$(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} \simeq Y_{0}$に対して $$y_{n+1} = y_{n} + h f ( x_{n} , y_{n} )$$

説明

オイラーメソッドは概念的に非常にシンプルな方法だが、数値解析の核心的なアイデアを示している。もちろん今では様々な問題が多いが、逆に言えば改善の余地も多い手法である。だからといって、オイラーメソッド自体が重要というわけではなく、導出過程をしっかりと理解することが必要だ。

導出

テイラー展開

$Y(x_{n+1} )$を$x_{n}$に対して$2$項までテイラー展開すると $$Y ( x_{n+1} ) = Y ( x_{n} ) + ( x_{n+1} - x_{n}) y ' ( x_{n} ) + {{( x_{n+1} - x_{n})^2} \over {2}} Y’’ ( \xi_{n} )$$ 整理すると $$Y ( x_{n+1} ) = Y ( x_{n} ) + h y ' ( x_{n} ) + {{h^2} \over {2}} Y’’ ( \xi_{n} )$$ ここで、誤差項$\displaystyle {{h^2} \over {2}} Y’’ ( \xi_{n} )$を無視すると $$y_{n+1} = y_{n} + h f ( x_{n} , y_{n} )$$

■

テイラー展開を用いた導出から直ちに分かることは、$h$が小さくなると誤差も減少するということだ。また、テイラー近似がどうして$2$次でなければならないという特別な理由はないため、次数を上げれば誤差が減少するということだ。

数値積分

$[x_{n} , x_{n+1}]$から$Y’(t) = f (t, Y(t))$の定積分は $$\int_{x_{n}}^{x_{n+1}} f(t,Y(t)) dt = Y(x_{n+1}) - Y(x_{n})$$ 少しの誤差はあるが、数値積分によって$\displaystyle \int_{x_{n}}^{x_{n+1}} f(t,Y(t)) dt \simeq h f(x_{n} , Y(x_{n}) )$であるため、 $$y_{n+1} - y_{n} = h f ( x_{n} , y_{n} )$$

■

分割積分の考え方からして、$h$が小さくなればなるほど誤差も小さくなると推測できる。図からわかるように、実際の積分値と数値的に計算された値の差はかなり大きいが、台形公式やシンプソン法を使えば誤差が減少するだろう。

実装

以下はRで書かれたコードだ。

Euler<-function(f,Y\_0,a,b,h=10^(-3))
{
  Y <- t(Y\_0)
  node <- seq(a,b,by=h)
  
  for(x in node)
  {
    Y<-rbind(Y,Y[length(Y[,1]),]+h*f(x,Y[length(Y[,1]),]))
  }
  
  return(Y)
}
 
f<-function(x,y) {y}
out<-Euler(f,seq(1,1,len=100),0,2)
out[,1]
 
g<-function(x,y) {1/(1+x^2) + - 2*(y^(2))}
out<-Euler(g,seq(0,0,len=100),0,2,h=0.2)
out[,1]