● 講義資料
▼ 講義予定
•
ニュートン法・逐次近似とその収束の判定● 前回の講義のまとめ
▼ 浮動小数点演算の減算に関する相殺
•
浮動小数点数x, y ∈ F
の減算に関しては,保護桁が1
桁あれば,相対誤差2ǫ
以内で結果を 求めることができる.•
しかし,a = 1.21 × 10 − 2 , b = 1.00, c = 1.23 × 10 − 2
とおき,二次方程式ax 2 + bx + c = 0
の 解を考えると,真の解は,− b + √
b 2 − 4ac
2a = − 0.0123018,
− b − √
b 2 − 4ac
2a = − 82.6323
である. これを, 10進3
桁浮動小数点数の体系で計算すると,− b + √
b 2 − 4ac
2a = − 4.07 × 10 − 2 ,
− b − √
b 2 − 4ac
2a = − 8.14 × 10 1
となり, それぞれの相対誤差は, 2.30846, 0.014913となる. すなわち,
− b + √ 2 b a
2− 4 ac
を正しく 計算していないことがわかる.この現象は,以下の相殺
(canceling)
と呼ばれるものである.•
一般にx ∈ F
を考えると,その有効桁の下の方の桁は「正しい値」を表しているとはいえな い.x ∈ R
の近似値としてF
の値を考えているので,最低でも「最後の桁」は正しい値では ないと考えるのが妥当である.この「正しい値」を表していると考えられない桁のことを「汚 染桁」と呼ぶ.• x, y ∈ F
が非常に近い値を持つとき,その減算を行うと,汚染桁と考えられる部分が,x ⊖ y
の上 位桁にあらわれる. たとえば, 10進3
桁浮動小数点数の体系でx = 1.02 × 10 0 , y = 1.01 × 10 0
としたとき,それぞれの最後の桁は汚染桁と考えて良い. この減算結果はx ⊖ y = 1.00 × 10 − 2
となり, 汚染桁から計算された「1」が最上位桁となり, このx ⊖ y
の値は「信用できない」値となっている.
•
一般に,浮動小数点数の計算において「値が近い数の減算」を行うと,その結果は大きな誤差 を生じることがあることがわかる.•
最初の二次方程式の解の公式については,b > 0
の時には,− b + √
b 2 − 4ac
2a = 2c
− b − √
b 2 − 4ac
と書き直すことにより,
2c
− b − √
b 2 − 4ac = − 1.21 × 10 − 2
と計算され,相対誤差は0.0164041
となり,許容範囲内となる.より一般には,二次方程式
ax 2 + bx + c = 0
の解を公式にしたがって求める際には,− b − √
b 2 − 4ac
2a , 2c
− b − √
b 2 − 4ac , b > 0,
− b + √
b 2 − 4ac
2a , 2c
− b + √
b 2 − 4ac , b < 0
と計算すべきであることがわかる.▼ 正多角形の周長による円周率の近似計算
•
円周率π
の定義として,「単位円周の長さの２倍」または「単位円の面積」を採用して,π
の 近似値を求めることを考える. （これは,厳密な定義ではないが,紀元前からπ
を計算する前 提とされていた.）• ℓ n
を単位円に内接する正n
角形の周長の1/2
倍,L n
を単位円に外接する正n
角形の周長 の1/2
倍とおくと,容易にℓ n < π < L n
であることがわかる. （ℓn < π
の証明には,「２点 を結ぶ最短線は直線である」との事実を使い,π < L n
の証明には,面積の比較を用いる.）•
簡単な図形的考察と三角関数の性質からℓ n = n sin( π
n ), L n = n tan( π n ), ℓ n < ℓ 2n < · · · < π < · · · < L 2n < L n ,
n lim →∞ ℓ n = lim
n →∞ L n = π
であることがわかる.•
また,テイラーの定理を用いれば,ℓ n = n sin( π
n ) = π − π 3
6n 2 + π 5
120n 4 + O(n − 6 ), L n = n tan( π
n ) = π + π 3
3n 2 + 2π 5
15n 4 + O(n − 6 ),
であることもわかる.•
半角の公式を用いると,ℓ 2 n = 2n s
1 − p
1 − (ℓ n /n) 2
2 ,
L 2 n = 2n 2 L n
p 1 + (L n /n) 2 − 1
が成り立つ. また,2 = 1
+ 1
•
この公式をn = 6, ℓ 6 = 3, L 6 = 2 √
3
を初期値として用いて計算すると, 10− 8
程度の精度し かπ
の近似値を得ることができない.これは, 上の式で
n
が大きくなったときに「近い値同士の減算」を行っていて, それにともなう
canceling
が原因であることがわかる.•
これを解消する方法として,ℓ 2n = ℓ n
s 2 1 + p
1 − (ℓ/n) 2
L 2n = 2L n
p 1 + (L n /n) 2 + 1
と書き換えることにより, 10
− 14
程度の精度でπ
を求めることが可能となる.•
また,より計算が単純な方法として,1 L 2 n
= 1 2
1 ℓ n + 1
L n
, ℓ 2 n = p
ℓ n L 2 n ,
を得ることができる.▼ 巾級数による
e
, π
の近似値の計算•
自然対数の底e,
円周率π
の近似値を求めるために広く使われる方法は,初等超越関数の特 殊値として求める方法である.• e x = 1 + x + x 2
2+ · · · + x n!
n+ · · ·
であり,この巾級数の収束半径は無限大であるので,x = 1
を代入することにより,e = 1 + 1 + 1
2 + · · · + 1 n! + · · ·
の右辺の無限級数の和を近似することで
e
の近似値を得ることができる.•
いま,e x =
n
X
k =0
x k
k! + ξ n +1
(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)!) − 1
で計算する」と いう言い方をする.•
円周率π
はπ
4 = arctan(1)
として
arctan(x)
の特殊値として得ることができる.•
したがって, arctan(x)のテイラー展開を考え,そこにx = 1
を代入すればよいように思える が,これは誤りである. （誤りの理由は２つある）• arctan(x)
のテイラー展開はarctan(x) = x − x 3 3 + x 5
5 + · · · + ( − 1) n x 2 n +1 2n + 1 + · · ·
であるが,右辺の巾級数の収束半径は
1
であるので,両辺に1
を代入することはできない. 右 辺の巾級数がx = 1
で収束すれば,その値がarctan(1)
に一致することはアーベルの定理の 主張である. すなわち, 右辺の巾級数がx = 1
で収束するかどうかは自明ではない. （これ が「誤りの第一の理由」である.）•
いま,初項1,
公比− x 2
の有限等比級数を考えると1 − x 2 + x 4 + · · · + ( − 1) n x 2 n = 1 − ( − 1) n x 2n+2
1 + x 2 = 1
1 + x 2 − ( − 1) n x 2n+2 1 + x 2
となる. この両辺を積分することにより（左辺は有限和であるので,項別に積分ができる）
x − x 3 3 + x 5
5 + · · · + ( − 1) n x 2 n +1
2n + 1 = arctan(x) − Z x
0
( − 1) n t 2 n +2 1 + t 2 dt
を得る. すなわち,x − x 3 3 + x 5
5 + · · · + ( − 1) n x 2n+1
2n + 1 − arctan(x) ≤
Z x
0
t 2n+2 1 + t 2 dt
≤ Z x
0
t 2 n +2 dt ≤
1
2n + 3 x = 1, x 2n+3
2n + 3 | x | < 1,
となる.•
すなわち, arctan(x) のテイラー展開による計算は, その打ち切り誤差がx = 1
の時にはO(n − 1 ), | x | < 1
の時にはO(x − n )
となり,収束の様子が本質的に異る. 別の言い方をすれば,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),
● 講義資料
▼ ニュートン法
•
以下の図は,ニュートン法をf (x) = x 2 − 2
に適用し,√
2
の近似値を求めた例である. （√ 2
と近似値の値との絶対誤差を表示している）1e-16 1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1
0 1 2 3 4 5
Newton method for sqrt(2)
•
以下の図は, ニュートン法をf k (x) = (x 2 − 2) k
に適用し,√
2
の近似値を求めた例である.（
√
2
と近似値の値との絶対誤差を表示している）1e-16 1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1
0 20 40 60 80 100 120
Newton method for sqrt(2)
k=1 k=2 k=3 k=4
•
以下の図は,各種の求根アルゴリズムをf k (x) = x 2 − 2
に適用して√
2
の近似値を求めた例 である. （√
2
と近似値の値との絶対誤差を表示している）1e-16 1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1
0 2 4 6 8 10 12
Some methods for sqrt(2)
newton secant bisection (1) bisection (2)
ここで,用いた方法は以下の通りである.
– “newton”
は「ニュートン法」を,a 0 = 2.0
で適用.– “secant”
は「割線法」を,a 0 = 2.0, a 1 = 2.0 − 0.01
で適用.– “bisection(1)”
は「二分法」を,I 0 = [0.0, 2.0]
で適用.– “bisection(2)”
は「改良した二分法」を,I 0 = [0.0, 2.0]
で適用.「改良した二分法」とは以下のものである. 通常の二分法では,
I n = [a n , b n ], f (a n ) < 0, f (b n ) > 0
に対して,c n +1 = (a n + b n )/2
とおき,f (c n +1 )
の符号を判定している. その 代りに, (an , f (a n )), (b n , f (b n ))
を結ぶ直線とy = 0
との交点をc n +1
とおき,f (c n +1 )
の符号を判定したものである.● 実習内容
（以下で「★」の数は推奨の程度を示します. 多いほど推奨の度合いが大きい. また,「†」は 難易度またはマニアックな程度を表します. 多いほど難しいかマニアックかです.）
