解答例(PDF)

コンピュータ工学基礎演習 2014
10 常微分方程式を解く・解答例
機能創成専攻生体工学領域井村誠孝
[email protected]
課題ex10-1 [必須]
講義資料: http://bit.ly/bpe-ecp14
 次の微分方程式をオイラー法を用いて解き，t=2の
ときのyの値y (2)を小数点以下6桁まで求めよ．時
間刻みは ∆t =10−4 とせよ．
y′ =− y + 2t exp(−t )
y (0) = 0
 解答
y (2) = 0.541350
2014/6/30
Exercises in Computer Programming
2 / 16
講義資料: http://bit.ly/bpe-ecp14
解のグラフ
2014/6/30
Exercises in Computer Programming
3 / 16
講義資料: http://bit.ly/bpe-ecp14
確認: 小さい値の表記
 C言語では a × 10b を aeb と表記できます(指数表
記)．
例
 12.3456 (=1.23456×101) → 1.23456e1
 あまり意味はないですが、123456e-4 と書くことも可能．
 1.0×10-1 → 1e-1
 0が沢山付く場合の表記に適している．
 0.00000002 → 2e-8
2014/6/30
Exercises in Computer Programming
4 / 16
課題ex10-2 [発展]
講義資料: http://bit.ly/bpe-ecp14
 前問と同じ常微分方程式の初期値問題について引
き続き考える．y (2)の数値解を，時間刻みを変え
て求め，時刻の増分 ∆t と，数値解と解析解との誤
差 ∆yとの関係について調査せよ．
 解析解は
=
y t 2 exp(−t )
2014/6/30
Exercises in Computer Programming
5 / 16
講義資料: http://bit.ly/bpe-ecp14
解のグラフ
2014/6/30
Exercises in Computer Programming
6 / 16
講義資料: http://bit.ly/bpe-ecp14
t=2付近の拡大
2014/6/30
Exercises in Computer Programming
7 / 16
講義資料: http://bit.ly/bpe-ecp14
時間の増分と誤差の関係
誤差大
刻みは小さければ
小さい方がよいとは
限らない．
DTを1e-1(10-1)から1e-9(10-9)
まで変え，真値との差を調べる．
誤差小
刻み小
2014/6/30
刻み大
Exercises in Computer Programming
8 / 16
講義資料: http://bit.ly/bpe-ecp14
時間刻みを小さくすると誤差が増える
 DTが非常に小さいことがどこで影響しているか?
y += dydt(t, y) * DT;
 yの大きさ: 1(100)のオーダー
 dydt(導関数の値)の大きさ: 1(100)のオーダー
 DT: 10-9のオーダー
 dydt*DT: 10-9のオーダー(有効数字は十分ある)
→doubleの精度は10進数で15桁程度であるから，
dydt*DT の上 6 桁程度しか y には加えられない．
dydt*DTの実質的な有効桁数は6桁(情報落ち)．
2014/6/30
Exercises in Computer Programming
9 / 16
課題ex10-3 [必須]
講義資料: http://bit.ly/bpe-ecp14
 次の微分方程式をオイラー法とルンゲ・クッタ法
を用いてそれぞれ解き， t=1 のときのyの値y (1)を
求めよ．同じ ∆t のときの，誤差について調べよ．
y′ = sin t
y (0) = 1
2014/6/30
Exercises in Computer Programming
10 / 16
講義資料: http://bit.ly/bpe-ecp14
解析解を求める
両辺をtで積分して，
dy
= sin t
dt
⇔ y=
− cos t + C
初期条件t=0でy=1より，
1 = C −1 ⇔ C = 2
よって
2014/6/30
y= 2 − cos t
2 − cos1 =
1.45969769413186...
y (1) =
Exercises in Computer Programming
11 / 16
講義資料: http://bit.ly/bpe-ecp14
ルンゲ・クッタ法の実装
オイラー法
y += dydt( t, y ) * DT;
ルンゲ・クッタ法
double k1, k2, k3, k4, dy;
k1 = dydt(t, y);
k2 = dydt(t + DT/2, y + k1*DT/2);
k3 = dydt(t + DT/2, y + k2*DT/2);
k4 = dydt(t + DT, y + k3*DT );
dy = (k1 + 2*k2 + 2*k3 + k4) / 6;
y += dy * DT;
2014/6/30
Exercises in Computer Programming
12 / 16
講義資料: http://bit.ly/bpe-ecp14
誤差の比較
DT
オイラー法
ルンゲ・クッタ法
10-1
10-2
10-3
4.25×10-2
4.21×10-3
4.21×10-4
1.60×10-8
1.60×10-12
<10-15
ルンゲ・クッタ法の方が圧倒的に誤差が小さい．
オイラー法: 誤差～DT
ルンゲ・クッタ法: 誤差～DT4
2014/6/30
Exercises in Computer Programming
13 / 16
課題ex10-4 [発展]
講義資料: http://bit.ly/bpe-ecp14
 前問の常微分方程式の数値解について，誤差が同
程度になる場合(つまり時間刻みが異なる)の，計
算回数(導関数の評価回数)について，どちらが有
利か考察せよ．
2014/6/30
Exercises in Computer Programming
14 / 16
講義資料: http://bit.ly/bpe-ecp14
計算回数の評価...の前に
 1.0 / DT すなわち繰り返し回数が整数にならない
場合，tにDTを整数回加える限りはt=1.0にはなら
ず、y(1)が求められません．
 例: DT=0.07 とすると，1.0/DT=14.28... なので14回の繰
り返しが行われる．このとき最終的なtの値は
t=0.07×14=0.98となるため，誤差が生じる．
 対処法
 繰り返し数を1回増やす．
 現在のtにDTを加えると1より大きくなる場合には，DTで
はなく1-tを加える(刻み幅を最後だけ可変に)．
2014/6/30
Exercises in Computer Programming
15 / 16
講義資料: http://bit.ly/bpe-ecp14
計算回数の評価
 1.0×10-8程度の誤差になるDTを調べる．
 オイラー法: DT=2.50×10-8
 ルンゲ・クッタ法: DT=8.96×10-2
 繰り返し回数は1.0/DT程度
 オイラー法: 4.0×107回
 ルンゲ・クッタ法: 12回
 dydt()が呼び出された回数
 オイラー法: 4.0×107回
計算効率: 約83万倍異なる．
 ルンゲ・クッタ法: 48回
(この事例に限れば)
2014/6/30
Exercises in Computer Programming
16 / 16

Download Report