● 講義資料

● 講義資料

▼ 講義予定

• 収束の加速

• １階常微分方程式の初期値問題の数値解法

● 前回の講義のまとめ

▼ ニュートン法について

• √

2に代表される,代数方程式の解の近似値を求める方法として,以下のニュートン法が広く 知られている.

• 区間[a, b]上で定義された実数値C²-級関数f(x)が, [a, b]上で単調増加かつ下に凸であり,

f(a)<0,f(b)>0を満たしていると仮定する. この時,x0=b, xn+1 =xn−f(xn)/f^′(xn) により定義された数列{xⁿ} は,単調減少かつ下に有界であり, [a, b]上のf(x) = 0のただひ とつの解αに収束する. さらに,ある定数M >0 が存在して,

|xⁿ+1−α| ≤M|xⁿ−α|² が成り立つ.

• ニュートン法は以下の「逐次近似」の一例と考えられる.

• F:R−→Rを C²-級関数とする. この時,

xn+1=F(xⁿ)

によって数列{xⁿ}を定義して,写像F の固定点(α=F(α)を満たす点)を近似する方法を 逐次近似と呼ぶ.

• 関数f(x)に対するニュートン法xn+1=xⁿ−f(xⁿ)/f^′(xⁿ)は, F(x) =x−f(x)/f^′(x)と 定義した逐次近似である.

• Rのある近傍U と正定数L <1が存在して,任意のx,y∈U に対して

|F(x)−F(y)| ≤L|x−y| を満たすとき,F はU 上で縮小写像であると言う.

• C²-級関数F:R−→Rが,U ⊂R上で縮小写像であるとき,F は U にただひとつの固定点 を持つ（縮小写像の原理）.

証明は,x0∈U,xⁿ+1=F(xⁿ)で定義した点列{xⁿ}がコーシー列となることを示し,その 収束先が固定点となることを示せばよい.

• C²-級関数F: R−→Rに対する逐次近似xⁿ+1 =F(xⁿ)の固定点αがあらかじめわかって いる時を考える. この時,|F^′(α)|<1が成り立てば,αの近傍でF は縮小写像となる. 特に, αの十分近くに初期値x0 をとった逐次近似xⁿ+1=F(xⁿ)は,固定点αに収束する.

• 逐次近似xn+1=F(xⁿ)がαに収束すると仮定する. この時,ǫⁿ=|xⁿ−α|は ǫn+1≤ |F^′(α)|ǫⁿ+1

2|F^′′(α)|ǫ²n+O(ǫ³n) をみたす. 証明は,αでF を３次までテイラー展開して計算すればよい.

• ニュートン法の収束の速さは,f(x) = 0の解が単根のとき, ǫn+1≤Cǫ²n

をみたし,m 重根の時

ǫⁿ+1≤ m−1 m ǫⁿ をみたす.

証明は,F(x) =x−f(x)/f^′(x)と考え,上の結果を適用すればよい.

• この結果から,ニュートン法は,解が単根の時と重根の時とでは,その収束の様子が質的に異 ることがわかる.

• さらに,ニュートン法の収束の判定には,次の結果を用いることができる.

1. ある値αの近似列{xⁿ}が|α−xn+1| ≤C|α−xⁿ|²を満たしていると仮定する. この時, 任意の0< C <1 に対して,|xn+1−xⁿ|< Cǫ|xⁿ|が成り立つならば,|α−xⁿ|< ǫ|α| が成り立つ.

2. ある値 α の近似列 {xⁿ} が, ある 0 < L < 1 が存在して, |α−xn+1| ≤ L|α−xⁿ| を満たしていると仮定する. この時, |xⁿ+1−xⁿ| < (1−L)ǫ|xⁿ| が成り立つならば,

|α−xⁿ|< ǫ|α|が成り立つ.

• たとえば,ニュートン法の単根の場合には,収束の判定条件として,「|xⁿ+1−xⁿ|> ǫ|xⁿ|が成 り立つ間,繰り返しを続ける」と書けば,そのループから脱出した時点では,|xn+1−xⁿ|< ǫ|xⁿ| が成り立っているので,近似値xn+1 は解αを相対誤差ǫで近似していると考えてよいこと がわかる.

▼ その他の求根アルゴリズム

• 初回に解説した二分法と今回のニュートン法のほかに,簡単な求根アルゴリズムとしては,以 下のものが知られている.

以下では, 関数f は区間[a, b]上で定義された, 単調増加関数であり, f(a)<0,f(b)>0 と 仮定する. また,区間[a, b]におけるf(x) = 0の一意的な解をαとする.

• 「はさみうち法」などと呼ばれる方法.

– 二分法では,繰り返しの際にcⁿ= (aⁿ+bⁿ)/2 におけるf(cⁿ)の値の正負によって,区 間を取り替えていた.

– これを(aⁿ, f(aⁿ)), (bⁿ, f(bⁿ))を結ぶ直線と x軸との交点をcⁿ とおいて, f(cⁿ)の値 の正負によって区間を取り替える方法である.

– 二分法と同じく,「一次収束」する.

– 二分法と同じく, 任意のnに対して,α∈[aⁿ, bⁿ] が成り立つ. したがって,求めた近似 値の精度保証が容易である.

• 「割線法」

– 関数f が区間[a, b]で下に凸であると仮定できるとき,ニュートン法で「接線」を取った

代わりに, (xⁿ, f(xⁿ)), (xⁿ+1, f(xⁿ+1))を通る直線とx軸の交点をxⁿ+2 とおく方法.

– 収束の次数は(1 +√

5)/2となる. （ニュートン法よりも収束が遅い）

– 関数の微分f^′ が複雑な式になっている場合や, f^′ を数値的にしか求めることができな

