第10回 微分方程式と差分法:Runge-Kutta法(12月6日)



*出席確認

<出席確認メール> 添付ファイルなし 受付:本日14:30〜16:15

本日のキーワードは授業中に連絡します。


宛先:miwamoto[at]riko.shimane-u.ac.jp

件名:MMM2 20161206

本文: 学生番号 名前 「本日のキーワード」


*研究室訪問について

明日12月7日(水)13時から解析系セミナー室(総合理工1号館8階)

見学する人は、13時に必ず集合してください。途中退出可です。



オイラー法による2階常微分方程式解法


2階常微分方程式をオイラー法を用いて解く方法を考える。

 

 

 


 これは単振動を記述する2階の常微分方程式である。2階の常微分方程式は次のように連立の1階の常微分方程式に帰着できる。

 ,    とおくと、


 

 

 

 


と書き換えることができる。

これにより、1階常微分方程式と同様、オイラー法を用いて解を求めることができる。




2階常微分方程式のオイラー法



<変数の宣言>

 y1, y2, y1_new, y2_newをfloatで用意する。


<setup関数で>

初期値の設定 y1 = a1, y2 = a2


<draw関数で>

    y1_new = y1 + h*f1(t, y1, y2)

    y2_new = y2 + h*f2(t, y1, y2)

    y1 = y1_new

    y2 = y2_new

 

練習1 次の常微分方程式の初期値問題をオイラー法をオイラー法で解き、グラフを描け。グラフは横軸が t であり、縦軸を y(t) とせよ。いろいろな異なる時間刻み幅 h について計算し、結果について考察せよ。



 

 





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


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



4次ルンゲ・クッタ法アルゴリズム


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


<draw関数で>


    t = frameCount*h

    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



練習2 4次ルンゲ・クッタ法を用いて、次の問題を実際に解け。



オイラー法とルンゲクッタ法の精度


練習2の微分方程式をオイラー法とルンゲクッタ法で解いた結果を比べてみる。

左図は、刻み幅h=0.01とした場合に、赤は解析解(y = exp(t))であり、青がオイラー法、緑がルンゲクッタ法である。ルンゲクッタ法の精度が高いことがわかる。

右図は、刻み幅をh=0.5とした場合である。オイラー法の精度がかなり悪いのがわかる。一方、ルンゲクッタ法は、h=0.5とかなり粗く時間刻み幅をとっているにも関わらず、良く解析解を近似できていることがわかる。




ルンゲクッタ法の導出


 ルンゲ・クッタ法はテイラー展開を用いて、オイラー法よりも近似精度をあげたものである。練習で扱ったのは4次ルンゲ・クッタ法(4段公式)である。

詳しくは、こちら



2階微分方程式とルンゲクッタ法


練習1で扱った



 

 



 ,    とおいて次のように連立の1階の常微分方程式に帰着したもの


 

 

 

 


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

ただし、 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 が第二式の右辺となる。


<draw関数で>


   t = frameCount*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



練習3 ルンゲ・クッタ法を用いて次の微分方程式を解け。



 

 

0≦t≦20の範囲において、h=0.0001,0.001, 0.01, 0.1それぞれの場合について微分解 y(t)=cos(t) とどの程度あっているか各自調べよ。


練習4 これまでは、横軸を t 、縦軸を y としたグラフを描いていたが、横軸を y、縦軸を y' としたグラフを描け。(練習3の微分方程式について描いてみよ。)


練習5(発展) 重力の作用のみを受けて運動する物体を考える(真空中でボールを投げるような状況)。ニュートンの運動の第二法則から物体の時刻 t での鉛直方向の位置を y(t) とするとその物体の鉛直方向での運動は次のような常微分方程式で表される(下向きを正とする)。

    (5)

m は物体の質量、g は重力加速度である。色々な初期値  y(0), y'(0) を与えてシミュレーションしてみよ。


練習6(発展) さらに、上の問題に空気の抵抗のようなものを加味してみましょう。物体にかかる抵抗がそのときの物体の速度に比例するとすると、運動方程式は、

   (6)

となる。ただし、r は抵抗の強さを表す。色々な初期値および m, r で計算してみよ。

*なお、練習5と6は解析的に解ける問題なので、実際に解いてみて、今回のシミュレーション結果と見比べて理解を深めてください。


練習7(発展) 1960年代初頭に気象学者エドワード・ローレンツは流体力学に基づく次のような単純な3変数系の常微分方程式を示した。

このモデルはいわゆるカオスを生じるモデルとして有名である。パラメータは P, R, b の3種類があるが、P=16, b = 4, R=35として、u-w平面に解の軌道をプロットすると次のような図ができあがる。(初期値 (u0, v0, w0) = (5, 5, 5), t=0 から t=100 までの計算結果を線でつないだもの)

 

これはストレンジアトラクター(奇妙なアトラクター)と呼ばれるもので、有名なので目にしたことがある人も多いと思います。このようなカオス的な解の挙動が天気予報があたらない理由となっているともいわれています。例えば P=16, b=4 として R を色々と変えてみると、解の挙動があれこれ変わります。色々と試してみると面白いでしょう。



本日の課題1

下の図のように壁に一端を固定されたバネの他端におもりがつけられている状況を考える。

このおもりを少しだけ右に引っ張ってから手をはなす.おもりの中心の時刻 t におけるつり合いの位置からのずれを y(t) とする。また、時刻 0 でのおもりの位置を y(0)=1、おもりの速度を y'(0)=0 とすると、練習1と全く同じ単振動の問題となる。

これではおもりに働く力はバネの復元力のみで、ある意味非現実的なので、今回は復元力以外におもりの速度 y' に比例する抵抗が働く場合を考えましょう(おもりと床との間の摩擦力のようなものをイメージ)。

改良されたモデルは次のようになる(2ky' の項が摩擦をあらわしている)。

 

 

 

ここで k,ω はそれぞれ抵抗力、復元力の強さを表す正の定数である。

この問題を連立1階常微分方程式に変換し、オイラー法を用いて解き、グラフにせよ。

k < ω でのグラフと k > ω でのグラフを見比べて違いに関して考察してみて下さい。




本日の課題2

課題1を4次ルンゲクッタ法を用いて解き、グラフにせよ。


本日の課題3

課題1の微分方程式を用いて、おもりが動くアニメーションを作成せよ。


このブラウザは埋め込まれたオーディオファイルを再生できません




このブラウザは埋め込まれたオーディオファイルを再生できません




このブラウザは埋め込まれたオーディオファイルを再生できません


本日の課題4 横軸を y、縦軸を y' としたグラフを描いてみよ。