計算数学演習第１３回

ルンゲ・クッタ法

これまでオイラー法を用いた数値解法の演習を行いました。
オイラー法は単純でよいのですが、実際の解に近い数値解を出すためには、時間刻み幅 h をかなり小さくとらねばならないことが前回の課題等によってわかりました。
そこで、より精度の高い解法が必要となりますが、今回はそのような解法として、ルンゲ・クッタ法を扱い、それを用いた数値解法の演習を行います。
ルンゲ・クッタ法は様々な次数（精度）の方法がありますが、1次の方法はオイラー法と一致します。他にも2次や4次の方法がよく使われますが、導出については参考書やネットなどを参照してください（期末レポートに関係）。この授業で扱う4次のルンゲ・クッタ法は次にしめすようなアルゴリズムで、オイラー法に比べれば多少複雑なものとなっています。

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

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

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

N を定める(N=(int)(T/h)+1)）
y = a
j = 0,1,....,N-1の順に
    t = j*h // 「現在の」時間
     (一定間隔毎に)t と y の表示 (printf関数)
    //以下, t, y から y_new を求めるのに, 計算精度を上げるためやや遠回り（4次のルンゲ・クッタ法）をします.
    //具体的には 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 は次ステップの時刻に対応していましたが、ルンゲ・クッタ法では f を多用し、かつ、+h/2などもするので、表記と説明を簡単化するために t を現在の時刻としています。
もちろん、tを次ステップの時刻として、書き換えることも可能で、実際にやっていることは同じです。

課題1

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

(1)   \begin{eqnarray*} \dot{y} &=& y \\ y(0) &=& 1 \end{eqnarray*}

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

結果は大体次のようになります。

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)   \begin{equation*} \left\{ \begin{array} $\ddot{y}= -y $\\ y(0) = 1 \\ \dot{y}(0) = 0 \end{array} \right. \end{equation*}

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

(3)   \begin{equation*} \left\{ \begin{array} $\dot{y}_1 = y_2$ \\ \dot{y}_2 = -y_1 \\ y_1(0) = 1, y_2(0) = 0 \end{array} \right. \end{equation*}

を再び考えます(\ddot{y} は y の時間に関する２階微分)。
このような連立の常微分方程式をルンゲ･クッタ法を用いて解く場合のアルゴリズムを次に示します。
ただし、
上の1階微分方程式のときと同様に、T を解を求める時刻の上限、h を時間刻み幅として設定し、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 が(3)式の第一式の右辺、 f2 が(3)式の第二式の右辺となります。

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

N を定める
y1 = a1
y2 = a2
j = 0,1,....,N-1の順に
    t = j*h
     (一定間隔毎に) t と y1, y2 を画面に表示 (printf関数)
    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
を繰り返す

（変数が多くて定義するだけでも面倒そうですね… 配列を使うと楽になるかも？→どうすれば良いか考えてみましょう。）

課題２

次の微分方程式

(4)   \begin{eqnarray*} \dot{y}&=& v \\ \dot{v}&=& -y \end{eqnarray*}

を
y(0) = 1、v(0)=0、時間の刻み幅h=0.2としてルンゲクッタ法でとき、(y(t), v(t))の軌道をgnuplotで描画せよ。
このとき、同じくh=0.2のもとオイラー法で計算した結果と比較すること（下図参照）。

課題３

次の微分方程式

(5)   \begin{eqnarray*} \dot{y}&=& v \\ \dot{v}&=& -v + y -y^3 \end{eqnarray*}

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

ここで、h=0.01として、(y(t+h),v(t+h))から距離0.0001以内に(y(t),v(t))が存在すれば（すなわち\sqrt{(y(t+h)-y(t))^2+(v(t+h)-v(t))^2}<0.0001なら）、収束しているとみなすことにします。
（ヒント：収束する点はyが正のもの、yが負のものの2点が見つかるはずです。）

様々な初期値に対する軌道を下に示します。下の図のような結果を得るにはどのようなプログラムを書けば良いのでしょうか？

応用問題

課題3で、
どのようなy(0),v(0)を選べば、y が正の点に収束するのか。
その初期値(y(0),v(0))の集合をAとして、
Aをy–v平面上に図示するためには、どのようなプログラムをである必要があるか考えてみよ。

課題４

次の微分方程式

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

をルンゲクッタ法で数値的に解いて、A = 0.1と固定して、w = 0.61, 0.95, 1.0, 1.05の場合の y(t) の時間変化を0\le t\le 400の範囲で調べよ。
ただしy(0) = 1.0 , v(0) = 0.0 とする。(h=0.01とし、各時刻データを出力せよ。)

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

課題５

次の微分方程式

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

を適当な初期値から数値的に解いて、B = 0, 1, 6, 12 の場合の y の時間変化を調べよ。

これは、バネが大きく伸び縮みすると急激に硬くなるような場合の、強制振動系です。特に y, v の初期値が少し違う場合の y の時間変化を比較してください。
課題３とどのような特徴の違いがあるでしょうか。
ハードディスクから原子炉まで、様々な機械の「振動」が、このような方程式を応用しながら解析されています。

課題６

次の連立微分方程式を適当な初期条件のもと、以下の場合についてオイラー法かルンゲ・クッタ法で解け。

(8)   \begin{eqnarray*} \dot{x}&=& 1 - (1+b)x + x^2 y \\ \dot{y}&=& 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*} \dot{x}&=& 3x-xy \\ \dot{y}&=& -y+xy \end{eqnarray*}

課題６の方程式同様、この方程式も振動解を見せますが、その特徴は少々異なります。
例えば初期条件を変えた場合に、これらの（課題６と課題７それぞれの）方程式系の解の振る舞いが、どのように変化していくか、見比べてください。
（オイラー法の場合は時間刻み幅hを非常に小さくする必要があります。）

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

三次元のグラフの描画

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

　splot x*x-y*y

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

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

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

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

サンプルファイルsin3d_data.datを使って実際に作図を行ってみます。

num_data.datの中身は、

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

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

Department of Mathematical and Life Sciences, Graduate School of Integrated Sciences for Life, Hiroshima University