● 前回の講義のまとめ

● 前回の講義のまとめ

★ ニュートン法と逐次近似

•

ニュートン法は, なめらかな関数

f:R−→R

の解を求めるために用いられる方法であり, 通常, 以下 の定理に基づく方法と考えられる.

なめらかな関数

f:R−→R

に対して,

f(a)<0, f(b)>0,f^′(x)>0 on [a, b],f^′′(x)<0 on (a, b)

をみたすとき,

xn+1=xn− f(xn)

f^′(xn), x0=b

によって定められる数列

{xn}

は,

f(x) = 0

の区間

(a, b)

での唯一の解

α

に収束する.

•

したがって,

f

の零点を求めるためにニュートン法を適用する場合, 零点

α

の近傍（片側の近傍だけ でもよい）で

f

の単調性と凸性が保証される必要がある. （逆に言えば, それらが保証されるような 近傍に初期値

x0

をとる必要がある）

•

より一般に, なめらかな関数

φ:R−→R

が, 次の条件を満たすと仮定する.

ある正の定数

L <1

が存在して, 任意の

x,y∈(a, b)

に対して

|φ(x)−φ(y)| ≤L|x−y| (1)

を満たす. このとき,

φ

は区間

(a, b)

上にただ一つの固定点

α(φ(α) =α)

を持つ. 特に,

xn+1=φ(xn)

によって定められる数列

{xn}

は

α

に収束する.

–

特に,

|φ^′(α)| <1

が成り立つとき,

α

を含む区間

(a, b)

を十分小さく取れば, (1) が成り立つ.

（講義では, この形で説明をした.）

–

このとき,

|xn+1−α| ≤ |φ^′(α)| |xn−α|+|φ^′′(α)|

2 |xn−α|²+O(|xn−α|³)

が成り立つ.

–

ニュートン法は,

φ(x) =x−_f^f(x)′(x)

と置いた逐次近似と考えられ,

φ^′(x) =−f(x)f^′′(x)

(f^′(x))² , φ^′′(x) =f^′′(x)

f^′(x) −2f(x)(f^′′(x))² (f^′(x))³

が成り立つ.

–

特に, 解が単根

(f^′(α)6= 0)

のとき,

φ^′(α) = 0, φ^′′(α) = f^′′(α) f^′(α)

となり,

|xn+1−α| ≤M|xn−α|²+O(|xn−α|³)

となり, 二次収束することがわかる.

–

また, 解が

m

重根（f

^′(α) =· · ·=f^(m−1)(α) = 0,f^(m)(α)6= 0

のとき）

|xn+1−α| ≤ m−1

m |xn−α|

となり, 一次収束する.

•

このことから, 求めたい解

α

の十分近くに初期値をとったニュートン法は, 解

α

に収束することがわ かる. しかし, 一般には「どのくらい近くに初期値をとればよいか」がわからない. そのため, ニュー トン法を利用する場合には, 以下の情報にしたがって適切に初期値をとる必要がある.

–

解のおおよその値. （たとえば, 正と負の２つの解を持つことを確かめ, 正の解を求めるなど）

–

解の重複度と単調性と凸性. （たとえば, 単根の場合, 単調増加かつ下に凸の範囲を求め, その範 囲内に初期値をとるなど）

•

ニュートン法では, 「解

α

の近似値を, 相対誤差

δ

で求める」ために, いつ繰り返しを停止させるかは明 らかではない. すなわち, 真の値

α

と「近似値」

xn

を直接比較できないため, 評価式

“|xn−α| ≤δ|α|”

を利用することができない. そのため, 上の評価式（二次収束の式, または, 一次収束の式）を利用す ることにより, 次のように評価することができる.

★ 単根の時

|xn+1−xn| ≤ ^δ₂|xn|

が成り立てば,

|xn−α| ≤δ|α|

が成り立つ.

★

m

重根の時

|xn+1−xn| ≤ _2m^δ |xn|

が成り立てば,

|xn−α| ≤δ|α|

が成り立つ.

★ 収束の加速

•

数列

{an}

が極限値

α

を持ち,

an=α+c1λⁿ₁+c2λⁿ₂+· · · , 1> λ1> λ2>· · ·>0

という振る舞いをしているとき,

a⁽¹⁾_n = an−λ1an−1

1−λ1

で定義された数列

{a⁽¹⁾n }

も同じ極限に収束し,

|a⁽¹⁾_n −α|=c^′₂λ²₂+· · ·

をみたす. これを

Romberg

加速と呼ぶ. （Romberg-Richardson 加速,

Richardson

加速,

Richiardson

外挿と呼ぶこともある）

•

数列

{an}

が極限値

α

を持ち,

(an+1−α) = (A+ǫn)(an−α), ǫn→0, |A|<1