● 講義資料

▼ 収束の加速

次の図は,単位円に接する正多角形の周長による円周率の近似計算を, Romberg加速とAitken 加速を使って加速した例である.

1e-16 1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1

10 100 1000 10000 100000 1e+06 1e+07 1e+08 1e+09

original 1 step 2 step Aitken

次の図は, 単位円に接する正多角形の周長による円周率の近似計算を, Romberg加速を用いて計算 している例であり,λの値を正しい値から取り替えて計算した例である.

1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01

10 100 1000 10000 100000 1e+06 1e+07 1e+08 1e+09

original 0.25 0.245 0.255 0.24 0.26 0.20 0.30

次の図は, arctan(x)の近似値のテイラー展開による近似計算を, Aitken 加速を使って加速した例 である.

1e-16 1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1

1 10 100 1000 10000

original x=1.00 Aitken x=1.00 original x = 0.99 Aitken x = 0.99 original x = 0.90 Aitken x = 0.90

▼ 題材となる常微分方程式の例

講義中での説明・実習問題の題材となる常微分方程式の（初期値問題の）例としては,以下のも のがある. 物理定数などはすべて1 ととってある.

1. (指数関数)

x(t)^′ =λx(t),λは定数.

2. (ロジスティック方程式) x(t)^′ = (1−x(t)²).

3. (変数分離形の別の例) x(t)^′ = sin(x(t)).

4. (単振動) x(t)^′′=−x(t).

5. (２階線形常微分方程式)

ax(t)^′′+bx(t)^′+cx(t) =f(t). a,b,cは定数, a6= 0,f はなめらかな関数.

6. (単振り子)

x(t)^′′=−sin(x(t)).

7. (van der Pol 方程式)

x(t)^′′−(1−x(t)²)x(t)^′+x(t) = 0.

8. (中心力場の運動) X(t)^′′=− X(t)

|X(t)|³. (X(t)∈R³).

9. (二重振り子)

x^′′1 = −2 sin(x1)−x^′2²sin(x1−x2)−x^′1²cos(x1−x2) sin(x1−x2) + cos(x1−x2) sin(x2

2−cos(x1−x2)²

x^′′2 = 2 cos(x1−x2) sin(x1) + 2x^′1²sin(x1−x2) +x^′2²cos(x1−x2) sin(x1−x2)−2 sin(x2)

▼ 前進オイラー法による数値計算

以下の図は,前進オイラー法で種々の常微分方程式の初期値問題の数値解を計算した例である.

（特に書いていない限りh= 0.01である.）

0 0.5 1 1.5 2 2.5 3

0 0.2 0.4 0.6 0.8 1

x’ = x

0 0.2 0.4 0.6 0.8 1 1.2 1.4

0 0.2 0.4 0.6 0.8 1

x’ = -x

x^′(t) =x(t),x(0) = 1 x^′(t) =−x(t),x(0) = 1

-0.2 0 0.2 0.4 0.6 0.8 1 1.2

0 0.5 1 1.5 2 2.5 3 3.5 4 x’ = (1-x^2)

0 0.5 1 1.5 2 2.5 3 3.5 4

0 2 4 6 8 10

x’ = sin(x)

x^′(t) = (1−x(t)²),x(0) = 0 x^′(t) = sin(x(t)),x(0) = 1

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

0 2 4 6 8 10 12 14 16 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 x’’ = -x

x^′′(t) =−x(t),x(0) = 1,x^′(0) = 0

-3 -2 -1 0 1 2 3

0 2 4 6 8 10 12 14 16

x’’ = -sin(x)

-2 -1 0 1 2

-2 -1 0 1 2

x’’ = -sin(x)

x^′′(t) =−sin(x(t)),x(0) = 2,x^′(0) = 0

-3 -2 -1 0 1 2 3

0 5 10 15 20 25

x’’ = (1-x^2)x’ - x

-3 -2 -1 0 1 2 3

-3 -2 -1 0 1 2 3 x’’ = (1-x^2)x’ - x

x^′′(t) = (1−x(t)²)x^′(t)−x(t),x(0) = 1.5,x^′(0) = 1

-10 -5 0 5 10

-10 -5 0 5 10

X’’ = -X/|X|^3

X^′′(t) =−X(t)/|X(t)|³,X(0) = (1,1,0),X^′(0) = (1,0,0)

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2

0 2 4 6 8 10

Double Pendulum

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 Double Pendulum

二重振り子,

● 実習内容

（以下で「★」の数は推奨の程度を示します. 多いほど推奨の度合いが大きい. また,「†」は 難易度またはマニアックな程度を表します. 多いほど難しいかマニアックかです.）

1. （★★★）単位円に接する多角形の周長による円周率の近似計算をRomberg 加速によって 計算しなさい. （この資料の最初の図のような図を書きなさい）

2. （★★★）単位円に接する多角形の周長による円周率の近似計算を２段階のRomberg加速 によって計算しなさい. （この資料の最初の図のような図を書きなさい）

3. （★★★）単位円に接する多角形の周長による円周率の近似計算をAitlen加速によって計算 しなさい. （この資料の最初の図のような図を書きなさい）

4. （★★★）テイラー展開によるarctan(1), arctan(0.99)の近似値の計算をAitlen加速によっ て計算しなさい. （この資料の２つめの図のような図を書きなさい）

5. （★★）f(x) = (x²−2)^k に対するニュートン法による √

2 の近似計算を, Aitken 加速に よって計算しなさい.

6. （★★†）f(x) = (x²−2)^k に対するニュートン法による √

2 の近似計算を, k≥2 の時に

Romberg加速によって計算しなさい.

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（出てきた結果が正しいか否かを判断するのが難しい）