Part.10

常微分方程式の数値解法（オイラー法）

内容

  • 1階常微分方程式の数値解法（オイラー法）
  • 2階常微分方程式の数値解法（オイラー法）

目標

  • 1階常微分方程式を解く方法を知る
  • n階常微分方程式についても連立1階常微分方程式に帰着できること知る

常微分方程式の数値解法

次のような常微分方程式の初期値問題を、オイラー法を用いて数値的に解き、gnuplot を用いて結果をグラフにしてもらいます．

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

求めたいのは、ｔ＞０の範囲の関数 ｙ（ｔ）であり、関数 f(t, y) および定数 a は与えられているとします．また、 f(t ,y) の値は t, y が与えられればいつでも計算できるものとします．（プログラム中では関数として定義するとよい．）

コンピューターでは数値を離散的に（とびとびに）しか扱うことができません．そこで、（１）式をコンピューターで扱えるようにする必要があります．（１）をコンピュータで解く方法には色々とありますが、この演習では代表的なオイラー法（とルンゲ・クッタ法）を用います．その前に、準備として差分商について説明します．

差分商

f(x) を x の関数としたとき（ここでの f(x) は上記(1)式中の f とは関係ありません）、f(x) の x = a における微分係数は

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

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

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

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

（５）式のように差分を含む方程式を差分方程式といいます．

次に、（１）式の y を Y に置き換えた初期条件

の下で（５）式を解くことを考えます．（５）式は、

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

となり、Y(0) から直ちに Y(h) が求まります（注意：Y(0)は初期条件として与えられているので既知）．同様にして（７）式を繰り返し用いると、

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

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

となり、結局（１）式の常微分方程式の初期値問題が、Yjに関する漸化式の問題（１０）式に置き換えられました．この方法をオイラー法と呼びます．

入力促進メッセージに関する重要な注意

以下にオイラー法のアルゴリズムを示します．計算結果を printf 文で数値を画面に表示します．また，後に gnuplot でそれらをグラフ化します．そのときに，入力促進メッセージの表示方法が問題となります．すなわち，プログラム名が prog であるときに，

./prog > some.data

として，printf 関数の出力を画面ではなく some.data に書き込むわけですが，このままでは入力促進メッセージ(「初期値を入力してください」といったメッセージ)も画面に表示されず，some.dataに書き込まれてしまいます．これではメッセージの意味がない上に，不要な文字列が some.data に書き込まれ，gnuplot でグラフを作成する際に問題となります．

そこで入力促進メッセージに関しては，通常の printf 関数は使わず，fprintf 関数を用いることで上記の問題に対応します．具体的には，入力促進メッセージについては，例えば，

fprintf(stderr, "a の値を入力してください");

というふうに fprintf(stderr, "文字列");

を用いてください．"文字列" の部分は通常の printf 文と同様に使うことができ，変数内の数値の表示などに使う %d 等の変換文字列も同様に使えます．stderr は文字列の表示先を「標準エラー出力」にせよという意味です．実は，通常の printf 文は，

fprintf(stdout, "文字列");

と同等です．stdout は「標準出力」と呼ばれ，通常これはターミナルの画面です．「標準エラー出力」も通常ターミナル出力ですが，プログラムを以下のように実行したさいには

./prog > some.data

「標準出力」の出力が some.data に書き込まれ「標準エラー出力」の出力は画面に表示されることになります．

ですから，以下のようにしてください．

入力促進メッセージは fprintf(stderr, "文字列"); を用いる

数値計算結果の出力には printf("文字列"); も用いる

オイラー法のアルゴリズム（その１）

実際には、後に示す（その２）が良いが、説明のために（その１）を示す．

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

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

N, j は整数型(int)

