オイラー法, テイラー級数法, ルンゲ＝クッタ法（R による練習）

オイラー法（Euler's Method）

オイラー法（Euler's Method）とは, 1階常微分方程式の数値解法の中でおそらくもっともかんたんなもの.

この方法は、数学的に理解しやすく、プログラム的にも簡単なので、数値解析の初歩的な学習問題としてよく取りあげられる。しかし、1階常微分方程式の数値解法としては誤差が蓄積されるため、精度が悪く、元の微分方程式によってはいかなるhをとっても元の方程式の解に収束しないこともある方法なので、学習目的以外であまり使われない。
オイラー法 - Wikipedia

とのこと.

今回は例として, リッカティの方程式（リッカチの微分方程式 - Wikipedia）

$\displaystyle y'=f(x,y)= x^2 + y^2$

を考える.

これは解析的には解けないので, オイラー法によって近似解を求める.

解は一意でない（積分したら任意の値をとる積分定数が出てくる）ので, 初期値 $F(x_0) =y(x_0) = y_0$ を指定しておく.

オイラーのアイデア

$h>0$ を選んで $x_0 \le x \le x_0 + h$ では, 解を接線,

$\displaystyle l(x) =y_0 +(x-x_0) \cdot f(x_0, y_0)$

で置き換える.

点 $x_1 = x_0 + h$ では, $y_1 =y_0 + f (x_0, y_0)$ となる.

この点でまた, 同じことを繰り返し,一連の値,

$\displaystyle x_{n+1} = x_n + h, \, y_{n+1}= y_n +h f(x_n,y_n)$

が得られる.

これがオイラー法である.

Rによる数値実験（オイラー法）

初期値を, $x_0=-1.5, y_0=-1.4$ と適当に選び, ステップの大きさを $h=1/4, 1/8, 1/16, 1/32$ と半分にしていく.

収束の様子を下に図示する.

f:id:abrahamcow:20150130033850p:plain

testEuler<-function(h){
f <- function(x,y){x^2 + y^2}
#h <- 1/4
x_old <- -1.5
y_old <- -1.4 

x_new <- numeric()
y_new <- numeric()
i=1
while(x_old <= 1.5){
x_new[i] <- x_old + h
y_new[i] <- y_old + h*f(x_old, y_old)

x_old <- x_new[i]
y_old <- y_new[i]
i=i+1
}
return(data.frame(x=x_new,y=y_new))
}

res1 <- lapply(1/2^(2:5),testEuler)

plot(c(-1.5,1.5),c(-1.5,1.5),type="n",ylab="y",xlab="x",main="Euler Method")
points(res1[[4]],type="b", col="royalblue",pch=20)
points(res1[[3]],type="b", col="tomato",pch=20)
points(res1[[2]],type="b",col="orange2",pch=20)
points(res1[[1]],type="b",col="forestgreen",pch=20)


legend("topleft",rev(c("h=1/32","h=1/16","h=1/8","h=1/4")),pch=20,
       col=rev(c("royalblue","tomato","orange2","forestgreen")))

テイラー級数法

オイラー法に引き続き, テイラー級数法を考える.

テイラー級数法はオイラー法を改良した常微分方程式の数値解法で, 修正オイラー法とも言う.

テイラー級数法の説明

説明のために関数 f(x) を $x=x_0$ のまわりでテイラー展開した式を書いておく.

$\displaystyle f(x) = f(x_0) + \frac{d}{dx} f(x-x_0) + \frac{1}{2}\frac{d^2}{dx^2} f(x-x_0)^2 + O( (x-x_0)^3 )$

$y(x+h)$ を x のまわりでテイラー展開（上式において, $x=x+h, x_0=x$ と置き換えればよい）すると,

$\displaystyle y_{x+h} =y_{t} + \frac{d}{dt}y(t) h + \frac{1}{2} \frac{d^2}{dt^2}y(t) h^2 + O(h)$

オイラー法がテイラー級数の最初の2項しか使わないことに着目し、3項使ってみると,

$\displaystyle y_{n+1} = y_n + h y_n’ + \frac{h^2}{2}y_n’’$ .

オイラー法と同様, 一連の値,

$\displaystyle x_{n+1} = x_n + h, \, y_{n+1}= y_n +h y_n’ + \frac{h^2}{2}y_n’’$

が得られる.

R による数値実験（テイラー級数法）

引き続きリッカティの方程式

$\displaystyle y'=f(x,y)= x^2 + y^2$

を例にとる.

y は x の関数なので,

$\displaystyle y''=\frac{d}{d x}f(x,y)= 2x+2yy’ = 2x+2x^3y +2y^3$ .

初期値を, $x_0=-1.5, y_0=-1.4$ と適当に選び, ステップの大きさを $h=1/4, 1/8, 1/16, 1/32$ と半分にしていく.

収束の様子を下に図示する. オイラー法より格段に収束が早い.

f:id:abrahamcow:20150130040329p:plain