という振る舞いをしているとき,

a⁽¹⁾_n =an− (an−an+1)² an+2−2an+1+an

で定義された数列

{a⁽¹⁾n }

も同じ極限に収束し,

a⁽¹⁾n −α an−α →0

の意味で加速される. これを

Aitken

加速と呼ぶ.

•

いずれの加速方法でも, 加速された数列

{a⁽¹⁾n }

に対して, 再び加速を施すことが可能である. Romberg 加速を用いるときには, 加速された数列

{a⁽¹⁾n }

は

a⁽¹⁾_n =α+c⁽¹⁾₂ λⁿ₂+c⁽¹⁾₃ λⁿ₃ +· · ·, 1> λ2> λ3>· · ·>0

となるので, 加速係数としては

λ2

を用いれば良い.

• Romberg

加速は「収束の主要項」

λ1

の値が分からなければ加速することができないが, Aitken 加

速は, それが分からなくても加速することが可能である. 特に１次より速い収束をしている列に対し ても加速が可能である. しかし, 一方では

Aitken

加速は「収束しない数列」を収束させてしまうこ とがあるので注意が必要である.

•

半径１の円に内接する正

n

角形の周長を

2ℓn

とおいた, 列

{ℓn}

は,

ℓk−π= π³

6 k⁻²− π⁵

120k⁻⁴+O(k⁻⁶)

をみたすので,

an=ℓn02ⁿ

とおくと,

an−π= π³

6n²₀(1/4)ⁿ− π⁵

120n⁴₀(1/16)ⁿ+· · ·

をみたす. したがって,

λ1= 1/4

として

Romberg

加速によって収束を加速することができる.

•

ニュートン法

(m

重根) の場合には, 各ステップでの誤差

ǫn=|an−α|

は

ǫn+1= ^m−_m¹ǫn

となるので,

an =α+C

m−1

m n

+o(

m−1

m n

)

と考えて

Romberg

加速を用いることができる.

● 講義資料

▼ 講義予定

•

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

–

数値解法とは何か

–

（前進）オイラー法と後退オイラー法

● 講義資料

「常微分方程式の数値計算」の授業で題材にする代表的な常微分方程式は以下の通り.

x^′(t) =λx(t), x(0) =x0, λ6= 0. (2)

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

x^′′(t) +ax^′(t) +bx(t) =f(x), x(0) =x0, x^′(0) =v0. (4) x^′′(t) + sin(x(t)) = 0, x(0) =x0, x^′(0) =v0. (5) x^′′(t) + (x(t)²−1)x^′(t) +x(t) = 0, x(0) =x0, x^′(0) =v0. (6)

• (3)

は

“logistic equation”

と呼ばれ, 初期条件が

−1< x0 <1

をみたすとき, 解は

t∈(−∞,∞)

上 で

−1< x(t)<1

をみたす.

• (4)

は

a= 0,b=k²,f(x)≡0

の時, “単振動” の方程式である.

• (5)

は

“単振り子”

の方程式である.

• (6)

は

“Val del Pol’s equation”

と呼ばれ, 非線形抵抗を持つ電気回路の方程式である.

これらの中で

(2), (3), (4)

は, 厳密な解（解析解）を容易に求めることができるので, 数値解がどのくらい 解析解を近似しているかなどの解析に用いることができる.

★ 計算結果

•

（前進）オイラー法による計算

1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8

0 0.2 0.4 0.6 0.8 1 1.2

x’=x, x(0) = 1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 0.5 1 1.5 2 2.5

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

(2)λ= 1,x0= 1, h= 0.01 (3)h= 0.01

-1.5 -1 -0.5 0 0.5 1 1.5

0 2 4 6 8 10 12

x’’+x=0, x(0) = 1, x’(0) = 0

-1.5 -1 -0.5 0 0.5 1 1.5

-1.5 -1 -0.5 0 0.5 1 1.5

x’’+x=0, x(0) = 1, x’(0) = 0

(4)a= 0,b= 1,x0= 1, v0= 0

-1.5 -1 -0.5 0 0.5 1 1.5

0 5 10 15 20 25

x’’+sin(x)=0, x(0) = 1, x’(0) = 0

-1 -0.5 0 0.5 1

-1 -0.5 0 0.5 1

x’’+sin(x)=0, x(0) = 1, x’(0) = 0

(5)x0= 1,v0= 0

-2 -1 0 1 2

0 10 20 30 40 50

x’’+(x^2-1)x’+x=0, x(0) = 1, x’(0) = 1

-3 -2 -1 0 1 2 3

-3 -2 -1 0 1 2 3

x’’+(x^2-1)x’+x=0, x(0) = 1, x’(0) = 1

