● 講義資料

● 講義資料

▼ 講義予定

• 収束の加速

• 常微分方程式の初期値問題の数値解法 – 数値解とは何か

– 前進オイラー法

● 前回の講義のまとめ

▼ 浮動小数点演算の減算に関する相殺

• 浮動小数点数x,y∈F の減算に関しては,保護桁が1桁あれば,相対誤差2ǫ以内で結果を 求めることができる.

• しかし,a= 1.21×10⁻²,b= 1.00,c= 1.23×10⁻²とおき,二次方程式ax²+bx+c= 0の 解を考えると,真の解は,

−b+√

b²−4ac

2a =−0.0123018,

−b−√

b²−4ac

2a =−82.6323 である. これを, 10進3 桁浮動小数点数の体系で計算すると,

−b+√

b²−4ac

2a =−4.07×10⁻²,

−b−√

b²−4ac

2a =−8.14×10¹

となり, それぞれの相対誤差は, 2.30846, 0.014913となる. すなわち, ^−b+√₂^ba^2−4ac を正しく 計算していないことがわかる.

この現象は,以下の相殺(canceling)と呼ばれるものである.

• 一般にx∈F を考えると,その有効桁の下の方の桁は「正しい値」を表しているとはいえな い. x∈Rの近似値としてF の値を考えているので,最低でも「最後の桁」は正しい値では ないと考えるのが妥当である.この「正しい値」を表していると考えられない桁のことを「汚 染桁」と呼ぶ.

• x,y∈F が非常に近い値を持つとき,その減算を行うと,汚染桁と考えられる部分が,x⊖yの上 位桁にあらわれる. たとえば, 10進3桁浮動小数点数の体系でx= 1.02×10⁰,y= 1.01×10⁰ としたとき,それぞれの最後の桁は汚染桁と考えて良い. この減算結果はx⊖y= 1.00×10⁻² となり, 汚染桁から計算された「1」が最上位桁となり, このx⊖y の値は「信用できない」

値となっている.

• 一般に,浮動小数点数の計算において「値が近い数の減算」を行うと,その結果は大きな誤差 を生じることがあることがわかる.

• 最初の二次方程式の解の公式については,b >0 の時には,

−b+√

b²−4ac

2a = 2c

−b−√

b²−4ac と書き直すことにより,

2c

−b−√

b²−4ac =−1.21×10⁻² と計算され,相対誤差は0.0164041となり,許容範囲内となる.

より一般には,二次方程式ax²+bx+c= 0の解を公式にしたがって求める際には,

−b−√

b²−4ac

2a , 2c

−b−√

b²−4ac, b >0,

−b+√

b²−4ac

2a , 2c

−b+√

b²−4ac, b <0 と計算すべきであることがわかる.

▼ 正多角形の周長による円周率の近似計算

• 円周率πの定義として,「単位円周の長さの２倍」または「単位円の面積」を採用して,πの 近似値を求めることを考える. （これは,厳密な定義ではないが,紀元前からπを計算する前 提とされていた.）

• ℓⁿ を単位円に内接する正 n 角形の周長の 1/2倍,Lⁿ を単位円に外接する正 n角形の周長 の1/2倍とおくと,容易にℓn< π < Ln であることがわかる. （ℓn< πの証明には,「２点 を結ぶ最短線は直線である」との事実を使い,π < Lⁿ の証明には,面積の比較を用いる.）

• 簡単な図形的考察と三角関数の性質から ℓⁿ=nsin(π

n), Lⁿ=nsin(π n), ℓⁿ< ℓ2n<· · ·< π <· · ·< L2n< Lⁿ,

nlim→∞ℓⁿ= lim

n→∞Lⁿ=π であることがわかる.

• また,テイラーの定理を用いれば, ℓⁿ =nsin(π

n) =π− π³

6n² + π⁵

120n⁴ +O(n⁻⁶), Lⁿ =ntan(π

n) =π+ π³

3n² + 2π⁵

15n⁴+O(n⁻⁶), であることもわかる.

• 半角の公式を用いると,

ℓ2n = 2n s

1−p

1−(ℓⁿ/n)²

2 ,

L2n=2n² Lⁿ

p1 + (Lⁿ/n)²−1

が成り立つ. また,

2 L2n

= 1 ℓⁿ + 1

Lⁿ が成り立つ.

• この公式をn= 6,ℓ6= 3,L6= 2√

3を初期値として用いて計算すると, 10⁻⁸ 程度の精度し かπの近似値を得ることができない.

これは, 上の式でnが大きくなったときに「近い値同士の減算」を行っていて, それにとも

なうcancelingが原因であることがわかる.

• これを解消する方法として,

ℓ2n=ℓⁿ

s 2 1 +p

1−(ℓ/n)² L2n= 2Ln

p1 + (Lⁿ/n)²+ 1

と書き換えることにより, 10⁻¹⁴程度の精度でπを求めることが可能となる.

• また,より計算が単純な方法として, 1 L2n

=1 2

1 ℓⁿ + 1

Lⁿ

, ℓ2n =p

ℓⁿL2n, を得ることができる.

▼ 巾級数による

e , π

• 自然対数の底e, 円周率πの近似値を求めるために広く使われる方法は,初等超越関数の特 殊値として求める方法である.

