Rで常微分方程式
自発的には揺れないボロビルの力学的振動は、単純化すると以下のように表されます。
が揺れ幅(変位)、が固有振動数で、が抵抗です。は時間を表します。これをプロットしてみましょう。
1. 閉じた解を求める
定係数二階同次方程式なので、数学的に解くことができます。細かい計算は省略します。
ほぼこのまま、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. 微分方程式のまま処理する
、と置いて、連立微分方程式に展開してみましょう。
この状態で、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)