第11回(2012年1月11日)


今日の演習


今日の目標


常微分方程式の数値解法


次のような常微分方程式の初期値問題を、オイラー法を用いて数値的に解き、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)式に置き換えられました。この方法をオイラー法と呼びます。見やすいように(10)を書き直すと(11)になります.

   (11)


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

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

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

N, j は整数型(int)

.(11)式の 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 + 3exp(-t), y(0) = 1

について.

(1) 手計算にて解を求め、gnuplot でグラフ化せよ.

(2) オイラー法を用いて 時刻 0≦ t ≦20 の範囲での解を求めるプログラムを作成せよ.時刻 0≦ t ≦20 の範囲を20000等分割して(N = 20000)、100ステップ毎に printf関数により、t および y の値を表示し,それを gnuplot でグラフ化せよ.また (1) の結果と比較し、正しいか確認せよ.

まずは(上のアルゴリズムを参考に、)自力でプログラムするようにしてください。 どうしてもできない場合は、ここにサンプルプログラムを置いときますので、上のアルゴリズムとよく見比べて、何をしているのかしっかり理解してください。(これまで紹介してきたC言語に関する内容の総復習になります。)


課題2

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

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

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

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

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


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

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

 

式(1)   (11)

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

式(2)   (12)

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

ただし、

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

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

f1(t, y1, y2) および f2(t, y1, y2) はそれぞれ関数で(12)式の場合には f1 が第一式の右辺、 f2 が第二式の右辺となります.


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

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)

    y1 = y1_new

    y2 = y2_new

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

を繰り返す

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


課題3

(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 でもうまく解を近似できるような解法が求められます.そのような要求に応えるルンゲ・クッタ法があります(後で触れます).


gnuplot 付け足し

課題3で載せているグラフには、タイトルと軸、各曲線に名前がついてます.このように名前をつけるには、次のようにします.(データファイル名はそれぞれ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 "名前" と付けることで、各曲線に自由な名前が付けられます


グラフの印刷の仕方

gnuplot で作成したグラフを印刷したい場合には,ps ファイルに保存すると便利です.その場合の手順は次のようになります.(例えば sin(x) のグラフが描かれた ps ファイル 「sin_graph.eps」 を作成する場合.)

gnuplot 起動中に

set term postscript eps
set output "sin_graph.eps"
plot sin(x)
quit

これによって,「sin_graph.eps」 というファイルができあがります.ここでは、一行目で出力形式の変更、2行目で出力先のファイル名の指定、をしています. ファイル名や plot の部分を適宜変更して下さい. (plot ... の前に set title ... 等 (上の「gnuplot 付け足し」参照.) をすると、タイトル、軸の名前付きのグラフができます.)

作成した epsファイルをファイルブラウザでダブルクリックすると、適当なpsファイルビューアが起動するので、それから印刷してください.(最寄りのプリンターから印刷されます.)ちなみに eps ファイルとは PostScriptファイルの略です.PostScriptとはプリンタの制御言語で,プリンタに絵を描かせる命令の集まりです.


課題4

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

バネとおもり

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

式(1)   (13)

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

式(2)   (14)

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


課題5

次の微分方程式

dy/dt = v

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

を、y, v に適当な初期値を与え数値的に(オイラー法、ルンゲ・クッタ法のどちらかで)解いて、例えば B = 0, 12 の場合の y の時間変化を調べて下さい(時刻の上限を少し大きめにして下さい.).また、y と v の関係(gnuplot で横軸を y, 縦軸を v として)も plot してみて下さい. これは課題3で紹介した摩擦抵抗のあるバネ振動子に、周期的な外力を加えたもので、強制振動系と呼ばれます.


課題6

次の微分方程式

dy/dt = v

dv/dt = - 0.3*v - y3 + B*cos(t)

を数値的に解いて、B = 0, 1, 6, 12 の場合の y の時間変化を調べてください.また課題3同様、y と v の関係(gnuplot で横軸を y, 縦軸を v として)も plot してみて下さい. これは、バネが大きく伸び縮みすると急激に硬くなるような場合の、強制振動系です.特に y, v の初期値が少し違う場合の. y の時間変化を比較してください. 課題3とどのような特徴の違いがあるでしょうか.

ハードディスクから原子炉まで、様々な機械の「振動」が、このような方程式を応用しながら解析されています.


今日学んだ事


戻る