コンピュータプログラミング
- 46 -
10. 高階常微分方程式や連立微分方程式の解法
10.1. 連立微分方程式の解法
一般的な形で書くのは大変でありわかりづらいので，ここでは２変数に限って議論する．以下 のような連立常微分方程式を考える．
この方程式は以下のようにルンゲ・クッタ法を用いて計算することができる．
ここで，
10.2. 高階常微分方程式を計算する
もとの高階常微分方程式を連立微分方程式にすることによって，高階常微分方程式の計算は行 うことができる．例えば以下のような2階常微分方程式
は，以下のように連立微分方程式に変更することができる．
これは と についての連立微分方程式なので，上に示した方法で数値的に解くことができる．
10.3. 例（減衰振動）
簡単だけど面白い例として減衰振動の方程式を考える．
この方程式は解析的に解くことができる．解は判別式 の値によって区別される．
のとき指数関数で表される減衰を示し， のときは振動しながら減衰する．この様子を
ルンゲ・クッタ法によって計算してみよう．
コンピュータプログラミング
- 47 -
パラメータ eta を変更して幾つかの解を得てグラフにまとめてみよう．
Sub calcDampedOsc()
Const initialX As Double = 1#
Const initialV As Double = 0#
Const deltaT As Double = 0.2, maxT As Double = 20 Dim x As Double, newX As Double, t As Double Dim v As Double, newV As Double
Dim kx1 As Double, kx2 As Double, kx3 As Double, kx4 As Double Dim kv1 As Double, kv2 As Double, kv3 As Double, kv4 As Double Dim counter As Integer
ActiveCell.Value = "Time"
ActiveCell.Offset(0, 1).Value = "x(t)"
ActiveCell.Offset(0, 2).Value = "v(t)"
ActiveCell.Offset(1, 0).Value = 0
ActiveCell.Offset(1, 1).Value = initialX ActiveCell.Offset(1, 2).Value = initialV x = initialX
v = initialV counter = 2
For t = deltaT To maxT Step deltaT kx1 = deltaT * doX(x, v)
kv1 = deltaT * doV(x, v)
kx2 = deltaT * doX(x + kx1 / 2, v + kv1 / 2) kv2 = deltaT * doV(x + kx1 / 2, v + kv1 / 2) kx3 = deltaT * doX(x + kx2 / 2, v + kv2 / 2) kv3 = deltaT * doV(x + kx2 / 2, v + kv2 / 2) kx4 = deltaT * doX(x + kx3, v + kv3)
kv4 = deltaT * doV(x + kx3, v + kv3)
newX = x + (kx1 + 2 * kx2 + 2 * kx3 + kx4) / 6 newV = v + (kv1 + 2 * kv2 + 2 * kv3 + kv4) / 6
ActiveCell.Offset(counter).Value = t ActiveCell.Offset(counter, 1).Value = newX ActiveCell.Offset(counter, 2).Value = newV x = newX
v = newV
counter = counter + 1 Next t
End Sub
Function doX(x As Double, v As Double) As Double doX = v
End Function
Function doV(x As Double, v As Double) As Double Const k As Double = 0.5
Const eta As Double = 2.5
doV = -eta * v - k * x End Function
コンピュータプログラミング
- 48 - 計算例を下に示した．
10.4. 練習問題
A) 上の減衰振動において，振動がなくなるギリギリの値を見つけなさい．
B) 他の微分方程式についても計算する方法を考えなさい．例えば以下の方程式ではどうなる だろうか？
1) （単振動）
2) （van der Pol 方程式）
3) （Lotka-Volterra 方程式）
