Hatena::ブログ(Diary)

My Life as a Mock Quant このページをアンテナに追加 RSSフィード Twitter

2015-05-31

常微分方程式をRで解く際の精度

| 12:26 | 常微分方程式をRで解く際の精度を含むブックマーク

 ¥frac{dy}{dx} = -¥frac{¥sin(x)}{y}, y(0) = 1

という常微分方程式の解は

 y(x) = ¥sqrt{2¥cos(x) -1}

で、これをRで数値的に確かめたい。

ついでに、数値積分手法(ここではEuler法とRunge-Kutta4次)によって、どのくらい厳密解とずれるのかも調べる。

Rで常微分方程式を解く際にはdeSolveパッケージを用いるのが楽。結果を算出するためのコードは以下。

library(deSolve)
#厳密解
exact <- function(t){sqrt(2*cos(t)-1)}
#deSolveパッケージのode関数に食わせるための常微分方程式
f <- function(t, y, p){list(-sin(t)/y)}
#始状態から順方向で+終状態から逆方向に+厳密解の結果を同時に計算する関数
solve <- function(method)
{
  t1 <- seq(0, 1, by = 0.05)
  t2 <- seq(1, 0, by = -0.05)
  rbind(
    data.frame(ode(c(y=exact(0)), t1, f, c(), method), type="F"),
    data.frame(ode(c(y=exact(1)), t2, f, c(), method), type="B"),
    data.frame(time=t1, y=exact(t1),type="E")
  )
}

結果をggplot2でEuler法 v.s. Runge-Kutta 4次で対比できるように描画すると・・・

library(ggplot2)
ggplot(solve("euler"), aes(x=time, y=y)) +
  geom_path(size=2, aes(color=type)) +
  ggtitle("Differential Equation Solve(Euler)")

ggplot(solve("rk4"), aes(x=time, y=y)) +
  geom_path(size=2, aes(color=type)) +
  ggtitle("Differential Equation Solve(Runge-Kutta 4th order)")

Euler法では

  • 始状態から前向きに解く(図中F)
  • 終状態から後向きに解く(図中B)
  • 厳密解(図中E)

が一致しないが、Runge-Kutta 4次だと同じ時間(x)刻みでほぼ答えが合う。

なお、始状態からの場合は厳密解の時点0の値と、終状態からの場合は厳密解の時点1の値を使って計算を開始している。

f:id:teramonagi:20150531122245p:image

f:id:teramonagi:20150531122246p:image

トラックバック - http://d.hatena.ne.jp/teramonagi/20150531/1433042797
リンク元