オイラー法(Euler's Method)
オイラー法(Euler's Method)とは, 1階常微分方程式の数値解法の中でおそらくもっともかんたんなもの.
この方法は、数学的に理解しやすく、プログラム的にも簡単なので、数値解析の初歩的な学習問題としてよく取りあげられる。しかし、1階常微分方程式の数値解法としては誤差が蓄積されるため、精度が悪く、元の微分方程式によってはいかなるhをとっても元の方程式の解に収束しないこともある方法なので、学習目的以外であまり使われない。
オイラー法 - Wikipedia
とのこと.
今回は例として, リッカティの方程式(リッカチの微分方程式 - Wikipedia)
を考える.
これは解析的には解けないので, オイラー法によって近似解を求める.
解は一意でない(積分したら任意の値をとる積分定数が出てくる)ので, 初期値 を指定しておく.
Rによる数値実験(オイラー法)
初期値を, と適当に選び, ステップの大きさを
と半分にしていく.
収束の様子を下に図示する.
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 のまわりでテイラー展開( 上式において,
と置き換えればよい)すると,
オイラー法がテイラー級数の最初の2項しか使わないことに着目し、3項使ってみると,
.
オイラー法と同様, 一連の値,
が得られる.
R による数値実験(テイラー級数法)
引き続きリッカティの方程式
を例にとる.
y は x の関数なので,
.
初期値を, と適当に選び, ステップの大きさを
と半分にしていく.
収束の様子を下に図示する. オイラー法より格段に収束が早い.
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)と呼ばれるものらしい
初期値,
に対して、
を与える.
ここで、
である. すなわち, 次の値( )は, 現在の値(
)に間隔( h)と推定された勾配の積を加えたものである. その勾配は次の4つの勾配の重み付け平均で求める.
R による数値実験(ルンゲ=クッタ法)
引き続き, リッカティの方程式 を例にとり, 初期値を,
と適当に選び, ステップの大きさを
と半分にしていく.
圧倒的に収束が早い.
しかも高階の導関数を求めなくてよい.
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")))
参考文献
- オイラー法 - Wikipedia
- 認証がかかっています
- 接線の方程式
- 修正オイラー法
- http://www-mmc.es.hokudai.ac.jp/else/cdrom/main_part1.pdf
- ルンゲ=クッタ法 - Wikipedia

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