連立常微分方程式

ルンゲ・クッタ法の連立常微分方程式への適用

ルンゲクッタ法は開始点の値（初期値）がわかれば多変数の連立常微分方程式にも適用できる。

逐次反応（「反応工学」より）

最も簡単な逐次反応は

A→R,

R→S,

である。各成分の濃度は以下の連立常微分方程式となる。

ただし、初期濃度 C_A=1, C_R=0, C_S=0 とする。

先ず解析的に解いてみる。はじめのの解は　すなわち。これを第２式に代入すると、

を乗ずると　であり、この式は次式に同じである。

この左辺は二つの関数の積の微分形になっているので、

となる。これは変数分離して積分ができる。

よって

（, として置換積分をおこなう）

初期条件　より　

よって。これを変形・整理して

最後にについては、この反応の両論係数から全成分の濃度の和は一定である：

ことよりとなる。

以上のように連立常微分方程式はごく簡単なものでも解析解を得るのはたいへんである。

【実習】上の問題を初期濃度 C_A=1, C_R=0, C_S=0 として、k₁, k₂をパラメータとして積分をおこなう。

表計算+VBAで連立常微分方程式をシミュレーションする。rx05.xls　雛形 rx05_temp.xls

結果を　k₁= k₂, k₁> k₂, k₁< k₂ の場合についてレポートしなさい。

Sub SUB2000()
'---------ﾊﾟﾗﾒｰﾀ
k1 = Cells(1, 2)
k2 = Cells(2, 2)
'---------ﾋﾞﾌﾞﾝﾎｳﾃｲｼｷ
Select Case NK
Case 0
' [ DYDX=Y(0)' ]
DYDX = -k1 * Y(0)
Case 1
' [ DYDX=Y(1)' ]
DYDX = k1 * Y(0) - k2 * Y(1)
Case 2
' [ DYDX=Y(2)' ]
DYDX = k2 * Y(1)
End Select
End Sub

【参考】生態系におけるロトカ・ヴォルテラモデル

ある地域のウサギの個体数をx ，これをエサとするキツネの個体数をｙとする。それらの変化は被食者と捕食者の生存闘争モデルとして有名なロトカ・ヴォルテラ（Lotka-Volterra）の微分方程式：

で表される。この式は，キツネがいなければウサギの個体数は微分方程式 dx/dy=ax にしたがって増え，ウサギがいなければキツネの個体数は微分方程式 dy/dt=-by にしたがって減少すると仮定し，生物種の相互作用の項を cxy としてモデル化したものである。

ウサギの個体数が多いとキツネの個体数も増加が可能となる→キツネが大量に増えてウサギの個体数は減少する→餌食となるウサギが減少すればキツネの数も減少する→キツネが少なくなるとウサギの個体数が増加に転じる、という周期的変化をシミュレートすることができる。

境界層方程式（「流体工学」より）

相似変数

により境界層方程式が次のように常微分化される。

（境界条件：　　および）

なお、,

上式は、と置くことで連立常微分方程式：

となる。各式をη=0から同時に数値積分すれば解は求められる。しかしη=0の位置でのy₂=f''の値が未定である。その代わりにη＝∞でのf'の値が境界条件として与えられている。このような問題を「境界値問題」と言う。一方，数値積分は積分計算開始点の値がすべて既知の問題すなわち「初期値問題」しか解けない。どうするか？それは簡単で，η=0での値を仮定してとりあえず積分を実行し，η=∞での条件に合うまでこの方法を繰り返すのである。これをShooting-Methodと言う。実際にはη＝∞はη＝10にとり，f''(0)(=ｙ₂(0))を0.3から0.4程度にとって積分して，　f’(10)＝１になるまで試行する。

【実習】境界層方程式に挑戦！fl03.xls 雛形　fl03_temp.xls

解を求めて速度分布を描く。

Select Case NK
Case 0
' [ DYDX=Y(0)' ]
DYDX = Y(1)
Case 1
' [ DYDX=Y(1)' ]
DYDX = Y(2)
Case 2
' [ DYDX=Y(2)' ]
DYDX = -0.5 * Y(0) * Y(2)
End Select
End Sub

振動の微分方程式

単振動の微分方程式は重力とフックの法則より、

である。（k: バネ定数）

この２階の常微分方程式の解析解は６章で述べた。ここではこれをルンゲ・クッタ法により解く。そのため上式を、

と置いて、

のように連立常微分方程式とする。初期条件は、

（距離１引っ張って離す場合）

または

（静止状態からおもりに速度を与える場合）

【実習】表計算+BASICで振動の微分方程式の解を求める。これはバネ振動のシミュレーションと言える。初期条件を変えて振動の様子をみる。<bane1.xls><bane1_temp.xls>

Sub SUB2000()
'---------ﾋﾞﾌﾞﾝﾎｳﾃｲｼｷ
OM = Cells(2, 2)
Select Case NK
Case 0
' [ DYDX=Y(0)' ]
DYDX = Y(1)
Case 1
' [ DYDX=Y(1)' ]
DYDX = -OM ^ 2 * Y(0)
End Select
End Sub

さらに、右のようにダンパを加えた系を考える。力のバランスは、
（質量Mによる慣性力）＋（ダンパーの制動力）＋（ばねの力）＝（外力）
より、
速度vを位置ｘで書き換えると、v =dx /dtより、

ここで、操作変数：外力f、被制御変数：位置x。

外力がない場合を考え、y₀=x, y₁=x'
として、次の連立常微分方程式となる。

Sub SUB2000()
'---------ﾋﾞﾌﾞﾝﾎｳﾃｲｼｷ
M = Cells(2, 2)
D = Cells(3, 2)
k = Cells(4, 2)
Select Case NK
Case 0
' [ DYDX=Y(0)' ]
DYDX = Y(1)
Case 1
' [ DYDX=Y(1)' ]
DYDX = (-D / M) * Y(1) + (-k / M) * Y(0)
End Select
End Sub

実習

上の系（バネマスダンパ系）につき、バネ定数k, ダンピング定数D, 質量Mを変えて振動の様子をレポートする。
Ａ：kの影響　kが大きいと振幅が（大きくなる、小さくなる）
Ｂ：Dの影響　Dが大きいと振動の減衰が（速くなる、遅くなる）
Ｃ：Ｍ影響　Mが大きいと振動周期が（短くなる、長くなる）