第5回、第6回(2011年5月18日、5月25日)


今日の演習


常微分方程式の数値解法


今回は次のような常微分方程式の初期値問題を、オイラー法と呼ばれる方法を用いて数値的に解き、GLSCを用いて結果をグラフにしてもらいます。オイラー法は計算数学演習で扱っていますが, おさらいをしておきます. 次の問題を考えます。

           (1)

求めたいのは、t>0の範囲のy(t)であり、関数 f(t, y) および定数 a は与えられているとします。また、 f(t ,y) の値は t, y が与えられればいつでも計算できるものとします。(関数として定義するとよい。)

コンピューターでは数値を離散的に(とびとびに)しか扱うことができません。そこで、(1)式をコンピューターで扱えるようにする必要があります。(1)をコンピュータで解く方法には色々とありますが、今回は代表的なものの一つであるオイラー法を用います。その前に、準備として差分商について説明します。


差分商

 f(x) を x の関数としたとき、f(x) の x=a における微分係数は

         (2)

で定義されます。ここで、定義中の極限操作を取り払い、h を有限にとどめた

         (3)

を考えると、h を十分小さな正の実数ととれば、(3)式は(2)式の近似となっていると考えられます。この D のように、関数のいくつかの点における値の差を用いてその関数の微分係数を近似することを差分近似といい、(3)式の右辺のような量を差分商といいます。いまの場合は1階微分係数を近似する1階差分商です。

では、このような差分商を用いて(1)式を離散化してみましょう。まず、(1)式の第一式の微分係数を上記の差分商で置き換えてみます。

         (4)

ただし、h は正の数とします。(4)式と(1)式の第一式は異なる方程式なので、(1)式の第一式の解 y(t) は一般に(4)式を満たしません。そこで、混乱を防ぐために(4)式の y(t) を Y(t) と書き換えましょう。

         (5)

(5)式のように差分を含む方程式を差分方程式といいます。

 次に、(1)式の y を Y に置き換えた初期条件

         (6)

の下で(5)式を解くことを考えます。(5)式は、

         (7)

と変形できるので、t = 0 を代入すると、

         (8)

となり、Y(0) から直ちに Y(h) が求まります。同様にして(7)式を繰り返し用いると、

         (9)

というように、j=1,2,3,.... とすると Y((j+1)h) を Y(jh) から計算できることがわかります。

 h をいったん決めると、 t=jh 以外の時刻の Y の値は(7)式からは求めることができません。このように、差分方程式を解くと従属変数はとびとびの時刻で値が定まります。そのようなとびとびの時刻を格子点と呼びます。Y は格子点でのみ意味があるので、そのことを明示するために Y(jh) を Yj と書き換え、 Tj = jh とすると、(5)式と初期条件は、

         (10)

となり、結局(1)式の常微分方程式の初期値問題が、Yjに関する漸化式の問題(10)式に置き換えられました。この方法をオイラー法と呼びます。見やすいように(10)を書き直すと(11)になります.

   (11)


オイラー法のアルゴリズム1

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

T, h, a 及び 関数 f は実数型(float)

N, j は整数型(int)

Y[N+1], t[N+1] は実数型配列(float)


