● 前回の講義のまとめ

● 前回の講義のまとめ

• 数列{an}が極限値 αを持ち,

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

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

1−λ1

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

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

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

(a_n+1−α) = (A+ε_n)(a_n−α), ε_n→0, |A|<1 という振る舞いをしているとき,

a⁽¹⁾_n =a_n− (a_n−a_n+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角形の周長を 2nとおいた,列 {n} は, k−π= π³

6 k⁻²− π⁵

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

an−π= π³

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

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

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

• ニュートン法 (m 重根)の場合には, 各ステップでの誤差 ε_n =|a_n−α| はε_n+1 = ^m−1_m ε_n となる ので,

a_n =α+C

m−1

m _n

+o(

m−1

m _n

)

と考えて Romberg加速を用いることができる.

● 講義資料

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< x₀<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 -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

-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

(5)x₀= 1,v₀= 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

(3)a= 0,b= 1,x₀= 1, v₀= 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

-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

(5)x₀= 1,v₀= 0

4. 前進オイラー法と後退オイラー法での (1)の収束の様子λ= 1,x₀= 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. 時間刻み幅,解を求める時間の範囲を適当に決めて, (4)の数値解を,後退オイラー法を用いて構成し なさい.

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 ; }

【２階常微分方程式の場合の出力と gnuplotの利用法】

２階常微分方程式で「相空間」の図を出力するには,x(t)に対応する値だけではなく,x(t)に対応す る値も出力する必要がある. この場合,第１列にt, 第２列にx(t),第３列にx(t)の値を出力したと すると, gnuplotにおいて,

set size ratio 1

plot ‘‘euler.out’’ using 2:3 with p と入力すればよい.

• １行目は描画する図の縦横の比率を1:1にする指定である.

• ２行目の“using 2:3”を指定することにより,データファイルのの第２列目と第３列目を入力

データとすることができる.