これをAitken 加速と呼ぶ

● 前回の講義のまとめ

【加速】

1. 数列{a_n} が極限値αを持ち,

|a_n−α|=c1λⁿ₁ +c2λ²₂+· · · , 1> λ1> λ2>· · ·>0 という振る舞いをしているとき,

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

1−λ1

で定義された数列{a⁽¹⁾n }も同じ極限に収束し,

|a⁽¹⁾_n −α|=c₂λ²₂+· · · をみたす. これをRomberg加速と呼ぶ.

2. 数列{a_n} が極限値αを持ち,

(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 加速と呼ぶ.

3. 半径１の円に内接する正n角形の周長を 2n とおいた,列{_n} は, k−π=π³

6 k⁻²− π⁵

120k⁻⁴+O(k⁻⁶) をみたすので,an=n02ⁿ とおくと,

an−π= π³

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

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

をみたす. したがって,λ1= 1/4 としてRomberg-Richardson加速によって収束を加速するこ とができる.

4. Romberg加速は「収束の主要項」λ1の値が分からなければ加速することができないが, Aitken

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

【巾級数による関数の計算】

1. e^x, sin(x), arctan(x), log(x)などの初等超越関数を数値的に計算する（与えられたxに対する 関数値を計算する）ためには,テイラー展開を用いる.

2. arctan(x)を計算するには,そのテイラー展開 arctan(x) =x−x³

3 +x⁵

5 +· · ·+(−1)ⁿx²ⁿ⁺¹ 2n+ 1 +· · · に対して,

αn(x) =x−x³ 3 +x⁵

5 +· · ·+(−1)ⁿx²ⁿ⁺¹ 2n+ 1

を計算すればよい. しかし, arctan(x)のx= 0におけるテイラー展開の収束域は|x|<1であ ることに注意が必要である. （講義で書いたテイラー級数は間違っていました. ごめんなさい.）

3. 実際,x= 1に対して,

αn(1) = 1−1 3 +1

5 +· · ·+ (−1)ⁿ 2n+ 1 は収束することが知られているが,その打ち切り誤差Rn(x)は

|R_n(x)|=

_x

0

t²ⁿ⁺¹ 1 +t²dt

≤

_x

0 t²ⁿ⁺¹dt= x²⁽ⁿ⁺¹⁾ 2(n+ 1) =

O(x²ⁿ) |x|<1 O(n⁻¹) x= 1

と評価できる. すなわち,収束半径内と収束円上では収束の速度が全く異なることがわかる.

● 講義資料

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

x(t) =λx(t), x(0) =x0, λ= 0. (1)

x(t) = (1−x(t)²), x(0) = 0. (2)

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

• (2)は“logistic equation”と呼ばれ,初期条件が−1< x0<1をみたすとき,解はt∈(−∞,∞) 上で−1< x(t)<1 をみたす.

• (3)は a= 0,b=k²,f(x)≡0の時, “単振動”の方程式である.

• (4)は “単振り子”の方程式である.

• (5)は “Val del Pol’s equation”と呼ばれ,非線形抵抗を持つ電気回路の方程式である.

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

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

(1)λ= 1,x0= 1, h= 0.01 (2)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

(3)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.5 -1 -0.5 0 0.5 1 1.5

-1.5 -1 -0.5 0 0.5 1 1.5

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

(4)x0= 1,v0= 0

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

0 10 20 30 40 50 60

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

-3 -2 -1 0 1 2 3

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

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

(5)x0= 1,v0= 0 3. 後退オイラー法による計算

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

(1)λ= 1,x0= 1, h= 0.01 (2)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

-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

(4)x0= 1,v0= 0

4. 前進オイラー法と後退オイラー法での (1)の収束の様子λ= 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

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

5. (1) (λ= 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. (1), (2), (3)の真の解（「解析解」）を求めなさい.

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

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

4. (3) でa= 0,b=k², f(x)≡0とすると, (x(t))²+ (x(t))は時間に寄らず一定であることを示しな さい. 前進オイラー法・後退オイラー法で(3)を解いたとき, (x(t))²+ (x(t))の値がどのように変化 するかを観察しなさい.

5. (4)では, (x(t))²−cos(x(t))が時間に寄らず一定の値をとる. このことを示しなさい. また,前進オ イラー法で(4)を解いたとき, ¹₂(x(t))²−cos(x(t))の値がどのように変化するかを観察しなさい.

6. 時間刻み幅,解を求める時間の範囲を適当に決めて, (4)の数値解を,後退オイラー法を用いて構成し なさい. また,このとき ¹₂(x(t))²−cos(x(t))の値がどのように変化するかを観察しなさい.

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

8. 電子メールで「今日の講義の感想や意見」を送ってください.