第12回(2012年1月18日)


今日の演習


今日の目標


一階常微分方程式の数値解法

(先週の内容です. 必要に応じて各自復習してください.)


次のような常微分方程式の初期値問題を、オイラー法を用いて数値的に解き、gnuplot を用いて結果をグラフにしてもらいます.オイラー法に関しては、講義で既に習っていますので、ここでの内容は講義の復習となっています.(多少、記号の書き方が違うかもしれません.)

次の初期値問題を考えます.

           (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) が求まります(注意:Y(0)は初期条件として与えられているので既知).同様にして(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)式に置き換えられました.この方法をオイラー法と呼びます.


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

解を求める時刻の上限を T (下限は 0 )とし、適当な自然数 N (格子点の数) を定めて h=T/(N-1) とする
(N...時間刻み数、h...時間刻み幅 ).

T, h, a, t(時刻),Y(ある時刻での値), Y_new(次の時刻での値)及び 関数 f は実数型(floatまたはdouble)

N, j は整数型(int)

.(10)式の Yj, Yj+1, Yj+2, ... を、 「Yj の値から Yj+1 を求める。」、「次に Yj+1 の値から Yj+2 を求める。」、「次に...。」... と求める代りに、,Y(ある時刻での値), Y_new(次の時刻での値)を使って

「Y から Y_new を求める。(例えば Y が Yj に, Y_new が Yj+1 に対応)」、「Y を Y_new で置き換える。」、
「Y から Y_new を求める。(今度は Y が Yj+1 に, Y_new が Yj+2 に対応)」、「Y を Y_new で置き換える。」、...、
という繰り返しで求める.


Y の初期値 a を設定 (#define 文)

時刻の上限 T(下限は 0 ), 時間刻み数(格子点数) N を設定

時間刻み幅 h = T/(N-1) を設定

(ここから計算開始)

初期値を与える Y = a

初期時刻(t = 0)と初期値 a を画面に表示(printf関数)

各格子点での値を j = 0,1,2,.......,N-1 の順に計算 (for 文)

t = j*h

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


一定間隔毎に)t+h と Y_new の値を画面に表示(printf関数)


次の時刻の計算のため、Y を新しく求まった値(Y_new)で置き換え

Y = Y_new

以上を繰り返す

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


(一定時間毎に)とあるのは,例えば j が100毎に、等という意味である.つまり,

if (j % 100 == 0)
{
    printf("%f %f\n", t+h, Y_new);
}

というふうにする (a % b = [a を b で割った時の余り]).これだと作成されるデータファイルの容量は1/100となる.これは,前々回で見たように、グラフはある程度の折れ線数であれば十分なめらかに見えるので,グラフを描く為のデータとして巨大になりすぎることを防ぐ為である.ここで100という数値は適当に与えた物であって,グラフが十分なめらかに見えかつデータ数が小さくなるような適当な値を探す必要がある.


課題1

前回の復習として、以下の常微分方程式の初期値問題を数値的に解きます.

dy / dt = -y + sin(t)

オイラー法を用いて 時刻 0≦ t ≦20 の範囲での解を求めるプログラムを作成してください.時刻 0≦ t ≦20 の範囲を20000等分割して(N = 20000)、100ステップ毎に printf関数により、t および y の値を表示し,それを gnuplot でグラフ化せよ.(ここにサンプルプログラムを置いときます.



(今日の本題)

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

つづいて、2階常微分方程式をオイラー法を用いて解いてもらいます.次の問題を考えます.

 

式(1)   (11)

これは単振動を記述する2階の常微分方程式です.ただし簡単の為、時間 t での微分を ' を用いてあらわしています.2階の常微分方程式は次のように連立の1階の常微分方程式に帰着できます.y1(t)=y(t), y2(t)=y'(t) とおくと、(11)式は

式(2)   (12)

となります.もう少し一般的な書き方をすれば


 y1' = f1(t, y1, y2) (今回は f1(t, y1, y2) = y2)

 y2' = f2(t, y1, y2) (今回は f2(t, y1, y2) = -y1)

という形になります. こうなれば先程と同様、オイラー法を用いて解を求めることができます. アルゴリズムは次のようになります.

ただし、

T を解を求める時刻の上限 (下限 0) とし、N を時間刻み数、h = T/(N-1) を時間刻み幅とします.

変数は y1, y2 として、a1, a2 を初期値とします.

(12)式の第一式の右辺を f1(t, y1, y2), 第二式の右辺を f2(t, y1, y2) とそれぞれ定義します.


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

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

N, j は整数型(int)


T,N,h を定める

y1 = a1, y2 = a2

初期時刻0と初期値a1, a2を画面に表示(printf関数)

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)

    一定間隔毎t+h, y1_new, y2_new を画面に表示(printf関数)

    y1 = y1_new

    y2 = y2_new

