Excel VBA プログラミング

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

Download Report