オイラー法による常微分方程式の数値計算
電気情報工学科 ∗ コンピュータシミュレーション
2015 年 7 月 24 日 (金)
概要
物理現象は常微分方程式で表現されることが多々ある．ここでは，常微分方程式
の近似解を，数値計算により求めることを考える．数値計算方法はいくつかあるが，
この講義では特に，オイラー法を取扱う．コンピュータシミュレーションを実践する
ためには，プログラミングのスキルよりも，まず計算原理をしっかり理解することが
重要である．数学が苦手だと言うことは何の言い訳にもならない．少なくとも高校 3
年生程度の数学を応用することでシンプルに取扱うことが可能である．そのため，こ
の講義ではプログラミングよりも計算原理の理解を重視する．
1 講義内容概略
数値計算による常微分方程式の近似解を求めてみる．今回はオイラー法の解説を行い，
数値計算を行なうためのプログラムを作成する．特に以下についてを目標とする．
• オイラー法の計算原理をしっかりと理解し，それらを説明することができる．
• アルゴリズムのフローチャートが書ける．
• C 言語でプログラムを作成し，実際に常微分方程式の近似解を得ることができる．
2 常微分方程式
2.1
偏微分方程式と常微分方程式
多変数関数の微分 (偏微分) を含む微分方程式を偏微分方程式 (partial differential equation) という．これに対して，一変数の関数の微分を含む方程式を常微分方程式 (ordinary
differential equation) という．ここでは常微分方程式，特に 1 階の場合の解の近似値を求め
∗
秋田工業高等専門学校
1
る方法を取扱う．偏微分方程式については，後期に取扱うこととする．ここで取扱う方程
式の形は
dy(x)
= f (x, y)
dx
(1)
であり，この微分方程式を満たす y(x) を求めることになる．話を進める前に，この方程式
が何を意味しているかを考えることにする．式 (1) の左辺は，解 y(x) の導関数となってい
る．即ち，解の曲線の接線を表す．導関数の値が座標 (x, y) の関数になっているので，座
標が決まれば，その場の曲線の傾きが決きまることになる．今，
dy(x)
= sin(x) cos(x) − y cos(x)
dx
(2)
のような常微分方程式の右辺を考える．先ほど述べたように，この式は任意の点における
接線の傾きを与える．場所ごとに接線の傾きが決まっているので，それを x − y 平面に図
示することができる．式 (2) の右辺の値である各座標の傾きを線の傾きで表すと，図 1 の
ようになる．この傾きを方向場と言う．この方向場から解の様子がわかる．
この微分方程式の解析解は，
y = sin(x) − 1 + c1 e− sin(x)
(3)
である．1 階の微分方程式なので，1 個の未知数を含む．この未知数の値が異なる 5 本の
曲線と，先ほどの方向場を重ね書きすると，図 2 のようになる．微分方程式の解である曲
線 y(x) が方向場に沿っているのがよくわかる．
2.2
数値計算のイメージ
まず，(1) 式の微分方程式を極限の d の代わりに有限な ∆ に置き換える．∆ が小さけれ
ば，元の微分方程式の良い近似になるはずである．すると，式 (1) の微分方程式は，
∆y = f (x, y)∆x
(4)
のように近似できる．これを用いて， xi から ∆x 離れた y の値 yi+1 を計算する．
yi+1 = y(xi + ∆x)
= yi + ∆y
= yi + f (xi , yi )∆x
(5)
初期値 x0 , y0 を与えるてこの式を解くと，次々に (x1 , y1 ), (x2 , y2 ), (x3 , y3 ), . . . が計算できる．
式 (5) は，
• 次の yi+1 値は，もとの yi に，そこでの傾き f (xi , yi ) に x の増分 ∆x を乗じたものを加
えることにより求められる．
と言っている．イメージにすると，図 3 のようになる．
2
y
1
0
-1
0
2
4
6
8
10
x
図 1: 微分方程式
dy
dx
= sin(x) cos(x) − y cos(x) の方向場
y
1
0
-1
0
2
4
6
x
図 2: 方向場と解曲線
3
8
10
yi+1
yi
xi
∆x
xi+1
図 3: 方向場と微分方程式の解 (xi , yi ) と (xi+1 , yi+1 )
3 オイラー法の計算原理
初期条件を含めて数値計算により解くべき方程式は以下である．
dy
= f (x, y)
初期条件 y(a) = b
(6)
dx
この微分方程式の解を y = y(x) とすると， xi について微少量 δx を考えた場合のテイラー
展開は，
dy 1 d2 y 1 d3 y 2
yi+1 = y(xi + ∆x) = y(xi ) +
∆x
+
∆x
+
∆x3 + . . .
(7)
dx x=xi
2 dx2 x=xi
6 dx3 x=xi
となる．この式の右辺第 2 項は，式 (6) から計算できる．したがって，テイラー展開は，
yi+1 = yi + f (xi , yi )∆x + O(∆x2 )
(8)
と表すことができる．
オイラー法での数値計算では，計算の刻み幅 ∆x は十分に小さいとして，
yi+1 = yi + f (xi , yi )∆x
(9)
を計算する．
オイラー法をまとめると，以下に示すように微分方程式は差分方程式に近似できる．
yi+1 − yi
= f (xi , yi )
∆x
dy
xi+1 = xi + ∆x
dx = f (x, y)
(10)
⇒
y(a) = b
x0 = a
y = b
0
4
これから，オイラー法での数値計算の漸化式
yi+1 = yi + f (xi , yi )∆x
xi+1 = xi + ∆x
(11)
となる．初期値 (x0 , y0 ) が決まれば，(x1 , y1 ), (x2 , y2 ), · · · が同じ手続きで，芋づる式に計算
できる．通常，初期値 (a, b) は問題（プログラムの最初）で与えられる．
実際にプログラムを行うときは，for や while を用いて繰り返し計算を行う (芋づる式
の部分)．そして，計算結果の xi と yi は，配列 x[i] や y[i] に格納する．
x[0]=a;
y[0]=b;
dx = きざみ幅の計算
for(計算終了条件){
x[i+1]=x[i]+dx;
y[i+1]=y[i]+f(x[i],y[i])*dx;
}
この方法の計算のイメージは，図 4 のようになる．
y(x)
x0
x1
x2
x4
x
図 4: オイラー法．ある区間での y の変化 ∆y は，計算の始めの点の傾きに区間の幅 ∆x を
乗じて，求めている．
5
4 練習問題
• オイラー法を用いて常微分方程式を解くプログラムを作成する．プログ
ラムが妥当かどうかを知るために，まずは手計算で厳密解を計算できる
微分方程式を考える．以下に示す微分方程式を x の範囲を 0 ∼ 2 として
オイラー法により解析解を求めなさい．
dy
= 2x
dx
(12)
• 上の式の数値解析について，理論値と計算値を同じグラフに書きなさい．
• きざみ幅が 0.001 の場合の誤差のグラフを描き，それを考察しなさい．横
軸:x，縦軸:誤差．ただし，誤差は理論値と計算値の差とする．
• 次に，実際の物理問題をシミュレートしてみる．
物体の鉛直投げ上げを考えると物体の速度は，
dy
= v0 − gt
dt
(13)
となり，物体の高さ y(t)[m] はこれを t[s] について積分して
1
y = v0 t − gt2
2
(14)
となる．これらの式から計算すべき t[s] の範囲を見極めて，計算を行な
いなさい．ただし，初期条件は t0 = 0[s],y(0) = 0[m],v0 = 10[m/s] とする．
• この場合，誤差がどのようになったかを理由を含めて考察しなさい．
• オイラー法は微分方程式の数値計算方法で最も単純な方法であるが，解
の精度は一次であるため誤差が大きい．高精度に常微分方程式の数値計
算法をいくつか調べ，その計算原理の概要を示しなさい．
4.1
プログラムの穴埋め問題
リスト 1: オイラー法のプログラム（穴埋め）
1
2
3
4
5
6
7
8
9
10
11
# i n c l u d e < s t d i o . h>
# i n c l u d e <math . h>
# define
# define
# define
# define
# define
NSTEPS 100
NDIM 30000
X MAX 2 . 0
G 9.8
V0 1 0 . 0
//
//
//
//
//
計算ステップ数
用意する配列の大きさ
計算するの最大値 t
重力加速度 [m / s ˆ 2 ]
初速度 [m / s ]
double f ( double x , double y ) ;
double f r ( double x , double yr ) ;
6
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
/ / ============================================================
/ / main 関数
/ / ============================================================
i n t main ( v o i d ) {
d o u b l e y [NDIM] , x [NDIM ] , y r [NDIM ] , e r r o r [NDIM ] ;
/ / 計算結果を入れる配列
d o u b l e dx ;
int i ;
/ / カウンタ
FILE ∗ o u t ;
/ / 計算結果を保存
dx = ( d o u b l e )X MAX/ NSTEPS ;
y [0] = 0;
yr [0] = 0;
x [0] = 0;
/ / t の計算のきざみ幅
/ / y の初期化
/ / 理論値用の配列初期化
/ / t の初期化
/ / −−−−−−− オイラー法の計算 −−−−−−−−
f o r ( i =0; i <NSTEPS ; i ++){
/ / ==========================
/ / ここにオイラー法の計算を書く
/ / ==========================
}
/ / −−−−−−− 理論値の計算 −−−−−−−−−−−
f o r ( i =0; i <NSTEPS ; i ++){
/ / ==========================
/ / ここに理論値の計算を書く
/ / ==========================
}
/ / −−−−−−− 誤差の計算 −−−−−−−−−−−−−
f o r ( i =0; i <NSTEPS ; i ++){
/ / ==========================
/ / ここに誤差の計算を書く
/ / ==========================
}
/ / −−−−−−− 計算結果をファイルに保存 −−−−−−−−
o u t = f o p e n ( ” r e s u l t . t x t ” , ”w” ) ;
f o r ( i =0; i <=NSTEPS ; i ++){
f p r i n t f ( o u t , ” %20.15 f \ t %20.15 f \ t %20.15 f \ t %20.15 f \ t %d \ n ” , x [ i ] , y [ i ] , y r [ i ] , e r r o
}
f c l o s e ( out ) ;
return 0;
}
7
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
/ / ============================================================
/ / dy / d t = f ( x , y ) の f ( x , y ) を計算する関数
/ / ============================================================
double f ( double x , double y ) {
double ans ;
/ / ==========================
/ / ここに f ( x , y ) の計算を書く
/ / ==========================
return ans ;
}
/ / ============================================================
/ / 理論値を計算する関数
/ / ============================================================
double f r ( double x , double yr ) {
double ans ;
/ / ========================
/ / ここに理論値の計算を書く
/ / ========================
return ans ;
}
8
© Copyright 2024 ExpyDoc
