So-net無料ブログ作成

ねこ騙し数学 ルンゲ=クッタ法で連立微分方程式を解く [ネコ騙し数学]

ねこ騙し数学 ルンゲ=クッタで連立常微分方程式を解く

 

ルンゲ=クッタ(Runge-Kutta)法で二階常微分方程式の初期値問題が少し人気のようなので、連立常微分方程式を解く方法をご紹介しますにゃ。

二階の常微分方程式の一般形は、

ルンゲクタで連立微分方程式_htm_m76489ee2.gif

ルンゲクタで連立微分方程式_htm_m56201d10.gifは初期値、計算の出発点だにゃ。

二階常微分方程式で紹介したのは、実は、この特殊なものなんだにゃ。

 

十進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

 

で、以前紹介した2階常微分方程式のプログラムを少し修正するだけで、簡単に出来てしまうだにゃ。

十進BASICをコンピュータにインストールしてあれば、このプログラムをコピペすれば、
ルンゲクタで連立微分方程式_htm_m76489ee2.gif

のタイプの連立常微分方程式は、解けてしまうんだんにゃ。

ここで解いているのは、

ルンゲクタで連立微分方程式_htm_m8df0aa.gif

なんだけれども、
これをこのプログラムで数値的に解くと、

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

 

となる。

この連立微分方程式の厳密解は、

ルンゲクタで連立微分方程式_htm_5752668c.gif

になるのですが、厳密解と比較すると、

ルンゲクタで連立微分方程式_htm_77e4c589.gif

 

 

赤い曲線がの厳密解で、青い曲線がの厳密解なんですが、
h = 0.1 
と非常に粗いものでも x = 2.2 くらいまでは厳密解とよく一致していることがわかるにゃ。

厳密解を見るとわかるのだけれど、x =√ 7 = 2.645・・・ で不連続、しかも、この点で∞に発散してしまう。

数値計算は、こういう不連続点近傍での計算が苦手。

不連続点の近くでは、どうしても誤差が大きくなってしまう。


ですが、

ルンゲクタで連立微分方程式_htm_5f6283dc.gif

といった連立微分方程式であれば

ルンゲクタで連立微分方程式_htm_m61f832ed.gif

といった具合に、高精度で計算できるにゃ。

=0.1 と非常に粗いのに、この精度、文句なしでしょ。

ちなみに、この厳密解は、

ルンゲクタで連立微分方程式_htm_m68367644.gif

プログラムの修正箇所は


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でも書け」
と書き込んでいただければ、
記事に付け足すにゃ。


nice!(0)  コメント(2)  トラックバック(0) 
共通テーマ:学問

nice! 0

コメント 2

k

Cもお願いします
by k (2018-03-10 05:41) 

nemurineko

コメント、ありがとうございます。

現在、このブログは休眠中なので、お返事が遅くなり、申し訳ありません。
C言語のプログラムは、コチラに出ています。

http://nemuneko-gensokyou.blog.so-net.ne.jp/2016-11-24-13

よろしかったら、ご覧になってください。

また、Web版ですと、タグクラウドに数値解析というタグがあり、それをクリックしていただけると、数値解析に関係した過去の記事一覧をご覧になれますので、よろしかったら、ご利用ください。
by nemurineko (2018-03-30 02:34) 

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

トラックバック 0

この広告は前回の更新から一定期間経過したブログに表示されています。更新すると自動で解除されます。