餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

Rで常微分方程式

自発的には揺れないボロビルの力学的振動は、単純化すると以下のように表されます。
\frac{d^2y}{d^2t} + k\frac{dy}{dt} + \omega^2 y = 0
yが揺れ幅(変位)、\omega固有振動数で、kが抵抗です。tは時間を表します。これをプロットしてみましょう。

1. 閉じた解を求める

定係数二階同次方程式なので、数学的に解くことができます。細かい計算は省略します。
y = C_1 e^{-\frac{k}{2} t} \cos{ \Bigg ( t \sqrt{ \omega^2 - \frac{k^2}{4}} \Bigg )} + C_2 e^{-\frac{k}{2}t} \sin{ \Bigg ( t \sqrt{ \omega^2 - \frac{k^2}{4}} \Bigg ) }
ほぼこのまま、curveで描画できます。

C1 <- 5; C2 <- 5; k <- 0.1; omega <- 1
curve(C1*exp(-k/2*t)*cos(t*sqrt(omega^2-k^2/4)) + C2*exp(-k/2*t)*sin(t*sqrt(omega^2-k^2/4)), 0, 100, n=1000, xname="t", ylab="y")

もっと複雑な式になると、辛い思いをします。

2. 微分方程式のまま処理する

y=y_1dy_1/dt=y_2と置いて、連立微分方程式に展開してみましょう。
\frac{dy_1}{dt} = y_2
\frac{dy_2}{dt} = - ky_2 -\omega^2 y_1
この状態で、deSolveパッケージに渡すことができます。

library("deSolve")

model <- function(t, Y, params){
# params["omega"]としなくていいようにwithを使う
  with(as.list(params),{
    dy1 = Y[2]
    dy2 = -omega^2*Y[1] - k*Y[2]
    list(c(dy1, dy2))
  })
}

# 変位Y1とその微分の初期値
Y <- c(y1=8, y2=0)
# t=0からt=100まで0.1刻みで計算する
times <- seq(0, 100, 0.1)
# 常微分方程式を計算する
out1 <- ode(Y, times, model, parms = c(omega=1, k=0.1, F0=0.35, beta=1))

# プロットする上下の限界
ylim <- c(-10, 10)

# プロットする。軸は書かない
plot(y1 ~ time, data=out1, type="l", ylim=ylim, axes=FALSE, ylab="揺れ幅(=変位)", xlab="時間", lwd=1, main="固有振動数と揺れ幅")

# 軸を描画
axis(1, labels=seq(0, 100, 25), at=seq(0, 100, 25), pos=ylim[1])
axis(2, labels=seq(-10, 10, 5), at=seq(-10, 10, 5), pos=0)

プロット結果は以下のようになります。

3. deSolveパッケージに関して

初期値問題、微分代数方程式、偏微分方程式も解ける優れモノです。詳しい説明は以下を参照してください。

  • Soetaert, Petzoldt, and Setzer (2010) "Solving Differential Equations in R: Package deSolve," Journal of Statistical Software, Vol.33(9)