第8回(2011年6月8日)


今日の演習


中間テストのお知らせ


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

微分方程式の自由度(変数の数)が多くなってくると、解の様子をみるにはアニメーションを作成すると便利です。今回は主にこのアニメーションを通じて様々な方程式系の振る舞いを観察してみてください。(アニメーション作成の方法は 2-4 回演習のページを参考にしてください.)


課題0

バネ定数kのバネの一端を壁に固定し、他端に質量mのおもりをつけた状況を考える。



このおもりを少しだけ右に引っ張ってから手をはなす。おもりの中心の時刻 t におけるつり合いの位置からのずれを y(t) とすると。y(t)が満たす微分方程式は

 y'' = -ω2 y 

であり, y1(t) = y(t), y2(t) = y'(t) とおくと上式は

 y1' = y2
 y2' = -ω2 y1

となる.
y(0) = 0 (初期の釣り合いの位置からのからのずれが 0), y'(0) = 1 (おもりの初期速度が 1), ω = 1 であるとして、おもりの位置 ( y(t) , 0 ) の変化をアニメーションで表せ.
(ヒント) 課題0アルゴリズム.

解を求める時刻の上限を T とし、適当な自然数 N を定めて h=T/N とする。

T, h, a1, a2, y1, y2, y1_new, y2_new, t 及び 関数 f1, f2 は実数型(float)

N, j は整数型(int)


T,N,h を定める

初期値を与える y1 = a1, y2 = a2

g_init() 関数, g_device() 関数, g_def_scale()関数 等で描画環境の準備

g_sel_scale() で描画環境(仮想座標系)の選択

t = 0 でのグラフの描画 (曲線の始点を与える g_move 関数)

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

    t = j*h

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

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

    y1 = y1_new

    y2 = y2_new

    一定時間間隔毎(例えば j が100の倍数である場合に)に、
       1. g_cls 関数で描画されているものを消す
       2. g_circle 関数等で ( y1 , 0 ) に丸を描画
       3. g_sleep 関数で少々停止
   を行う

を繰り返す

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


課題1

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

X'' = - sin(2πX)

 という運動方程式を満たす。この方程式を適当な初期条件、X(0) = 0, X'(0) = a を与えオイラー法で解き(注意:時間刻み幅を十分小さくすること)、GLSC を用いた振り子の位置( sin(2πX) , -cos(2πX) ) )の変化の動画を作成せよ.(デモをお見せします.)


課題2

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 ) の変化をアニメーションで表示せよ.


課題3

課題1の振り子が2つ相互作用する場合を考える. 例えば, 2つの振り子の角度を 2πX1, 2πX2 とすると X1, X2, が

 X1'' = - sin(2π(X1-X2)) - sin(2πX1)

 X2'' = - sin(2π(X2-X1)) - sin(2πX2)

 という運動方程式を満たすとする(二つの振り子が引き合っている)。この方程式を適当な初期条件、X1(0), X1'(0), X2(0), X2'(0) を与えオイラー法で解き、GLSC を用いた振り子のアニメーションを作成せよ. (注意:時間刻み幅を十分小さくすること)

また振り子の数を3つにした場合

 X1'' = - ( sin(2π(X1-X2)) + sin(2π(X1-X3)) ) - sin(2πX1)

 X2'' = - ( sin(2π(X2-X3)) + sin(2π(X2-X1)) ) - sin(2πX2)

 X3'' = - ( sin(2π(X3-X1)) + sin(2π(X3-X2)) ) - sin(2πX3)


さらに A,B,C は適当な定数, これらをいろいろな値にした場合

 X1'' = - ( A sin(2π(X1-X2)) + C sin(2π(X1-X3)) ) - sin(2πX1)

 X2'' = - ( B sin(2π(X2-X3)) + A sin(2π(X2-X1)) ) - sin(2πX2)

 X3'' = - ( C sin(2π(X3-X1)) + B sin(2π(X3-X2)) ) - sin(2πX3)

や、もっと多くした場合も考察せよ.


課題4

次の「ブリュッセレーター(ブラッセレーター)モデル(方程式)(これが何を表すものか分からない人は、「BZ反応」や「Brusselator」等でインターネット等で検索してみよ.)」

 du/dt = 1 - (1 + b)*u + u2v
 dv/dt = b*u - u2v 

を、様々な b に対してオイラー法を用いて解き、

(1) u, v について、横軸 t、縦軸U, V としたグラフを作成せよ。
(2) u-v 空間(横軸 u, 縦軸 v) での座標(u,v)の時間変化のアニメーションを作成せよ。

(3) また、

 du/dt = 1 - (1 + b)*u + u2v
 dv/dt = b*u - u2v + d*(x-v) 
 dw/dt = 1 - (1 + b)*w + w2x
 dx/dt = b*w - w2x + d*(v-x)

を、様々な b, d に対してオイラー法で解き、座標(u,v),(w,x)の時間変化を同一平面上のアニメーションとして描け. b,d 変えた場合に、解の挙動((u,v),(w,x) の動き方 )がどのように変わるか考察せよ.

課題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)

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


課題6

 上記のルンゲ・クッタ法を用いて、次の問題を実際に解き、ルンゲ・クッタ法がオイラー法と比べてどの程度性能が良いか比べよ。次の初期値問題を 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) は関数として定義するとよい.


課題7

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

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

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


戻る