を繰り返す

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

注意

各ステップで, y1_new, y2_new が両方求まってから y1 = y1_new, y2 = y2_new と値を更新するようにします. ここの順序が変わってくると, y1 と y2 で時間がずれてしまうので, 正確な計算になりません.


課題2

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

色々な h に対する計算結果は次のようになる.この結果から、h は十分小さい値でないとまずいことが分かります.

グラフ

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


課題3

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

バネとおもり

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

式(1)   (13)

これではおもりに働く力はバネの復元力のみで、ある意味非現実的です.今回は復元力以外におもりの速度 y' に比例する抵抗が働く場合を考えましょう(おもりと床との間の摩擦力のようなものをイメージ).改良されたモデルは次のようになります(2ky' の項が新たに加わったもので、摩擦のようなものをあらわしている).

式(2)   (14)

ここで k,ω はそれぞれ抵抗力、復元力の強さを表す正の定数です.課題は、(14)式の問題を連立1階常微分方程式に変換し、オイラー法を用いて解き、gnuplotをつかってグラフにしてください.ω = 1 として、k < ω でのグラフと k > ω でのグラフを見比べて違いに関して考察してみて下さい.h を色々と変化させると数値解はどう変化するか? (実行例はこちら)


課題4

(1) 次の微分方程式

dy/dt = v

dv/dt = - 0.3*v - y + cos(w*t)

を、幾つか適当な w の場合について, y, v に適当な初期値を与え数値的に解き, y の時間変化を調べよ(時刻の刻み幅を変えずに上限を少し大きめ(100程度)にして下さい.).実行例はこちら

. この方程式系は、課題3で紹介した摩擦抵抗のあるバネ振動子に、周期的な外力を加えた(振動を与えた)もので、強制振動系と呼ばれます.

( * この方程式では、右辺に4つの変数(y, v, t, w)があるので、関数の定義では
f1(float time, float Y, float V, float W),
f2(float time, float Y, float V, float W)
等とするとよいでしょう. )


(2) w = 0.1, 0.2, 0.3, ... 2.0 の各 w 毎について(1)を数値的に解くプログラムを作成し、下のような、バネの振動の振幅 y の w 依存性を表すグラフを作成せよ。
( 方法:例えば w = 0.1*i ( i は整数で 0 ≦ i ≦ 20 ) とし, i についての繰り返し文を用いて w を変化させつつ, 各 w 毎について(1) を数値的に解くようにプログラムを作成する. )
下の結果は 横軸を w, 縦軸を y ( if文を使って t=70 から t=100 での 時間一定間隔毎の y の値のみを表示 )とした場合、以下のような図が描かれる (データは点でプロットしています.). 各自確認する事. 

グラフ

. この図から、バネを揺らす際には、与える揺れの周期を適度に(バネ自身が持っている固有の周期に近く)調整することで、より大きな振幅での揺れを実現できる事を意味しています。これは共振(共鳴)とをばれる現象で, ブランコで体験できたり、吊り橋が壊れたり、地震での建物の揺れで問題になってたりしています。(普通、吊り橋や建物では、共振が起こらないように設計、材料選択がなされています。)

  


課題5

「ブリュッセレーター(ブラッセレーター)モデル(方程式)」

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

について、(1) b < 2 の場合、(2) b > 2 の場合を、適当な初期条件 (x,y)(幾つか試す事)から数値計算し、解の様子を比較せよ。

課題6

及び「ロトカ・ボルテラモデル(方程式)」

 dx/dt = 3x - xy
 dy/dt = -y + xy

について、適当な初期条件 (x,y)(幾つか試す事)から数値計算し、解の様子を比較せよ。また「ブリュッセレーター(ブラッセレーター)モデル(方程式)」との違いも考察せよ。(ロトカ・ボルテラモデルの計算では、刻み幅を十分小さくする事。)

課題7

「ブリュッセレーター(ブラッセレーター)モデル(方程式)」、「ロトカ・ボルテラモデル(方程式)」が何者かを各自調べてみてください。

gnuplot 付け足し

課題2で載せているグラフには、タイトルと軸、各曲線に名前がついてます.このように名前をつけるには、次のようにします.(データファイル名はそれぞれdata1、data2、data3、data4、とした.)

gnuplot を起動

set title "Euler's method"
set xlabel "t"
set ylabel "y"
plot cos(x), "data1" w l title "h=0.0001", "data2" w l title "h=0.001", "data3" w l title "h=0.01", "data4" w l title "h=0.1"

ここでは一行目でグラフの名前を、2行目で横軸(x軸)の名前を、3行目で縦軸(y軸)の名前を設定しています. また描写の際、後ろに title "名前" と付けることで、各曲線に自由な名前が付けられます


今日学んだ事


戻る