modEuler<-function(h){
  f <- function(x,y){x^2 + y^2}
  f2 <- function(x,y){2*x + (2*x^2)*y + 2*y^3}
  #h <- 1/4
  x_old <- -1.5
  y_old <- -1.4 
  
  x_new <- numeric()
  y_new <- numeric()
  i=1
  while(x_old <= 1.5){
    x_new[i] <- x_old + h
    y_new[i] <- y_old + h*f(x_old, y_old) + ((h^2)/2)*f2(x_old,y_old)
    
    x_old <- x_new[i]
    y_old <- y_new[i]
    i=i+1
  }
  return(data.frame(x=x_new,y=y_new))
}

res1 <- lapply(1/2^(2:5),modEuler)

plot(c(-1.5,1.5),c(-1.5,1.5),type="n",ylab="y",xlab="x", main="modified Euler method")
points(res1[[4]],type="b", col="royalblue",pch=20)
points(res1[[3]],type="b", col="tomato",pch=20)
points(res1[[2]],type="b",col="orange2",pch=20)
points(res1[[1]],type="b",col="forestgreen",pch=20)

legend("topleft",c("h=1/32","h=1/16","h=1/8","h=1/4"),pch=20,
       col=c("royalblue","tomato","orange2","forestgreen"))

ノート

もちろんテイラー級数の4項目, 5項目, ... , まで使うことも可能である.

ルンゲ＝クッタ法（Runge-Kutta method）

古典的な4次ルンゲ＝クッタ法

とりあえず Wikipedia（ルンゲ＝クッタ法 - Wikipedia）にある通りに実装して動かしてみる.

一般に用いられているルンゲ＝クッタ法は4次のルンゲ＝クッタ法（RK4）と呼ばれるものらしい

初期値,

$y' = f(t, y), \quad y(t_0) = y_0$

に対して、

$\displaystyle y_{n+1} = y_n + {h \over 6} (k_1 + 2k_2 + 2k_3 + k_4)$

を与える.

ここで、

$\displaystyle k_1 = f \left( t_n, y_n \right)$

$\displaystyle k_2 = f \left( t_n + {h \over 2}, y_n + {h \over 2} k_1 \right)$

$\displaystyle k_3 = f \left( t_n + {h \over 2}, y_n + {h \over 2} k_2 \right)$

$\displaystyle k_4 = f \left( t_n + h, y_n + hk_3 \right)$

である. すなわち, 次の値（ $y_{n+1}$ ）は, 現在の値（ $y_n$ ）に間隔（ h）と推定された勾配の積を加えたものである. その勾配は次の4つの勾配の重み付け平均で求める.

$k_1$ は初期値における勾配
$k_2$ は区間の中央における勾配
$k_3$ は区間の中央における勾配
$k_4$ は区間の最後における勾配

R による数値実験（ルンゲ＝クッタ法）

引き続き, リッカティの方程式 $y'=f(x,y)= x^2 + y^2$ を例にとり, 初期値を, $x_0=-1.5, y_0=-1.4$ と適当に選び, ステップの大きさを $h=1/4, 1/8, 1/16, 1/32$ と半分にしていく.

f:id:abrahamcow:20170510074214p:plain

圧倒的に収束が早い.

しかも高階の導関数を求めなくてよい.

testRG<-function(h){
  f <- function(x,y){x^2 + y^2}
  #h <- 1/4
  x_old <- -1.5
  y_old <- -1.4 
  
  x_new <- numeric()
  y_new <- numeric()
  i=1
  while(x_old <= 1.5){
    x_new[i] <- x_old + h
    k1 <- f(x_old, y_old)
    k2 <- f(x_old+h/2, y_old+(h/2)*k1)
    k3 <- f(x_old+h/2, y_old+(h/2)*k2)
    k4 <- f(x_old+h, y_old+(h/2)*k3)
    y_new[i] <- y_old + (h/6)*(k1+2*k2+2*k3+k4)
    
    x_old <- x_new[i]
    y_old <- y_new[i]
    i=i+1
  }
  return(data.frame(x=x_new,y=y_new))
}


res1 <- lapply(1/2^(2:5),testRG)
plot(c(-1.5,1.5),c(-1.5,1.5),type="n",ylab="y",xlab="x",main="Runge-Kutta method")
points(res1[[4]],type="b", col="royalblue",pch=20)
points(res1[[3]],type="b", col="tomato",pch=20)
points(res1[[2]],type="b",col="orange2",pch=20)
points(res1[[1]],type="b",col="forestgreen",pch=20)

legend("topleft",rev(c("h=1/16","h=1/8","h=1/4","h=1/2")),pch=20,
       col=rev(c("royalblue","tomato","orange2","forestgreen")))

参考文献

解析教程・上新装版

作者: 蟹江幸博
出版社/メーカー: 丸善出版
発売日: 2012/04/20
メディア: 単行本
クリック: 2回
この商品を含むブログを見る

（pp.183-187 を参照）

廿TT

譬如水怙牛過窓櫺　頭角四蹄都過了　因甚麼尾巴過不得