読者です 読者をやめる 読者になる 読者になる

廿TT

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

オイラー法, テイラー級数法, ルンゲ=クッタ法(R による練習)

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
(点と点の間を線分でつないでいるが, 本当は2次曲線でつなぐのが適切.)

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:20150130033239p: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")))

オイラー法, テイラー級数法, ルンゲ=クッタ法それぞれの誤差についてはまた後に考える.