● 前回の講義のまとめ
•
（代数）方程式の解を求めるための方法として,以下を考察した.1.
区間縮小法（二分法）2.
ニュートン法（逐次近似）•
区間縮小法は「中間値の定理」を応用したものであり,連続な関数f : R −→ R
に対して,f ( a ) < 0, f ( b ) > 0, f ( x ) > 0 on [ a, b ]
をみたすとき,f ( x ) = 0
の区間( a, b )
での唯一の解α
を求める方法で ある.•
ニュートン法は,以下の定理に基づく方法である.なめらかな関数
f : R −→ R
に対して,f ( a ) < 0, f ( b ) > 0, f ( x ) > 0 on [ a, b ], f ( x ) < 0 on ( a, b )
をみたすとき,x n+1 = x n − f ( x n )
f ( x n ) , x 0 = b
によって定められる数列
{x n }
は,f ( x ) = 0
の区間( a, b )
での唯一の解α
に収束する.•
より一般に,なめらかな関数ϕ : R −→ R
が,次の条件を満たすと仮定する.ある正の定数
L < 1
が存在して,任意のx , y ∈ ( a, b )
に対して|ϕ ( x ) − ϕ ( y ) | ≤ L|x − y| (1)
を満たす. このとき,
ϕ
は区間( a, b )
上にただ一つの固定点α ( ϕ ( α ) = α )
を持つ. 特に,x n+1 = ϕ ( x n )
によって定められる数列
{x n }
はα
に収束する.–
特に,|ϕ ( α ) | < 1
が成り立つとき,α
を含む区間( a, b )
を十分小さく取れば, (1) が成り立つ.–
このとき,|x n+1 − α| ≤ |ϕ ( α ) | |x n − α| + |ϕ ( α ) |
2 |x n − α| 2 + O ( |x n − α| 3 )
が成り立つ.–
ニュートン法は,ϕ ( x ) = x − f f(x)(x)
と置いた逐次近似と考えられ,
ϕ ( x ) = − f ( x ) f ( x )
( f ( x )) 2 , ϕ ( x ) = f ( x )
f ( x ) − 2 f ( x )( f ( x )) 2 ( f ( x )) 3 が成り立つ.
–
特に,解が単根( f ( α ) = 0)
のとき,ϕ ( α ) = 0 , ϕ ( α ) = f ( α ) f ( α )
となり,|x n+1 − α| ≤ M |x n − α| 2 + O ( |x n − α| 3 )
–
また,解がm
重根（f ( α ) = · · · = f (m) ( α ) = 0
のとき）|x n+1 − α| ≤ m − 1
m |x n − α|
となり, 一次収束する.
•
ニュートン法では,「解α
の近似値を,相対誤差δ
で求める」ために,いつ繰り返しを停止させるかは明 らかではない. すなわち,真の値α
と「近似値」x nを直接比較できないため,評価式“ |x n −α| ≤ δ|α| ”
を利用することができない. そのため,上の評価式（二次収束の式,または,一次収束の式）を利用す
ることにより, 次のように評価することができる. （以下では, ε n = |x n − α| , e n = |x n+1 − x n |
と
おく.）
★ 単根の時
|x n+1 − x n | ≤ 2 δ
が成り立てば,|x n − α| ≤ δ|x n |
が成り立つ.【証明】 三角不等式より,
ε n = |x n − α| ≤ |x n+1 − x n | + |x n+1 − α| = e n + ε n+1 が成り立つ. 一方,二次収束の式から,
ε n+1 ≤ Mε 2 n が成り立つので,
(1 − Mε n ) ε n ≤ e n
が成り立つ. ここで,
n
を十分大きくとれば, 1− Mε n > 1 / 2
と取ることができる. した がってε n ≤ 2 e n , |x n − α| ≤ 2 |x n+1 − x n |
が成り立つ.いま,
|x n+1 − x n | ≤ (1 / 2) δ|x n |
と取れば,|x n − α| ≤ 2 |x n+1 − x n | ≤ δ|x n |
となり,このとき
x n は相対誤差δ
でα
を近似していることがわかる.
★
m
重根の時|x n+1 − x n | ≤ mδ
が成り立てば,|x n − α| ≤ δ|x n |
が成り立つ.【証明】 この場合は一次収束の式
ε n+1 ≤ m − 1 m ε n
が成り立つので,三角不等式より
ε n ≤ me n , |x n − α| ≤ m|x n+1 − x n |
が成り立つ.いま,
|x n+1 − x n | ≤ (1 /m ) δ|x n |
と取れば,|x n − α| ≤ m|x n+1 − x n | ≤ δ|x n |
となり,このとき
x n は相対誤差δ
でα
を近似していることがわかる.
★ プログラムサンプル
【区間縮小法で
√
2
の近似値を計算する】√ 2
の近似値を絶対誤差10 −10と相対誤差10 −10で求める. 変数r
は区間の右端点を,l
は左端点を
表す. 常にr > l > 0
が成り立つので,相対誤差の評価はl − r < δr
を用いればよい
r
は区間の右端点を,l
は左端点を 表す. 常にr > l > 0
が成り立つので,相対誤差の評価はl − r < δr
を用いればよい/*
区間縮小法でSQRT_2
を計算する*
*
初期条件でf(l) < 0, f(r) > 0
を満たす*/
#include <stdio.h>
#include <math.h>
#define D 1.0E-10 double f(double) ;
int main(int argc, char **argv) {
double l, r, c ;
/*
絶対誤差D
で求める*/
l = 0.0 ; r = 2.0 ; while(r - l > D) {
c = (r+l)/2.0 ;
/*
この出力行はDEBUG
用*/
/* printf("%.16f\t%+.16e\n", c, c - M_SQRT2) ; */
if (f(c) < 0) l = c ; else r = c ;
}
printf("SQRT2 = %.16f\n", c) ; /*
相対誤差D
で求める*/
l = 0.0 ; r = 2.0 ; while(r - l > D*r) {
c = (r+l)/2.0 ;
/*
この出力行はDEBUG
用*/
/* printf("%.16f\t%+.16e\n", c, (c - M_SQRT2)/M_SQRT2) ; */
if (f(c) < 0) l = c ; else r = c ;
}
printf("SQRT2 = %.16f\n", c) ; return 0 ;
}
double f(double x) {
return x*x - 2.0 ;
}
● 講義資料
【多角形近似による
π
の計算のRomberg
加速】1 3.000000000000
2 3.105828541230 3.141561970632
3 3.132628613281 3.141590732969 3.141592650458 4 3.139350203047 3.141592533505 3.141592653541 5 3.141031950891 3.141592646084 3.141592653589 6 3.141452472285 3.141592653121 3.141592653590 7 3.141557607912 3.141592653560 3.141592653590 8 3.141583892148 3.141592653588 3.141592653590 9 3.141590463228 3.141592653590 3.141592653590 10 3.141592105999 3.141592653590 3.141592653590 11 3.141592516692 3.141592653590 3.141592653590 12 3.141592619365 3.141592653590 3.141592653590 13 3.141592645034 3.141592653590 3.141592653590 14 3.141592651451 3.141592653590 3.141592653590 15 3.141592653055 3.141592653590 3.141592653590 16 3.141592653456 3.141592653590 3.141592653590 17 3.141592653556 3.141592653590 3.141592653590 18 3.141592653581 3.141592653590 3.141592653590 19 3.141592653588 3.141592653590 3.141592653590 20 3.141592653589 3.141592653590 3.141592653590 21 3.141592653590 3.141592653590 3.141592653590
1e-16 1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1
0 5 10 15 20 25
0 1 2
【多角形近似による
π
の計算のAitken
加速】1 3.000000000000 2 3.105828541230 3 3.132628613281
4 3.139350203047 3.141593134327 5 3.141031950891 3.141592683620
6 3.141452472285 3.141592655466 3.141592653591 7 3.141557607912 3.141592653707 3.141592653590 8 3.141583892148 3.141592653597 3.141592653590 9 3.141590463228 3.141592653590 3.141592653590 10 3.141592105999 3.141592653590 3.141592653590 11 3.141592516692 3.141592653590 3.141592653590 12 3.141592619365 3.141592653590 3.141592653590 13 3.141592645034 3.141592653590 3.141592653590
1e-16 1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1
0 2 4 6 8 10 12
0 1 2
【テイラー展開による
arctan( x )
の近似】1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1
1 10 100 1000 10000 100000
"arctan_06.result_1"
x = 0 . 900 , 0 . 910 , . . ., 1 . 000
● 実習内容
1.
多項式( x − 1) k ( x + 1), k = 2, 3, 4
の解x = 1
に対するニュートン法による近似列を, Romberg 加 速を用いて加速計算しなさい.2.
多項式( x − 1) k ( x + 1), k = 2, 3, 4
の解x = 1
に対するニュートン法による近似列を, Aitken加速 を用いて加速計算しなさい.3.
多角形近似による円周率π
の近似値を求める計算をRomberg
加速を用いて加速計算しなさい.4.
多角形近似による円周率π
の近似値を求める計算をAitken
加速を用いて加速計算しなさい.5.
マチンの公式π = 16 arctan(1 / 5) − 4 arctan(1 / 239)
を用いて,
π
の値の近似値を絶対誤差10 −10 で求めなさい. また,π = 4 arctan(1)
を用いて,π
の値
の近似値を絶対誤差 10 −5 で求めなさい. また,これらの計算を加速計算しなさい.
【注意】
