第7回(2009年6月10日)


今日の演習


1階常微分方程式の数値解法(ルンゲ・クッタ法)

今回は次のような常微分方程式の初期値問題を、ルンゲ・クッタ法と呼ばれる方法を用いて数値的に解き、GLSCを用いて結果をグラフにしてもらいます。

           (1)

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

 ルンゲ・クッタ法の導出については、このページの最後を参照してください(詳細は参考書を参照してください)。ルンゲ・クッタ法は次にしめすようなアルゴリズムで、オイラー法に比べれば多少複雑なものとなっています. 比較の為に両方のアルゴリズムを載せておきます(赤字の部分が異なります). よく比較してください.


< オイラー法 >

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

変数は y とし、a を初期値とします。y_new は実数型の変数です。


T,N,h を定める

a の入力を受ける (printf 関数 + scanf 関数)

y = a

t = 0 でのグラフの描画 (g_move 関数)

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

    t = i*h

    y_new = y + h*f(t,y)

    y = y_new

    一定間隔毎にグラフを描画 (g_plot 関数)

を繰り返す

f(t,y) は関数として定義するとよい.


< ルンゲ・クッタ法 >

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

変数は y とし、a を初期値とします。y_new, r_1, r_2, r_3, r_4 はそれぞれ実数型の変数です。


T,N,h を定める

a の入力を受ける (printf 関数 + scanf 関数)

y = a

t = 0 でのグラフの描画 (g_move 関数)

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

    t = i*h

    r_1 = f(t, y)

    r_2 = f(t + h/2, y + (h/2)*r_1)

    r_3 = f(t + h/2, y + (h/2)*r_2)

    r_4 = f(t + h, y + h*r_3)

    y_new = y + (h/6)*(r_1 + 2*r_2 + 2*r_3 + r_4)

    y = y_new

    一定間隔毎にグラフを描画 (g_plot 関数)

を繰り返す

f(t,y) は関数として定義するとよい.


* オイラー法では, (とびとびの)時刻 t = jh での値 y から 次の(とびとびの)時刻 t = (j+1)h での値 y_new を直接計算しています。本来微分方程式では時刻は連続なのですが、この場合は時刻 t = jh 〜 (j+1)h の間の事は一切考慮されておらず、その分誤差が現れてきます.

* ルンゲ・クッタ法では, この t = jh 〜 (j+1)h の間の情報を (r_1, 〜 r_4 を使って ) いくらか取り入れる事で、よりよい近似計算が可能になっています. (勿論誤差は現れるのですが、オイラー法と比べて雲泥の差が現れます. 以下の課題で実感できるのではと思います.)


課題1

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

 オイラー法を用いた解法プログラムを ~/work/euler.c とし、同問題をルンゲ・クッタ法を用いて解くプログラムを ~/work/rk.c とする。解 y(t) = exp(t) とどの程度あっているか、一目でわかる GLSC プログラムにせよ。

 結果は大体次のようになります(このグラフは GLSC でかいたものではありませんが大体同じようなグラフをGLSCで描いてください)。h=0.5とかなり粗く時間刻み幅をとっているにも関わらず、ルンゲ・クッタ法は良く微分解を近似できていることがわかります。それに比べてオイラー法は微分解から大きく離れてしまっています。これでは正しい計算とは言えないでしょう。

グラフ


2階(n階)常微分方程式の数値解法(ルンゲ・クッタ法)

前回の課題で扱った




          ・・・ (1)


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





          ・・・ (2)


を再び考えます。このような連立の常微分方程式をルンゲ・クッタ法を用いて解く場合のアルゴリズムを次に示します。

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

変数は y1, y2 とし、それぞれの初期値を a1, a2 とします。r1_1, r1_2, r1_3, r1_4, r2_1, r2_2, r2_3, r2_4 および y1_new, y2_new はそれぞれ実数型の変数です。

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

(今回, ω = 1 とします.)


T,N,h を定める

a1, a2 の入力を受ける (printf 関数 + scanf 関数)

y1 = a1

y2 = a2

t = 0 でのグラフの描画 (g_move 関数)

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

   t = i*h

    r1_1 = f1(t, y1, y2)
    r2_1 = f2(t, y1, y2)

    r1_2 = f1(t + h/2, y1 + (h/2)*r1_1, y2 + (h/2)*r2_1)
    r2_2 = f2(t + h/2, y1 + (h/2)*r1_1, y2 + (h/2)*r2_1)

    r1_3 = f1(t + h/2, y1 + (h/2)*r1_2, y2 + (h/2)*r2_2)
    r2_3 = f2(t + h/2, y1 + (h/2)*r1_2, y2 + (h/2)*r2_2)

    r1_4 = f1(t + h, y1 + h*r1_3, y2 + h*r2_3)
    r2_4 = f2(t + h, y1 + h*r1_3, y2 + h*r2_3)

    y1_new = y1 + (h/6)*(r1_1 + 2*r1_2 + 2*r1_3 + r1_4)
    y2_new = y2 + (h/6)*(r2_1 + 2*r2_2 + 2*r2_3 + r2_4)

    y1 = y1_new
    y2 = y2_new

    一定間隔毎にグラフを描画 (g_plot 関数)