(6)x0= 1,v0= 0

•

後退オイラー法による計算

1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8

0 0.2 0.4 0.6 0.8 1 1.2

x’=x, x(0) = 1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 0.5 1 1.5 2 2.5

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

(2)λ= 1,x0= 1, h= 0.01 (3)h= 0.01

-1.5 -1 -0.5 0 0.5 1 1.5

0 2 4 6 8 10 12

x’’+x=0, x(0) = 1, x’(0) = 0

-1.5 -1 -0.5 0 0.5 1 1.5

-1.5 -1 -0.5 0 0.5 1 1.5

x’’+x=0, x(0) = 1, x’(0) = 0

(4)a= 0,b= 1,x0= 1, v0= 0

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

0 5 10 15 20 25

x’’+sin(x)=0, x(0) = 1, x’(0) = 0

-1 -0.5 0 0.5 1

-1 -0.5 0 0.5 1

x’’+sin(x)=0, x(0) = 1, x’(0) = 0

(5)x0= 1,v0= 0

-2 -1 0 1 2

0 10 20 30 40 50

x’’+(x^2-1)x’+x=0, x(0) = 1, x’(0) = 1

-3 -2 -1 0 1 2 3

-3 -2 -1 0 1 2 3

x’’+(x^2-1)x’+x=0, x(0) = 1, x’(0) = 1

(6)x0= 1,v0= 0

•

前進オイラー法と後退オイラー法での

(2)

の収束の様子

λ= 1,x0= 1,h= 0.01

1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3

0 0.2 0.4 0.6 0.8 1 1.2 1

1.5 2 2.5 3 3.5 4

0 0.2 0.4 0.6 0.8 1 1.2

前進オイラー法 後退オイラー法

• (2) (λ= 1)

の

h

を動かしたときの誤差.

t= 1.0

における

e

との誤差. 横軸は

1/h.

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

10 100 1000 10000 100000 1e+06 1e+07 1e+08 1e+09 1e+10 1e+11 x’=x, x(0) = 1 error

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

10 100 1000 10000 100000 1e+06 1e+07 1e+08 1e+09 1e+10 1e+11 x’=x, x(0) = 1 error

前進オイラー法 後退オイラー法

● 実習内容

1. (2), (3), (4)

の真の解（「解析解」）を求めなさい.

2.

時間刻み幅, 解を求める時間の範囲を適当に決めて, (2), (3), (4), (5), (6) の数値解を, （前進）オイ ラー法を用いて構成しなさい. また, 解析解を計算できる方程式に関しては, 時間刻み幅を変化させ たとき, 解析解との誤差がどうなるかを観察しなさい.

3.

時間刻み幅, 解を求める時間の範囲を適当に決めて, (2), (3), (4) の数値解を, 後退オイラー法を用い て構成しなさい. また, 解析解を計算できる方程式に関しては, 時間刻み幅を変化させたとき, 解析解 との誤差がどうなるかを観察しなさい.

4. (4)

で

a= 0,b=k²,f(x)≡0

とすると, (x

^′(t))²+ (x(t))²

は時間に寄らず一定であることを示しな さい. 前進オイラー法・後退オイラー法で

(4)

を解いたとき, (x

^′(t))²+ (x(t))²

の値がどのように変 化するかを観察しなさい.

5. (5)

では,

¹₂(x^′(t))²−cos(x(t))

が時間に寄らず一定の値をとることを示しなさい. また, 前進オイラー 法で

(5)

を解いたとき,

¹₂(x^′(t))²−cos(x(t))

の値がどのように変化するかを観察しなさい.

6.

時間刻み幅, 解を求める時間の範囲を適当に決めて, (5) の数値解を, 後退オイラー法を用いて構成し なさい. また, このとき

¹₂(x^′(t))²−cos(x(t))

の値がどのように変化するかを観察しなさい.

7.

時間刻み幅, 解を求める時間の範囲を適当に決めて, (5) の数値解を, 後退オイラー法を用いて構成し なさい.

8.

後退オイラー法の局所離散化誤差, 大域的離散化誤差, 丸め誤差を評価しなさい.

★ サンプル

【オイラー法のサンプルプログラム】

x^′=f(x)

をオイラー法で解くプログラムの一つの例.

/* 前進 Euler 法 */

#include <stdio.h>

#define X0 (1.0) /* 初期値 */

#define H (0.01) /* 時間刻み幅 */

#define T (1.0) /* 終了時間 */

int main(int argc, char **argv) {

double t=0.0 ;

double x, x0 ;

x0 = X0 ; while(t<(T+H)) {

printf("%.16e\t%.16e\n", t, x0) ; x = x0 + H*f(x0) ;

x0 = x ; t += H ; }

return 0 ; }