数値解析におけるオイラー法
メソッド 1
$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]
Atkinson. (1989). 数値解析入門(第2版): p341. ↩︎