を繰り返す

f1(t,y1,y2), f2(t, y1,y2) は関数として定義するとよい.


課題2

 ルンゲ・クッタ法を用いて(1)を解け。前回(第6回)の課題2同様のグラフを、0≦t≦20の範囲において、h=0.0001,0.001, 0.01, 0.1それぞれの場合について解 y(t)=cos(t) とどの程度あっているか各自調べよ。


課題3

 これまでは、横軸を t 、縦軸を y としたグラフを描いてもらいました。例えば横軸を y、 縦軸を y' としたグラフを描いてみて下さい。


レポート課題

今回の課題(1〜8)で作成したルンゲ・クッタ法を用いたプログラムのうちの一つを rep0610.c というファイル名で report フォルダにおいてください. (プログラムの先頭にコメント文で名前と学籍番号、どちらの課題のプログラムかを書く事. またプログラムの説明もコメント文として適宜記入する事.)

自動回収日時:2009年6月16日(火)午後4時


前回の課題4(ルンゲ・クッタ法でもやってみてください. とくに h をあまり大きくしない場合で, オイラー法の時と比較してみてください. )

下図のような棒の長さが一定の振り子の、重力下での運動について、振り子の角度を 2πX とすると、X は

X'' = -sin(2πX)

 という運動方程式を満たす(摩擦は無視しています)。この方程式を適当な初期条件、X(0) = 0, X'(0) = a でシミュレーションし、GLSC を用いた振り子の動画を作成せよ.(デモをお見せします.)


課題5

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

dx/dt = 1 - (1 + b)*x + x2y

dy/dt = b*x - x2y

(1) b < 2 の場合
(2) b > 2 の場合

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


課題6

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

dx/dt = 3x - xy

dy/dt = -y + xy

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

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


課題7

 惑星の運動を考えます。万有引力の法則から2物体(質量はそれぞれ M,m)に働く力は

wpe3.jpg (1439 バイト)

となります。rは2物体間の距離、Gは万有引力定数で

wpe4.jpg (2817 バイト)

です。原点に太陽(質量 M)があり、惑星(質量 m)の位置を

wpe1.jpg (3009 バイト)

とします。惑星に働く力は

wpe6.jpg (5166 バイト)

となり、

wpe7.jpg (1519 バイト)

より

wpe8.jpg (1269 バイト)

とおくと、2階の常微分方程式の初期値問題

wpe9.jpg (4810 バイト)

が得られます。これを1階の常微分方程式系に直し、数値計算にて解をもとめて下さい。
(まずは r3 を考えずに, r1 - r2 平面上の運動方程式として考えてもよろしいです.)


課題8

水星は太陽に非常に近く、その効果が公転に影響するため、その運動は以下のような方程式に従うようになる。

課題7同様、数値的に解をもとめてください。ただし h << g とします。


ルンゲ・クッタ法について

 簡単にルンゲ・クッタ法について説明します。演習で取り扱ったルンゲ・クッタ法は4段公式とよばれるものですが、ここでは簡単の為、2段公式を取り扱います。

 次の常微分方程式の初期値問題の数値解を求めたいとしましょう。

式(1)   (R1)

2段公式とは、 で(R1)の近似解 が求められたとして、時間ステップ幅 h だけ進んだ における値 を次の形で計算する公式です。

式(2)    (R2)

つまり、からへ進む1ステップにおいて方程式(R1)の第一式の右辺のを2回計算するので、2段公式と呼びます。この2段公式には、パラメータ が含まれますが、これらは方程式(R1)の真の解 x(t) のテイラー展開

    (R3)

と(R2)式の第一式の同様の展開とが h のなるべく高い次数の項まで一致するように定めます。つまり、と仮定して(R2)式の の式の右辺を 周りでテイラー展開すると、

     (R4)

となります。これと(R3)がの項まで一致するように

式(5)     (R5)

とおく。公式(R2)で定められる近似解のテイラー展開が h の2次の項まで真の解のテイラー展開に一致するという意味で、この公式を2次の公式と呼びます。この公式は2段公式でもあるので、2段2次の公式とも呼ばれます。パラメータを定める方程式(R5)には自由度が2残っているのでその定め方によって色々な公式が導かれます。(改良オイラー法、ホイン法、修正オイラー法など)

演習で用いた4段4次のルンゲ・クッタ法の公式も同様にして求めることが出来ます。


戻る