第12回

オイラー法

今日の演習

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

今日の目標

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

常微分方程式の数値解法

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

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

:11_files:1.gif （１）

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

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

差分商

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

（２）

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

（３）

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

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

:11_files:4.gif （４）

ただし、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に関する漸化式の問題（１０）式に置き換えられました．この方法をオイラー法と呼びます．

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

解を求める時刻の上限を 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 を設定(#define 文)

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

（ここから計算開始）

初期値を与える 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) は関数として定義するとよい．

ただし，#defineを使ってしまうと，aの値の入力をキーボードから受けてとることができない．

その場合は，aをdefineせずに

また，上のプログラムは「時刻の上限Tを明示的に指定しているが，時間刻み幅hはT/(N-1)により間接的に与えている」．

hは計算の精度に直結するので，hを直接指定したい場合もあるだろう．

その場合は，「時刻の上限TはT*hと間接的に指定しているが，時間刻み幅hは直接値を指定する」

ということになる．どちらも一長一短であるので，好きに使い分けると良い．（アルゴリズムとしては同等である）

後者の記法で，なおかつ初期値をキーボードから受け取る例を下に示す．

（例２）

時刻の下限を0 , 時間刻み数(格子点数) N を設定(#define)．

時間刻み幅 h を指定(陽に0.01などと指定する）(#define)

（ここから計算開始）

aの入力を促す

初期値を与える 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("%lf %lf\n", t+h, Y_new);
}

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

課題１

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

　dy/dt = -y + 3exp(-t),　y(0) = 1

を解け．

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

　　解答例へのリンク

課題２

次の常微分方程式の初期値問題をオイラー法を用いて解いて、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 - 2y³,　y(0) = 0.1 もしくは　y(0) = -0.1　(時刻 0≦ t ≦20 の範囲で時間刻み幅h=0.001、100ステップ毎に値を表示)

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

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

（１1）

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

（1２）

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

ただし、

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

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

f1(t, y1, y2) および f2(t, y1, y2) はそれぞれ関数で（１２）式の場合には 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) は関数として定義するとよい．

課題３

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

:11_files:graph1.gif 　

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

レポート課題

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

:11_files:omori.gif

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

:11_files:a1.gif （１３）

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

（１４）

ここで k,ω はそれぞれ抵抗力、復元力の強さを表す正の定数です．課題は、（１４）式の問題を連立１階常微分方程式に変換し、オイラー法を用いて解くプログラムrep20170120.cを作成してください．

ただし，ω=1.0, h=0.0001とし，出力はt, y(t), y'(t)の３成分で，

0.0 y(0.0) y'(0.0)

0.01 y(0.01) y'(0.01)

0.02 y(0.02) y'(0.02)

...

49.99 y(49.99) y'(49.99)

50.0 y(50.0) y'(50.0)

というふうに，tは0.01刻みで出力するようにしてください．（ヒント:100回に１回出力するので，%を使うと良いでしょう）

kの値について：y(t)のグラフが下図のようになるkの値を見つけ，そのkの値を使ってプログラムを作成してください．

rep

レポート提出ページ（提出期限2017年1月27日午前10時）

今日学んだ事

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