ルンゲ＝クッタで連立微分方程式を解く [数値解析]
ねこ騙し数学 ルンゲ＝クッタで連立常微分方程式を解く
ルンゲ＝クッタ（Runge-Kutta）法で二階常微分方程式の初期値問題が少し人気のようなので、連立常微分方程式を解く方法をご紹介しますにゃ。
二階の常微分方程式の一般形は、
は初期値、計算の出発点だにゃ。
二階常微分方程式で紹介したのは、実は、この特殊なものなんだにゃ。
十進BASICによるプログラムは、
REM 関数の設定
DEF f(x,y,z) = x*y*z
DEF g(x,y,z) = x*y/z
REM 初期値の設定
LET x0 = 1
LET y0 = 1/3
LET z0 = 1
REM
LET n = 16
LET h = 0.1
REM
LET x = x0
LET y = y0
LET z = z0
PRINT x, y, z
REM ルンゲ＝クッタの計算本体
FOR i = 1 TO n
LET k11 = h*f(x,y,z)
LET k21 = h*g(x,y,z)
LET k12 = h*f(x+h/2, y+k11/2, z+k21/2)
LET k22 = h*g(x+h/2, y+k11/2, z+k21/2)
LET k13 = h*f(x+h/2, y+k12/2, z+k22/2)
LET k23 = h*g(x+h/2, y+k12/2, z+k22/2)
LET k14 = h*f(x+h, y+k13, z+k23)
LET k24 = h*g(x+h, y+k23, z+k23)
LET x = x + h
LET y = y + (k11 + 2*k12 + 2*k13 + k14)/6
LET z = z + (k21 + 2*k22 + 2*k23 + k24)/6
PRINT x, y, z
NEXT i
END
で、以前紹介した２階常微分方程式のプログラムを少し修正するだけで、簡単に出来てしまうだにゃ。
十進BASICをコンピュータにインストールしてあれば、このプログラムをコピペすれば、
のタイプの連立常微分方程式は、解けてしまうんだんにゃ。
ここで解いているのは、
なんだけれども、
これをこのプログラムで数値的に解くと、
x y z
1 .333333333333333 1
1.1 .370934142767238 1.03624581918336
1.2 .418896708802566 1.07901936703824
1.3 .480885237687925 1.12961617057169
1.4 .56236002979378 1.18975185867056
1.5 .671707630892471 1.26173957102924
1.6 .822277218503532 1.34876663208168
1.7 1.03623809430333 1.4553415016314
1.8 1.35233413329635 1.58804657579305
1.9 1.84263146266847 1.75687124427831
2 2.65203570314821 1.97772117936245
2.1 4.10260093960148 2.2775157352983
2.2 7.01258316226914 2.70560740397124
2.3 13.8929551819076 3.36293533637003
2.4 35.0309001510378 4.49152623721934
2.5 140.370860126822 6.83926320667073
2.6 1828.56941770406 14.1230967733292
となる。
この連立微分方程式の厳密解は、
になるのですが、厳密解と比較すると、
赤い曲線がy の厳密解で、青い曲線がz の厳密解なんですが、
h = 0.1 と非常に粗いものでも x = 2.2 くらいまでは厳密解とよく一致していることがわかるにゃ。
厳密解を見るとわかるのだけれど、x =√ 7 =2.645・・・で不連続、しかも、この点で∞に発散してしまう。
数値計算は、こういう不連続点近傍での微分方程式の振る舞いを再現するのが苦手。
不連続点の近くでは、どうしても誤差が大きくなってしまう。
ですが、
といった連立微分方程式であれば
、
といった具合に、高精度で計算できるにゃ。
ｈ=0.1 と非常に粗いのに、この精度、文句なしでしょ。
ちなみに、この厳密解は、
DEF f(x,y,z) = y-z+1
DEF g(x,y,z) = y+z+1
LET x0 = 0
LET y0 = 0
LET z0 = -1
です。
もちろん、もっと複雑な連立微分方程式であっても高精度で解けるにゃ。
不定積分で解けない微分方程式であっても、解ける！！
C言語やFortranなどのプログラムも欲しければ、コメントにでも
「⑨ネコ、CやFortranでも書け」
と書き込んでいただければ、
記事に付け足すにゃ。
タグ：微分積分
コメント 0