オイラー法を使って,微分方程式を計算するプログラムを作成せよ.以下,詳細に説明し
ているので,これを良く読めばプログラムの作成ができるはずである.
次の微分方程式
をオイラー法により計算する.ただし,初期条件は
![$ x=0$](img3.png)
の時
![$ y=0$](img42.png)
とする.当然,これは
数値計算するまでもなく,解は
![$\displaystyle y=\sin x$](img43.png) |
(2) |
と分かっている.分かっているが,ここではコンピュータープログラムにより計算する.
ここでの内容を良く理解すれば,通常解けない微分方程式でも,コンピューターにより計
算できることが分かるだろう.
コンピューターで微分方程式を計算する方法を示す.微分方程式の解を
とする.
すると,微分--正しくは導関数--は,
![$\displaystyle \frac{\mathrm{d}y}{\mathrm{d}x}\simeq\frac{\Delta y}{\Delta x}=\frac{f(x+\Delta x)-f(x)}{\Delta x}$](img45.png) |
(3) |
と近似することができる.右辺の
![$ \Delta x\rightarrow 0$](img46.png)
の極限が微分
の値となる.したがって,元の微分方程式は
![$\displaystyle \frac{f(x+\Delta x)-f(x)}{\Delta x}\simeq\cos x$](img47.png) |
(4) |
と書くことができる.これと,式(
1)から,
![$\displaystyle f(x+\Delta x)\simeq f(x)+\Delta x\cos x$](img48.png) |
(5) |
である.これがオイラー法で微分方程式を解くときの基本の式である.
![$ \Delta x$](img49.png)
の値を
小さくすると,この近似の精度が上がる.初期条件を上手に使うと,微分方程式の解の値
が芋づる式にわかる.仮りに,
![$ \Delta x=0.001$](img50.png)
とすると,次のような手順で計算する.
|
(1) |
![$ f(0.000)=0$](img51.png) |
これは初期条件である. |
|
(2) |
![$ f(0.001)=f(0.000)+0.001\times\cos(0.000)$](img52.png) |
右辺は(1)の結果を利用する |
|
(3) |
![$ f(0.002)=f(0.001)+0.001\times\cos(0.001)$](img53.png) |
右辺は(2)の結果を利用する |
|
(4) |
![$ f(0.003)=f(0.002)+0.001\times\cos(0.002)$](img54.png) |
右辺は(3)の結果を利用する |
|
(5) |
![$ f(0.004)=f(0.003)+0.001\times\cos(0.003)$](img55.png) |
右辺は(4)の結果を利用する |
|
|
![$ \vdots$](img56.png) |
これを繰り返す |
(1)は初期条件なので計算するまでもない.そして,(1)の結果を利用すると,(2)の右辺
が計算でき,その左辺の値を決めることができる.同様に(2)の結果を利用すると,(3)の右辺
が計算でき,その左辺の値が決められる.これを繰り返すのである.すると,
![$ x=0$](img3.png)
の初
期条件からはじめて,任意の
![$ x$](img5.png)
まで
![$ f(x)$](img1.png)
の値を求めることができる.元の微分方程式
(
1)の近似解が求まったことになる.
![$ \Delta x=0.001$](img50.png)
は
計算のステップ幅と呼ばれ,これを小さくするとさらに計算時間は必要になるが,計算精度が上がる.
例えば,これを4000回この計算を繰り返すと,微分方程式(1)
の解が
から
まで計算できる.すなわち,0〜
(229.1[度])までの値が
(0.057[度])きざみで分かるのである.4000
回の計算なんか,コンピューターに取っ手は大したことはない.私のPCでは,計算と表示
に用した時間は,0.15秒である.コンピューターの計算速度には,本当に驚かされる.
計算の原理が分かったので,プログラミング方法を示す.計算のフローチャートは図
2のようになる.これにしたがって,プログラムを書けば良い.難
しいことは何もない.
まずは,計算のステップ幅
を決めなくてはならない.C言語では
dx=4.0/4000;
と書く.プログラムではギリシャ文字は使えないので,ステップ幅を
dxとしている.
解となる計算結果は後で利用することを考えて,配列に格納する方が良い.配列の先頭に
は,初期条件を格納する.すなわち,
x[0]=0.0;
y[0]=0.0;;
とするのである.次に
iの値を
1〜
4000まで変えて,
x[i]=i*dx;
y[i]=y[i-1]+dx*cos(x[i-1));
を繰り返し計算する.このように計算回数が決まっている繰り返しには,
for文を使
うのが一般的である.
for(i=1;i<=4000;i++){
ここに繰り返したい内容を書く.
}
ホームページ:
Yamamoto's laboratory著者:
山本昌志
Yamamoto Masashi
平成18年7月18日