第10回(2011年6月22日)


今日の演習


常微分方程式の解のアニメーション (3)

微分方程式の自由度(変数の数)が多くなってくると、解の様子をみるにはアニメーションを作成すると便利です。今回は実際に、バネで結合した多数のおもりの運動を例に、自由度の多い微分方程式系の振る舞いを観察してみます。


課題0 中間テストのおさらい

中間テストの問題1、2を解答できなかった方は、まずその復習から始めてください.

問題はこちら.

問題1のサンプルプログラムはこちら.


課題1(以前やっていなかった人はやってみてください)

x軸上にバネを2つ並べ、おもりでつなげた場合どのような運動をするか、数値的に計算してみましょう。ただし、全てのバネの自然長R(伸びも縮みもしていないときの長さ)とバネ定数が同じで、おもりの質量も同じとします。




左のおもりを 1, 右のおもりを 2, とします. どのバネも伸びも縮みしていない状態であれば、m番目のおもりの位置はバネの自然長Rを用いて、m R (m = 1 or 2) と書けます。ここで、おもり 1 の R からのズレを y, おもり 2 の 2R からのズレを z とします。

(1) バネの堅さ ω = 1 として, 右のおもりが z(t) = sin( ωo t ) と変化する場合, おもり1の R からのズレ y(t) は以下の式に従って変化します。

  y'' = - 2 y + sin( ωo t )

y(0) = 0, y'(0) = 1 とし, 適当なωoo = 0.5 ~ 2) について解き、2つのおもりの位置 ( y(t)+R , 0 ), ( z(t)+2R , 0 ) の変化をアニメーションで表示せよ.
(共鳴(または共振 ... タコマナローズ橋、マイクのハウリングで起きた(起きる)事)を体験してみてください. )

(2) おもり2も自由に振動する場合, y, z の変化は次の式に従います。

  y''= -ω2 y + ω2( z - y )

  z''= -ω2( z - y )

y(0) = 0, z(0) = 0, とし 適当な y'(0), z'(0), ω を与えて解き, 2つのおもりの位置 ( y(t) + R , 0 ), ( z(t) + 2R , 0 ) の変化をアニメーションで表示せよ.


課題2 (先に課題5以降をやっても良いかもしれません)

x軸上にバネを多数並べ、おもりでつなげた場合どのような運動をするか、数値的に計算してみましょう。ただし、全てのバネの自然長R(伸びも縮みもしていないときの長さ)とバネ定数が同じで、おもりの質量も同じとします。


左のおもりから順に 0, 1, ..., M, M+1 と番号を振ります。全てのバネが伸びも縮みもしていない状態であれば、m番目のおもりの位置はバネの自然長Rを用いて、m R と書けます。これを用いると、m番目(m≠0,m≠M+1)のおもりの運動は、m R からのズレ ym は次のように書けます。



・・・ (1)

0 番目と M+1 番目のおもりは固定されて動かないとします。つまり、



          ・・・ (2)

とします。


プログラムを書く際の注意

これまでオイラー法を用いて2階の微分方程式をを解く際に

y'' = f(t, y, y')

に対し、y1 = y, y2 = y' などと置き、

y1' = y2
y2' = f( t , y1, y2 )

として解いてきました。しかしおもりの位置そのものが y1 , y2 , ... , と添字つきで表される場合、やや混乱が生じるおそれがあります。このような場合、例えば課題0のサンプルプログラムでもしていますように、y1' = v1 , y2' = v2 , ... などと置き、

yj' = vj
vj' = f( yj+1 , yj , yj-1 )

などとした方が見易くなります。

また今回考える方程式系では、おもり j が従う式の右辺は yj+1 , yj , yj-1 のみの関数なので, そのように置きます.



(1) M=3、ω=1、R=1 として、式(1)と式(2)を解くプログラムを作成せよ。初期条件は、y1=0.25、y2=y3=0、y'm=0.0 とする。適度なステップ毎に、おもりの位置(m番目のおもりの位置は ( m R + ym , 0 ) を g_circle 等で表示させること。次図は M = 1 の場合で, グラフにする横軸の範囲は[0,2]とした。




(2) M=30、ω=1、R=1 として式(1)と式(2)を解き、横軸 mR , 縦軸 ym としたアニメーションを作成せよ。また、初期条件を様々に変えて、その様子を観察せよ

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

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

           (1)

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

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


< オイラー法 >

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

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


T,N,h を定める

y = a

t = 0 でのグラフの描画

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

    t = i*h

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

    y = y_new

    一定間隔毎にグラフを描画

を繰り返す

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


< ルンゲ・クッタ法 >

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

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


T,N,h を定める

y = a

t = 0 でのグラフの描画

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

    一定間隔毎にグラフを描画

を繰り返す

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 を使って ) いくらか取り入れる事で、よりよい近似計算が可能になっています. (勿論誤差は現れるのですが、オイラー法と比べて雲泥の差が現れます. 以下の課題で実感できるのではと思います.)


課題3

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

 オイラー法を用いた解法プログラムを euler.c とし、同問題をルンゲ・クッタ法を用いて解くプログラムを 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 を定める

y1 = a1

y2 = a2

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

    一定間隔毎にグラフを描画

を繰り返す

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


課題4

上記課題0、課題1をルンゲ・クッタ法を用いて実施せよ。

課題5

(ルンゲクッタ法がお勧めですが、オイラー法でみもちろんできます。)
地球の周りの人工衛星の運動について考えます。人工衛星と地球は互いに引力が働くが、(現代の)人工衛星にくらべ地球は非常に重いので、事実上人工衛星は下図のように引っ張られながら運動する。その起動が下図ののような2次元平面上を動くとすると(この仮定は実際に正しい)、その位置( x , y )は以下のような方程式

  x'' = - x / r3
  y'' = - y / r3
  ただし r = ( x2 + y2 )1/2

に従って変化する。この方程式を初期条件 x(0) = 0 , y(0) = 1 , x'(0) = a , y'(0) = b , とし、いろいろな a , b について、人工衛星の(楕円)軌道(下図(b))をアニメーションで観察せよ。

(a)


(b)

課題6

(ルンゲクッタ法がお勧めですが、オイラー法でみもちろんできます。)
地球が2つありそれらの周りを人工衛星が運動する場合を考えます(下図(a))。1つめの地球が原点、2つめの地球が位置 (2, 0) にいる場合、人工衛星の位置( x , y )は以下のような方程式

  x'' = - x / r13 - (x - 2) / r23
  y'' = - y / r13 - y / r23

  ただし r1 = ( x2 + y2 )1/2
      r2 = ( (x - 2)2 + y2 )1/2

に従って変化する。この方程式を初期条件 x(0) = 0 , y(0) = 1 , x'(0) = a , y'(0) = b , とし、いろいろな a , b について、人工衛星の軌道をアニメーションで観察せよ。

(a)


(b)

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

 簡単にルンゲ・クッタ法について説明します。演習で取り扱ったルンゲ・クッタ法は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次のルンゲ・クッタ法の公式も同様にして求めることが出来ます。


戻る