• e^x= 1 +x+^x₂² +· · ·+^x_n!ⁿ+· · · であり,この巾級数の収束半径は無限大であるので,x= 1 を代入することにより,

e= 1 + 1 +1

2+· · ·+ 1 n!+· · ·

の右辺の無限級数の和を近似することでeの近似値を得ることができる.

• いま,

e^x=

n

X

k=0

x^k

k! + ξⁿ⁺¹

(n+ 1)!, 0< ξ <1 であることを使えば,

e−

n

X

k=0

1 k!

≤ 1

(n+ 1)!

が成り立つ. そこで, 1

(n+ 1)! < ǫが成り立つnまで計算すれば, 1

(n+ 1)! < ǫ < ǫeとなり, eを相対誤差ǫで求めることができる.

• このように,ほしい精度で計算するために,無限の計算を適切な有限の計算で打ち切ったこと による誤差を打ち切り誤差とよび,上の計算は「打ち切り誤差((n+ 1)!)⁻¹ で計算する」と いう言い方をする.

• 円周率πは

π

4 = arctan(1)

としてarctan(x)の特殊値として得ることができる.

• したがって, arctan(x)のテイラー展開を考え,そこにx= 1を代入すればよいように思える が,これは誤りである. （誤りの理由は２つある）

• arctan(x)のテイラー展開は

arctan(x) =x−x³ 3 +x⁵

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

であるが,右辺の巾級数の収束半径は1であるので,両辺に1を代入することはできない. 右 辺の巾級数がx= 1 で収束すれば,その値がarctan(1)に一致することはアーベルの定理の 主張である. すなわち, 右辺の巾級数が x= 1で収束するかどうかは自明ではない. （これ が「誤りの第一の理由」である.）

• いま,初項1,公比−x² の有限等比級数を考えると 1−x²+x⁴+· · ·+ (−1)ⁿx²ⁿ= 1−(−1)ⁿx²ⁿ⁺²

1 +x² = 1

1 +x² −(−1)ⁿx²ⁿ⁺² 1 +x²

となる. この両辺を積分することにより（左辺は有限和であるので,項別に積分ができる）

x−x³ 3 +x⁵

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

2n+ 1 = arctan(x)− Z ^x

0

(−1)ⁿt²ⁿ⁺² 1 +t² dt を得る. すなわち,

x−x³ 3 +x⁵

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

2n+ 1 −arctan(x) ≤

Z ^x

0

t²ⁿ⁺² 1 +t²dt

≤ Z ^x

0

t²ⁿ⁺²dt≤









 1

2n+ 3 x= 1, x²ⁿ⁺³

2n+ 3 |x|<1, となる.

• すなわち, arctan(x) のテイラー展開による計算は, その打ち切り誤差が x = 1 の時には O(n⁻¹),|x|<1の時にはO(x⁻ⁿ)となり,収束の様子が本質的に異る. 別の言い方をすれば,

arctan(1) =π/4の近似値を求めるために, arctan(1)のテイラー展開（これが収束すること

は上記の計算でわかっている）を計算することは望ましくない. （これが「誤りの第二の理 由」である.）

• arctan(1) =π/4の値を計算する高速な方法としては,たとえば,以下のものが知られている.

π

4 = 4 arctan(1/5)−arctan(1/239), (マチンの公式) π

4 = arctan(1/2) + arctan(1/3), , (ガウス) π

4 = 12 arctan(1/49) + 32 arctan(1/57)−5 arctan(1/239) + 12 arctan(1/110443),

（他にも山ほどあるし,他の方法でπを高速に計算することは多くの方法が知られている.）

● 講義資料

▼ 収束の加速

次の図は, 単位円に接する正多角形の周長による円周率の近似計算を, Ronberg加速と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

次の図は, 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) 2−cos(x1−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)

-3 -2 -1 0 1 2 3 4 5 6 7

0 5 10 15 20

Double Pendulum

-3 -2 -1 0 1 2 3 4 5 6 7

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

二重振り子,x1(0) = 2,x2(0) = 0,x^′₁(0) =x^′₂(0) = 0

以下の図は,

x^′(t) = 4

sign(x(t))p

|x(t)|+ max

0, t−|x(t)| t

cos

πlog(t) log(2)

, x(0) = 0

を hを取り替えて計算した. h= 2^−k ととり計算すると,kの偶奇により,解の収束先が異る.

-4 -3 -2 -1 0 1 2 3 4

0 0.2 0.4 0.6 0.8 1

実際には,t= 0,y= 0からは計算できないので,t=h, y= 0から計算を開始している.

▼ サンプルプログラム

#include <stdio.h>

#define T (1.0) /* ここで指定した T まで計算する */

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

#define X0 (1.0) /* 初期条件 */

double f(double t, double x) ;

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

double t, x, new_x ;

t = 0.0 ; x = X0 ;

printf("%f %.15e\n", t, x) ; while(t < T) {

new_x = x + H*f(t, x) ; t += H ;

printf("%f %.15e\n", t, new_x) ; x = new_x ;

}

return 0 ; }

/* 微分方程式を定義する */

double f(double t, double x) {

return x ; }

● 実習内容

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

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