計算数学演習第１３回

ルンゲ・クッタ法

これまでオイラー法を用いた数値解法の演習を行いました．
オイラー法は単純でよいのですが、実際の解に近い数値解を出すためには、時間刻み幅 h をかなり小さくとらねばならないことが前回の課題等によってわかりました．
そこで、より精度の高い解法が必要となりますが、今回はそのような解法として、ルンゲ・クッタ法を扱い、それを用いた数値解法の演習を行います．

ルンゲ・クッタ法の導出については、このページの最後を参照してください（詳細は参考書を参照してください）．ルンゲ・クッタ法は次にしめすようなアルゴリズムで、オイラー法に比べれば多少複雑なものとなっています（赤字で示したところが一階微分方程式に関するオイラー法からの変更箇所です.）．

一階常微分方程式に対する ルンゲ・クッタ法

T を解を求める時刻の上限とし、h を時間刻み幅、N を時間刻み数とします．
また、変数は y とし、a を初期値とします．r1, r2, r3, r4 はそれぞれ実数型の変数です．

T,N を設定 (#define 文）

初期値 a を設定 (#define 文）

h を定める

y = a

t = 0 での y の値の表示 (printf 関数)

j = 0,1,2,....,N-1の順に

    t = j*h

    //以下, t, y から y_new を求めるのに, 計算精度を上げるためやや遠回りをします.
    //具体的には t, y から r1 を求め, t, y, r1 から r2 を求め, ... t, y, r3 から r4 を求め, 最後に　y, r1, r2, r3, r4 から y_new を求める.
    //この遠回りのおかげで、計算精度が非常に良くなります. )

    r1 = f(t, y)

    r2 = f(t + h/2, y + (h/2)*r1)

    r3 = f(t + h/2, y + (h/2)*r2)

    r4 = f(t + h, y + h*r3)

    y_new = y + (h/6)*(r1 + 2*r2 + 2*r3 + r4)

    y = y_new

    (一定間隔毎に)t+h と y_new の表示 (printf関数)

を繰り返す

（二階常微分方程式に対する ルンゲ・クッタ法）

先の課題で扱った

(1)   \begin{equation*} \left\{ \begin{array} $y''= -y $\\ y(0) = 1 \\ y'(0) = 0 \end{array} \right. \end{equation*}

をy1(t)=y(t), y2(t)=y'(t) とおいて次のように連立の１階の常微分方程式に帰着したもの

(2)   \begin{equation*} \left\{ \begin{array} $y_1' = y_2$ \\ y_2' = -y_1 \\ y_1(0) = 1, y_2(0) = 0 \end{array} \right. \end{equation*}

を再び考えます(y” は y の時間に関する２階微分)．
このような連立の常微分方程式をルンゲ･クッタ法を用いて解く場合のアルゴリズムを次に示します．

ただし、

T を解を求める時間の上限とし、h を時間刻み幅、N を時間刻み数とする.

変数は y1, y2 とし、それぞれの初期値を a1, a2 とします．r11, r21, r31, r41, r12, r22, r32, r42 および y1_new, y2_new はそれぞれ実数型の変数とする．

f1(t, y1, y2) および f2(t, y1, y2) はそれぞれ実数型の関数で（２）式の場合には f1 が第一式の右辺、 f2 が第二式の右辺となります．

T,N を設定 (#define 文）

初期値 a1, a2 を設定 (#define 文）

h を定める

y1 = a1

y2 = a2

t = 0 での y1, y2 の表示 (printf 関数)

i = 0,1,2,....,N-1の順に

    t = i*h

    r11 = f1(t, y1, y2)

    r12 = f2(t, y1, y2) 

    r21 = f1(t + h/2, y1 + (h/2)*r11, y2 + (h/2)*r12)

    r22 = f2(t + h/2, y1 + (h/2)*r11, y2 + (h/2)*r12)

    r31 = f1(t + h/2, y1 + (h/2)*r21, y2 + (h/2)*r22)

    r32 = f2(t + h/2, y1 + (h/2)*r21, y2 + (h/2)*r22)

    r41 = f1(t + h, y1 + h*r31, y2 + h*r32)

    r42 = f2(t + h, y1 + h*r31, y2 + h*r32)

    y1_new = y1 + (h/6)*(r11 + 2*r21 + 2*r31 + r41)

    y2_new = y2 + (h/6)*(r12 + 2*r22 + 2*r32 + r42)

    y1 = y1_new

    y2 = y2_new

    (一定間隔毎に) t+h と y1_new, y2_new を画面に表示 (printf関数)

を繰り返す

 

練習課題

課題1

上記のルンゲ･クッタ法を用いて、次の問題を実際に解き、ルンゲ･クッタ法がオイラー法と比べてどの程度性能が良いか比べよ．
次の初期値問題を 時刻 0≦t≦5 の範囲で h=0.5 とし、実際にオイラー法とルンゲ･クッタ法を用いて解け．

(3)   \begin{eqnarray*} \frac{d}{dt}y &=& y \\ y(0) &=& 1 \end{eqnarray*}

オイラー法を用いた結果と、同じ問にルンゲ･クッタ法を用いた結果が、微分解 y(t) = exp(t) とどの程度あっているか、gunuplot を用いて確認せよ．

結果は大体次のようになります．
h=0.5とかなり粗く時間刻み幅をとっているにも関わらず、ルンゲ･クッタ法は良く微分解を近似できていることが わかります．
それに比べてオイラー法は微分解から大きく離れてしまっています．
これではオイラー法の結果は正しい計算とは言えないでしょう (もちろん h をもっと小くすれば、オイラー法でもよい近似を与えます.)．

 

ここで紹介したものは4次のルンゲ・クッタ法と呼ばれるもので，４次の精度をもっています．
つまり，y(t)=Y(t)として，
ルンゲクッタ法で時間刻み幅hだけ計算したものをY(t+h)とすると，
微分方程式の解y(t+h)をテイラー展開したものと，このY(t+h)はhの４次の項までは一致するということを意味しています．
オイラー法は，hについて１次項までは一致しているので解の精度は1次です．

ルンゲ・クッタ法とは，そもそもなにをやっているのでしょうか？
このルンゲ・クッタ法の仕組みを細かく計算して本当に４次の精度になるか確認したい人は是非やってみてください。
(ニューメリカルレシピ in Cなどが参考になると思います）
２次の精度のルンゲ・クッタ法ともいうべき方法は，中点法などが知られています．

課題2

次の微分方程式

(4)   \begin{eqnarray*} \frac{dy}{dt}&=& v \\ \frac{dv}{dt}&=& -y \end{eqnarray*}

を

y(0) = 1，v(0)=0，時間の刻み幅h=0.2としてルンゲクッタ法でとき，(y(t), v(t))の軌道をgnuplotで描画せよ．

このとき，同じくh=0.2のもとオイラー法で計算した結果と比較すること．

参考までに，原点を中心とした半径1の円と，h=0.2のもとオイラー法で計算した結果を下図に示す．

課題3

次の微分方程式

(5)   \begin{eqnarray*} \frac{dy}{dt}&=& v \\ \frac{dv}{dt}&=& -v + y -y^3 \end{eqnarray*}

をルンゲクッタ法でとき，初期値y(0),v(0)によって，最終的に落ち着く点（収束する点）が変わることを調べよ．

ただし，距離0.0001以内に(y(t),v(t))が存在すれば、(y_fix, v_fix)に収束しているとみなすことにする。

（ヒント：収束する先の点はyが正のもの，yが負のものの二点が見つかるはずである．）

応用問題

課題3で，

どのようなy(0),v(0)を選べば，yが正の点に収束するのか．

その初期値(y(0),v(0))の集合をAとして，

Aをyv平面上に図示するためには，どのようなプログラムをつくる必要があるか考えてみよ．

課題4

次の微分方程式

(6)   \begin{eqnarray*} \frac{dy}{dt}&=& v \\ \frac{dv}{dt}&=& -y + A*\cos(w*t) \end{eqnarray*}

をルンゲクッタ法で数値的に解いて、A = 0.1と固定して、w = 0.61, 0.95, 1.0, 1.05の場合の y(t) の時間変化を0≦t≦400の範囲で調べよ．

ただしy(0) = 1.0 , v(0) = 0.0 とする。(時間の刻み幅hは0.01となるようにせよ)

また課題３同様、y と v の関係（gnuplot で横軸を y, 縦軸を v として）も plot してみて下さい.

（w = 1.0の場合、系自身が持っている振動数と外力の振動数が完全に一致します。このとき、共振と呼ばれる現象がおこり、振幅がt→∞で発散することがあります。確認してみてください。）

課題5

次の微分方程式

(7)   \begin{eqnarray*} \frac{dy}{dt}&=& v \\ \frac{dv}{dt}&=& -0.3*v - y^3 + B*\cos(t) \end{eqnarray*}

を数値的に解いて、B = 0, 1, 6, 12 の場合の y の時間変化を調べてください．
また課題３同様、y と v の関係（gnuplot で横軸を y, 縦軸を v として）も plot してみて下さい.
これは、バネが大きく伸び縮みすると急激に硬くなるような場合の、強制振動系です.特に y, v の初期値が少し違う場合の. y の時間変化を比較してください.
課題３とどのような特徴の違いがあるでしょうか.

ハードディスクから原子炉まで、様々な機械の「振動」が、このような方程式を応用しながら解析されています.

課題６

次の連立微分方程式を適当な初期条件のもと、オイラー法かルンゲ・クッタ法で解いてください。

(8)   \begin{eqnarray*} \frac{dx}{dt}&=& 1 - (1+b)x + x^2 y \\ \frac{dy}{dt}&=& bx - x^2y \end{eqnarray*}

(1) b < 2 の場合

(2) b > 2 の場合

この方程式は、ブラッセレーター方程式といわれるもので、BZ反応と呼ばれる、時間的に化学成分の濃度が振動するような化学反応の振動をモデル化した方程式です。
(BZ 反応とは何でしょう？ http://www.chem.scphys.kyoto-u.ac.jp/nonnonWWW/kitahata/bz_1.html あたりが参考になります。)

課題７

次の連立微分方程式を適当な初期条件のもと、オイラー法かルンゲ・クッタ法で解いてください。

(9)   \begin{eqnarray*} \frac{dx}{dt}&=& 3x-xy \\ \frac{dy}{dt}&=& -y+xy \end{eqnarray*}

課題６の方程式同様、この方程式も振動解を見せますが、その特徴は少々違っております。
例えば初期条件を変えた場合に、これらの（課題６と課題７それぞれの）方程式系の解の振る舞いが、どのように変化していくか、見比べてください。

この方程式は、ロトカ・ボルテラ方程式といわれるもので、生態系内のある被食者と補食者の数 (x, y) 増減関係を簡単にモデル化した方程式です。（各時刻、被食者はある割合（3x）で増えて、補食者に出会うと食べられ減る(-xy)。同様に捕食者は被食者 に出会うとそれを食べて増えて(xy)、ある割合で(-y)で寿命等で減っていく。）

三次元のグラフの描画

guplotで三次元のグラフを描画するためには、splotというコマンドを使います。
Gnuplotを起動して、以下のコマンドを打ってみてください。

　splot x*x-y*y

グラフが描画されたとおもいます。

このままでは、どちらの軸がx軸、y軸なのかが分かりにくいので、

各軸にlabelをつけておくと便利です。(set xlabel “x” 等)

また、数値データの入ったファイルを描画することも可能です。

このサンプルファイルnum_data.txtを使って実際に作図を行ってみる。

num_data.txtの中身は、

0.000000 1.000000 0.000000 0.000000

0.300000 0.955336 0.295520 0.564642

0.600000 0.825336 0.564642 0.932039

…

といった４列のデータ(0.3*i , cos(0.3*i), sin(0.3*i), sin(0.6*i) )が入っている。

　splot “num_dat.dat” using 1:2:3 w l

　splot “num_dat.dat” using 2:1:3 w lp

　splot “num_dat.dat” using 1:2:4 w lp

と打ち込んでみて、結果を確認してみましょう。