T, N, aを設定 (#define文で定数として定義すれば良い)

配列 Y[N+1], t[N+1] を用意する

h = T/N(TおよびN双方が整数型でないか注意)

Y[0] = a (初期値の設定)

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

t[0] = 0

j = 0,1,2,.......,N-1 の順に (for 文)

t[j+1] = (j+1)*h

Y[j+1] = Y[j] + h*f(t[j], Y[j])

以上を繰り返す

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

点(t[0], Y[0]), (t[1], Y[1]), ... , (t[N], Y[N]) を直線で結ぶことにより数値解をグラフとして描画 (g_move, g_plot 関数)

f(t,Y) は関数(上で復習してます)として定義するとよい.

計算数学演習での「関数の定義」の解説とサンプルプログラムはこちらにあります.


以下のアルゴリズム2では,数値解を求めながらグラフを描いています.(アルゴリズム1では,数値解を全て求めた後グラフを描いている.)


オイラー法のアルゴリズム2

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

T, h, a, Y, Y_new, t 及び 関数 f は実数型(float)

N, j は整数型(int)


T,N,aを設定

h = T/N

Y = a

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

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

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

j = 0,1,2,.......,N-1 の順に (for 文)

t = j*h

Y_new = Y + h*f(t, Y)

一定間隔毎にグラフを描画 (g_plot 関数を用いて (t+h, Y_new) に線をつなげる )

Y = Y_new ( Y の値を更新 )

以上を繰り返す

f(t,Y) は関数(上で復習してます)として定義するとよい.

計算数学演習での「関数の定義」の解説とサンプルプログラムはこちらにあります.


簡単な常微分方程式とその数値解のグラフ化


さて,実際にオイラー法をつかって常微分方程式を解いてみましょう.

まず簡単ですが非常に重要で、講義でも扱っている次のような常微分方程式を用いて演習を進めます。

 (P)

この微分方程式は解析的にも簡単に解けるのですが,今回は(P)をオイラー法を用いて数値計算で解いてもらい,得られた数値解をグラフ化してもらいます.


計算数学演習での「オイラー法の解説とサンプルプログラム」はこちらにあります ので 参考にしてください.
上記サンプルと、前回までのGLSCでのグラフの描画をうまく組み合わせて、以下の「課題1」に取り組んでください.
以下の課題では,なるべくアルゴリズム2を用いるようにしてください,(分かりにくい人はアルゴリズム1を用いても構いません.)


課題1

(P)をオイラー法で解き,得られた数値解をグラフとして表示するプログラムを作成せよ.ただし,

Y0=1, λ=1, 0≦t≦10

とする。例えば次のようになる。(この例では、0≦t≦10 を 100等分している。つまり時間刻み幅を h と書くことにすると、 h = 10.0/100 = 0.1 である。)(縦軸がN(t)となっているがy(t)と読み替えてください.)

解析的に求めた解をグラフとして表示したものは以下であるが,殆ど違いがない事に注意.(課題2でみるが,微妙に異なっている)


課題2

解析的に求めた解(S)とオイラー法で求めた数値解を重ねて描くプログラムを作成せよ.例えば,次のようになる.この例では,解(S)を緑色で描き,数値解を赤色で描いている.(この例では、0≦t≦10 を 100等分している)

刻み幅を小さくすると両者の差は縮まる.(次の例では、0≦t≦10 を 300等分している)


課題3

次の常微分方程式の初期値問題をオイラー法を用いて解いて、t と y の関係を GLSC でグラフ化せよ.

 (1) dy/dt = -y + sin(t), y(0) = 1  (時刻 0≦ t ≦20 の範囲を20000等分割、100ステップ毎に値を表示)

 (2) dy/dt = y(5-y), y(0) = 1  (時刻 0≦ t ≦20 の範囲を20000等分割、50ステップ毎に値を表示)

 (3) dy/dt = y - 2y3, y(0) = 0.1 もしくは y(0) = -0.1 (時刻 0≦ t ≦20 の範囲を20000等分割、100ステップ毎に値を表示)


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


2階常微分方程式もオイラー法を用いて解くことができます。次の問題を考えます。



          ・・・ (7)


講義の説明の通り(下の枠中)、これは単振動を記述する2階の常微分方程式です。ただし簡単の為、時間 t での微分を ' を用いてあらわしています。

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



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

          ・・・ (8)


ただし、
          ・・・ (9)


とした。

2階の常微分方程式は次のように連立の1階の常微分方程式に帰着できます。y1(t)=y(t), y2(t)=y'(t) とおくと、式(7)は、




          ・・・ (10)


となります。この形になれば、前回紹介したオイラー法を用いて解を求めることができます。アルゴリズムは次のようになります。但し、 T を解を求める時間の上限とし、h を時間刻み幅、N を時間刻み数とします。また、変数は y1, y2 として、a1, a2 を初期値とします。また、f1(t, y1, y2) および f2(t, y1, y2) はそれぞれ関数で式(10)の場合には f1 が第一式の右辺、 f2 が第二式の右辺となります。


オイラー法のアルゴリズム

解を求める時刻の上限を 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

    一定間隔毎にグラフを描画 (g_plot 関数を用いてに線をつなげる)

を繰り返す

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


課題4

式(7)を解析的に解け。


課題5

ω=1として式(10) をオイラー法で解き、グラフを描け。グラフは横軸が t であり、縦軸を y(t) とせよ。また、時間刻み幅 h を変えたときに結果がどのようにかわるか観察せよ。0≦t≦20の範囲において、h=0.0001,0.001, 0.01, 0.1それぞれの場合について、解析解とどの程度あっているか各自調べよ。

色々な h に対する計算結果は次のようになる(緑はh=0.1、青はh=0.01、ピンクはh=0.001、黄色はh=0.0001のときの結果。赤は解析解のグラフ)。この結果から、h は十分小さい値でないとまずいことがよく分かる。




オイラー法では h を十分小さくとらないと本当の解をうまく近似できません。h を小さく取るということは、計算時間が多くかかることを意味します。そこでより高精度で大きな h でもうまく解を近似できるような解法が求められます。そのような要求に応える解法としては、ルンゲ・クッタ法があります。



課題6

課題5のプログラムを用いて、


          ・・・ (11)


を計算せよ。結果は、横軸をt、縦軸を E とするグラフで確認すること。この E は m=k=1 のときの、バネの力学エネルギーに対応する。


課題7

実際のバネではバネの復元力以外におもりの速度 y' に比例する抵抗が働く。(おもりと床との間の摩擦力のようなものをイメージ)抵抗を考慮すると、式(8)ではなく次の式になる。(2γy' の項が新たに加わった。)



          ・・・ (12)


ここで γは抵抗力の強さを表す正の定数です。式(12)の問題を連立1階常微分方程式に変換し、オイラー法を用いて解き、GLSCをつかって「課題5のようなグラフ」を作成せよ. 各自 γ < ω での様子と γ > ω での様子を見比べること。

レポート課題(5/25分)

次の2階の微分方程式

X'' = - sin(2πX)

を、初期条件として X(0) = 0, X'(0) = 0.7 を与えオイラー法で解き、X と t について「課題1のようなグラフ」を作成するプログラムを作成せよ。

ファイル名を rep0525.c とし、以下の提出用ページより提出する事。

回収用ページ:提出期限2011年5月31日(火)午後0時


課題8(現実の動きを動画で再現。)

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

X'' = - sin(2πX)

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


課題9

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

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

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

(1) u, v について「課題5」のようなグラフを作成せよ。
(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) の動き方 )がどのように変わるか考察せよ.

課題10

課題8の振り子が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)

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


戻る