T,Nを設定
h = T/(N-1)
a の入力を受ける (fprintf 関数 + scanf 関数）
Y = a
初期時刻0と初期値aを画面に表示（printf関数）
j = 0,1,2,.......,N-1 の順に (for 文）
  t = j*h
  Y_new = Y + h*f(t, Y)
  Y = Y_new
  （一定間隔毎に）t+h と Y_new の値を画面に表示（printf関数）
以上を繰り返す

（一定時間毎に）とあるのは，例えば j が100毎にという意味である．つまり，

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

というふうにする．これだと作成されるデータファイルの容量は1/100となる．これは，先に見たようにある程度の折れ線数で得あれば十分なめらかに見えるので，グラフを描く為のデータとして巨大になりすぎることを防ぐ為である．ここで100という数値は適当に与えた物であって，グラフが十分なめらかに見えかつデータが小さくなるような適当な値を探す必要がある．

オイラー法のアルゴリズム（その２）

その１中の Y_new 変数は実際には不要であるので消去した場合．

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

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

N, j は整数型(int)

T,Nを設定
h = T/(N-1)
a の入力を受ける (fprintf 関数 + scanf 関数）
Y = a
初期時刻0と初期値aを画面に表示（printf関数）
j = 0,1,2,.......,N-1 の順に (for 文）
  t = j*h
  Y = Y + h*f(t, Y)
  （一定間隔毎に）t+hとYの値を画面に表示（printf関数）
以上を繰り返す

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

課題１

常微分方程式の初期値問題

を解け．

オイラー法を用いて 0≦ t ≦20 の範囲での解を求めるプログラムを作成せよ．0≦ t ≦20 の範囲を20000等分割して、100ステップ毎に printf関数により、t および y の値を表示し，それを gnuplot でグラフ化せよ．また手計算にて解を求め、数値計算結果が正しいか確認せよ．

課題２

ロジスティック方程式の解を変数分離にて求め、レポート課題の数値計算により得られた解と比較せよ．

課題３

次の常微分方程式の初期値問題をオイラー法を用いて解け．いろいろな異なる時間刻み幅 h について計算し、結果について考察せよ．

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

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

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

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

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

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

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関数）
を繰り返す

課題４

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

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

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

課題５

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

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

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

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

課題６

次の、ロトカ・ボルテラ競争系について、図書館の書籍、WEB上の情報などをもちいて調べまとめよ．どういった現象を表すモデルであるか？

dx/dt = r1 x(1 - (x + ay)/K1)

dy/dt = r2 y(1 - (bx + y)/K2)

課題７

課題６中のロトカ・ボルテラ競争系モデルを数値計算を用いて解け．

課題８

ある種の生物の増殖に関する大変有効なモデルとして次のロジスティック方程式がある．

dx/dt = rx(1 - x/K) (L)

ここに、x(t) は時刻 t における、生物集団のサイズであり，K > 0 を環境収容力と呼び，r > 0 を内的自然増加率と呼ぶ．

今、 r = 1, K = 100 とし、初期条件として、次の３種類のものを考える．

  1. x(0) = 0
  2. x(0) = 10
  3. x(0) = 150

問題１：

このとき、それぞれの初期値に対する (L) の解を 0≦ t ≦10 についてオイラー法で求め，結果をgnuplotでグラフ化せよ．ただし、N = 10000 （時間刻み数）とする．データは100毎に出せば十分である．グラフは、横軸を t 、縦軸を x(t) とする．次に示す画像を参考にせよ．

問題２：

ロジスティック方程式(L)の解を解析的に求めよ．また得られた解が，確かに(L)の解であることも確かめよ．

ちなみに，解は以下のようになる．

x(t) = K/(1 + C *exp(-r*t))

ここに C は，初期値 x(0) から決まる定数で，

C = K/x(0) - 1

で与えられる．

問題３：

オイラー法で求めた近似解と問題２で求めた真の解を比べたい．h を小さくすると近似解は真の解に近づくであろうがどのように近づくであろうか．問題２で求めた真の解を x(t)，刻み幅 h によってオイラー法で求めた近似解を X(t;h) と書くことにする．今，初期値は 10 すなわち x(0) = 10, X(0;h) = 10 とする．次のような尺度をもって近似解の正確さはかることにする．

E(t;h) = |x(t) - X(t;h)|

E(t;h) は刻み幅が h である時の時刻 t における真の解と近似解の差を与える．ここでは，t = 2 における E を様々な h について調べることにする．すなわち，h = 0.1, 0.01, 0.001, 0.0001 それぞれについて，E(2; 0.1), E(2; 0.01), E(2; 0.001), E(2; 0.0001) を求め次の表を完成させよ．

注意１：gnuplot でも整数同士の演算は整数の範囲で行われるので注意．例えば，plot 100/150*x としてみれば何を言いたいのかわかると思う．正しいグラフを得るには，plot 100/150.0*x とする必要がある．

注意２：E(2;h) の表示には %f よりも %e を用いた方がよいかもしれない．実数型変数内の数値が大変大きいまたは小さい場合には %e の出力の方が正確である．例えば，%e による表示が，1.234e-5 は 1.234*10-5 であり，2.345e12 は 2.345*1012 を意味する．gnuplot のデータファイルにもこの表記は使える．

問題４：

問題３で作成した表を数値データとして打ち込み，gnuplot でグラフ化せよ．また，問題3の結果と問題4で作成したグラフから，近似解と真の解の差は h に関してどのように減少するか考えよ．

今日学んだ事

  • オイラー法による常微分方程式に数値解法