講義資料前回の講義のまとめ

2014 年度・後期・数理解析・計算機数学３・第６回
1
● 講義資料
▼ 講義予定
• 収束の加速
• １階常微分方程式の初期値問題の数値解法
● 前回の講義のまとめ
▼ ニュートン法について
•
√
2 に代表される, 代数方程式の解の近似値を求める方法として, 以下のニュートン法が広く
知られている.
• 区間 [a, b] 上で定義された実数値 C 2 -級関数 f (x) が, [a, b] 上で単調増加かつ下に凸であり,
f (a) < 0, f (b) > 0 を満たしていると仮定する. この時, x0 = b, xn+1 = xn − f (xn )/f ′ (xn )
により定義された数列 {xn } は, 単調減少かつ下に有界であり, [a, b] 上の f (x) = 0 のただひ
とつの解 α に収束する. さらに, ある定数 M > 0 が存在して,
|xn+1 − α| ≤ M |xn − α|2
が成り立つ.
• ニュートン法は以下の「逐次近似」の一例と考えられる.
• F : R −→ R を C 2 -級関数とする. この時,
xn+1 = F (xn )
によって数列 {xn } を定義して, 写像 F の固定点 (α = F (α) を満たす点) を近似する方法を
逐次近似と呼ぶ.
• 関数 f (x) に対するニュートン法 xn+1 = xn − f (xn )/f ′ (xn ) は, F (x) = x − f (x)/f ′ (x) と
定義した逐次近似である.
• R のある近傍 U と正定数 L < 1 が存在して, 任意の x, y ∈ U に対して
|F (x) − F (y)| ≤ L|x − y|
を満たすとき, F は U 上で縮小写像であると言う.
• C 2 -級関数 F : R −→ R が, U ⊂ R 上で縮小写像であるとき, F は U にただひとつの固定点
を持つ（縮小写像の原理）.
証明は, x0 ∈ U , xn+1 = F (xn ) で定義した点列 {xn } がコーシー列となることを示し, その
収束先が固定点となることを示せばよい.
• C 2 -級関数 F : R −→ R に対する逐次近似 xn+1 = F (xn ) の固定点 α があらかじめわかって
いる時を考える. この時, |F ′ (α)| < 1 が成り立てば, α の近傍で F は縮小写像となる. 特に,
α の十分近くに初期値 x0 をとった逐次近似 xn+1 = F (xn ) は, 固定点 α に収束する.
証明は, x の近傍で平均値の定理を用いればよい.
Nov. 05, 2014, Version: 1.0.1
[email protected]
2014 年度・後期・数理解析・計算機数学３・第６回
2
• 逐次近似 xn+1 = F (xn ) が α に収束すると仮定する. この時, ǫn = |xn − α| は
1
ǫn+1 ≤ |F ′ (α)|ǫn + |F ′′ (α)|ǫ2n + O(ǫ3n )
2
をみたす. 証明は, α で F を３次までテイラー展開して計算すればよい.
• ニュートン法の収束の速さは, f (x) = 0 の解が単根のとき,
ǫn+1 ≤ Cǫ2n
をみたし, m 重根の時
ǫn+1 ≤
m−1
ǫn
m
をみたす.
証明は, F (x) = x − f (x)/f ′ (x) と考え, 上の結果を適用すればよい.
• この結果から, ニュートン法は, 解が単根の時と重根の時とでは, その収束の様子が質的に異
ることがわかる.
• さらに, ニュートン法の収束の判定には, 次の結果を用いることができる.
1. ある値 α の近似列 {xn } が |α − xn+1 | ≤ C|α − xn |2 を満たしていると仮定する. この時,
任意の 0 < C < 1 に対して, |xn+1 − xn | < Cǫ|xn | が成り立つならば, |α − xn | < ǫ|α|
が成り立つ.
2. ある値 α の近似列 {xn } が, ある 0 < L < 1 が存在して, |α − xn+1 | ≤ L|α − xn |
を満たしていると仮定する. この時, |xn+1 − xn | < (1 − L)ǫ|xn | が成り立つならば,
|α − xn | < ǫ|α| が成り立つ.
• たとえば, ニュートン法の単根の場合には, 収束の判定条件として, 「|xn+1 − xn | > ǫ|xn | が成
り立つ間, 繰り返しを続ける」と書けば, そのループから脱出した時点では, |xn+1 −xn | < ǫ|xn |
が成り立っているので, 近似値 xn+1 は解 α を相対誤差 ǫ で近似していると考えてよいこと
がわかる.
Nov. 05, 2014, Version: 1.0.1
[email protected]
2014 年度・後期・数理解析・計算機数学３・第６回
3
▼ その他の求根アルゴリズム
• 初回に解説した二分法と今回のニュートン法のほかに, 簡単な求根アルゴリズムとしては, 以
下のものが知られている.
以下では, 関数 f は区間 [a, b] 上で定義された, 単調増加関数であり, f (a) < 0, f (b) > 0 と
仮定する. また, 区間 [a, b] における f (x) = 0 の一意的な解を α とする.
• 「はさみうち法」などと呼ばれる方法.
– 二分法では, 繰り返しの際に cn = (an + bn )/2 における f (cn ) の値の正負によって, 区
間を取り替えていた.
– これを (an , f (an )), (bn , f (bn )) を結ぶ直線と x 軸との交点を cn とおいて, f (cn ) の値
の正負によって区間を取り替える方法である.
– 二分法と同じく, 「一次収束」する.
– 二分法と同じく, 任意の n に対して, α ∈ [an , bn ] が成り立つ. したがって, 求めた近似
値の精度保証が容易である.
• 「割線法」
– 関数 f が区間 [a, b] で下に凸であると仮定できるとき, ニュートン法で「接線」を取った
代わりに, (xn , f (xn )), (xn+1 , f (xn+1 )) を通る直線と x 軸の交点を xn+2 とおく方法.
√
– 収束の次数は (1 + 5)/2 となる. （ニュートン法よりも収束が遅い）
– 関数の微分 f ′ が複雑な式になっている場合や, f ′ を数値的にしか求めることができな
い場合には, ニュートン法よりも誤差を小さくできる可能性がある.
Nov. 05, 2014, Version: 1.0.1
[email protected]
2014 年度・後期・数理解析・計算機数学３・第６回
4
● 講義資料
▼ 収束の加速
次の図は, 単位円に接する正多角形の周長による円周率の近似計算を, Romberg 加速と Aitken
加速を使って加速した例である.
1
original
1 step
2 step
Aitken
0.01
0.0001
1e-06
1e-08
1e-10
1e-12
1e-14
1e-16
10
100
1000
10000
100000
1e+06
1e+07
1e+08
1e+09
次の図は, 単位円に接する正多角形の周長による円周率の近似計算を, Romberg 加速を用いて計算
している例であり, λ の値を正しい値から取り替えて計算した例である.
original
0.25
0.245
0.255
0.24
0.26
0.20
0.30
0.01
0.0001
1e-06
1e-08
1e-10
1e-12
1e-14
10
100
1000
10000
100000
1e+06
1e+07
1e+08
1e+09
次の図は, arctan(x) の近似値のテイラー展開による近似計算を, Aitken 加速を使って加速した例
である.
1
original x=1.00
Aitken x=1.00
original x = 0.99
Aitken x = 0.99
original x = 0.90
Aitken x = 0.90
0.01
0.0001
1e-06
1e-08
1e-10
1e-12
1e-14
1e-16
1
Nov. 05, 2014, Version: 1.0.1
10
100
1000
10000
[email protected]
2014 年度・後期・数理解析・計算機数学３・第６回
5
▼ 題材となる常微分方程式の例
講義中での説明・実習問題の題材となる常微分方程式の（初期値問題の）例としては, 以下のも
のがある. 物理定数などはすべて 1 ととってある.
1. (指数関数)
x(t)′ = λx(t), λ は定数.
2. (ロジスティック方程式)
x(t)′ = (1 − x(t)2 ).
3. (変数分離形の別の例)
x(t)′ = sin(x(t)).
4. (単振動)
x(t)′′ = −x(t).
5. (２階線形常微分方程式)
ax(t)′′ + bx(t)′ + cx(t) = f (t). a, b, c は定数, a 6= 0, f はなめらかな関数.
6. (単振り子)
x(t)′′ = − sin(x(t)).
7. (van der Pol 方程式)
x(t)′′ − (1 − x(t)2 )x(t)′ + x(t) = 0.
8. (中心力場の運動)
X(t)
X(t)′′ = −
. (X(t) ∈ R3 ).
|X(t)|3
9. (二重振り子)
x′′1 =
x′′2 =
′2
−2 sin(x1 ) − x′2
2 sin(x1 − x2 ) − x1 cos(x1 − x2 ) sin(x1 − x2 ) + cos(x1 − x2 ) sin(x2
2 − cos(x1 − x2 )2
′2
2 cos(x1 − x2 ) sin(x1 ) + 2x′2
1 sin(x1 − x2 ) + x2 cos(x1 − x2 ) sin(x1 − x2 ) − 2 sin(x2 )
2 − cos(x1 − x2 )2
Nov. 05, 2014, Version: 1.0.1
[email protected]
2014 年度・後期・数理解析・計算機数学３・第６回
6
▼ 前進オイラー法による数値計算
以下の図は, 前進オイラー法で種々の常微分方程式の初期値問題の数値解を計算した例である.
（特に書いていない限り h = 0.01 である.）
x’ = x
x’ = -x
3
1.4
1.2
1
0.8
0.6
0.4
0.2
0
2.5
2
1.5
1
0.5
0
0
0.2
0.4
0.6
0.8
1
0
x′ (t) = x(t), x(0) = 1
0.2
0.4
0.6
0.8
1
x′ (t) = −x(t), x(0) = 1
x’ = (1-x^2)
x’ = sin(x)
1.2
4
3.5
3
2.5
2
1.5
1
0.5
0
1
0.8
0.6
0.4
0.2
0
-0.2
0
0.5
1
1.5
2
2.5
3
3.5
4
x′ (t) = (1 − x(t)2 ), x(0) = 0
0
2
4
6
8
10
x′ (t) = sin(x(t)), x(0) = 1
x’’ = -x
x’’ = -x
2
1.5
1
0.5
0
-0.5
-1
-1.5
-2
2
1.5
1
0.5
0
-0.5
-1
-1.5
-2
0
2
4
6
8
10
12
14
16
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
x′′ (t) = −x(t), x(0) = 1, x′ (0) = 0
Nov. 05, 2014, Version: 1.0.1
[email protected]
2014 年度・後期・数理解析・計算機数学３・第６回
7
x’’ = -sin(x)
x’’ = -sin(x)
3
2
2
1
1
0
0
-1
-1
-2
-2
-3
0
2
4
6
8
10
12
14
16
-2
-1
0
1
2
x′′ (t) = − sin(x(t)), x(0) = 2, x′ (0) = 0
x’’ = (1-x^2)x’ - x
x’’ = (1-x^2)x’ - x
3
3
2
2
1
1
0
0
-1
-1
-2
-2
-3
-3
0
5
10
15
20
25
-3 -2 -1
0
1
2
3
x′′ (t) = (1 − x(t)2 )x′ (t) − x(t), x(0) = 1.5, x′ (0) = 1
X’’ = -X/|X|^3
10
5
0
-5
-10
-10
-5
0
5
10
X ′′ (t) = −X(t)/|X(t)|3 , X(0) = (1, 1, 0), X ′ (0) = (1, 0, 0)
Double Pendulum
Double Pendulum
2
1.5
1
0.5
0
-0.5
-1
-1.5
-2
-2.5
1.5
1
0.5
0
-0.5
-1
-1.5
-2
-2.5
0
2
4
6
8
10
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
二重振り子, x1 (0) = 2, x2 (0) = 0, x′1 (0) = x′2 (0) = 0
Nov. 05, 2014, Version: 1.0.1
[email protected]
2014 年度・後期・数理解析・計算機数学３・第６回
8
● 実習内容
（以下で「★」の数は推奨の程度を示します. 多いほど推奨の度合いが大きい. また, 「†」は
難易度またはマニアックな程度を表します. 多いほど難しいかマニアックかです.）
1. （★★★）単位円に接する多角形の周長による円周率の近似計算を Romberg 加速によって
計算しなさい. （この資料の最初の図のような図を書きなさい）
2. （★★★）単位円に接する多角形の周長による円周率の近似計算を２段階の Romberg 加速
によって計算しなさい. （この資料の最初の図のような図を書きなさい）
3. （★★★）単位円に接する多角形の周長による円周率の近似計算を Aitlen 加速によって計算
しなさい. （この資料の最初の図のような図を書きなさい）
4. （★★★）テイラー展開による arctan(1), arctan(0.99) の近似値の計算を Aitlen 加速によっ
て計算しなさい. （この資料の２つめの図のような図を書きなさい）
5. （★★）f (x) = (x2 − 2)k に対するニュートン法による
よって計算しなさい.
√
6. （★★†）f (x) = (x2 − 2)k に対するニュートン法による
Romberg 加速によって計算しなさい.
2 の近似計算を, Aitken 加速に
√
2 の近似計算を, k ≥ 2 の時に
7. 上記 1 から 9 の常微分方程式の初期値問題の数値解を前進オイラー法で計算しなさい.
• 特に指定していない限り x(t) ∈ R である.
• 刻み幅は（とりあえず） h = 0.01 とする. (h を取り替えて計算してみるとよい)
• 初期時刻 t = 0 から, 解の様子が良くわかる程度の t で計算すること. （あまり短くて
も解の様子がわからないし, あまり長くても余計にわけがわからなくなることがある）
• 初期条件は, 適切にとること.
• 1 次元 2 階常微分方程式に関しては, 相空間での解の様子の図もつくること.
• 自由度が高いものについては, 解の様子がわかる図をつくること.
• 厳密解がわかるものについては, 厳密解との比較をするとよい.
• 1, 5 は適切に定数などを設定して解くとよい.
なお, 難易度と推奨の度合いは以下の通り.
• （★★★） 1, 2, 3, 4, 5, 6, 7
• （★★★†） 8 （ちょっと面倒だから）
• （★★★†††） 9 （出てきた結果が正しいか否かを判断するのが難しい）
Nov. 05, 2014, Version: 1.0.1
[email protected]