計算機では整数の四則計算の組み合わせで,
より複雑な計算をしているとおもってほぼまちがいない.
たとえば の近似計算を考えてみよう.
この数を近似計算するにはいろいろな方法があるが, 一つの方法は,
は
と
の交点をもとめることである.
一般に
を連続微分可能関数とし,
点 における
の接線の方程式は,
たとえば図 6.2 を見られたい.
この図は,
に対し,
からスタートした
Newton 法の挙動を示している. この場合,
は 1 と 2 という値を
繰り返しとることになる. これは極端な例だが, 極値の周辺で,
不安定な状況が起こることは容易に想像がつくであろう.
初期値さえ適切ならば, Newton 法は高速に の
解を求めることができる.
Newton 法で
の近似値を求めるプログラムは以下のとおり.
def sqrtByNewton(A) { Epsilon = 0.0001; P = 0.0; Q = deval(A); while (!( (Q-P > -Epsilon) && (Q-P < Epsilon))) { P = Q; print(Q); Q = P-(P*P-A)/(2.0*P); } return(Q); }
|
たとえば, ![]() [422] sqrtByNewton(2); 2 1.5 1.41667 1.41422 1.41421 |
Newton 法ほど早くないが, 安定性のある方法として, 2 分法 (bisection method) がある.
が連続関数とするとき,
かつ
なら
となる根
が
区間
に存在するという事実を思いだそう.
この事実をもちいて, 根の存在範囲を狭めていくのが,
2 分法である.
def bisection(A,B,F) { Epsilon = 0.00001; A = deval(A); B=deval(B); if (subst(F,x,A) >= 0 || subst(F,x,B) <= 0) error("F(A)<0 and F(B)>0 must hold."); while (true) { C = (A+B)/2; if ((A-B > -Epsilon) && (A-B < Epsilon)) return(C); FC = subst(F,x,C); if (FC > 0) B = C; else if (FC < 0) A = C; else if (FC > -Epsilon && FC < Epsilon) return(C); } }
|
左のプログラムでは, 根の存在区間の幅が Epsilon よりも
小さくなった時にループを抜けて, その区間の中点の値を返す.
実行すると次のようになる.
[205] bisection(1,2,x^2-2); 1.41421 |
上の問題で, 収束の様子を調べたが, Newton 法でもとまった近似解にどの程度誤差があるかは, 次のようにして理論的に調べることが可能である.
根を ,
をニュートン法にあらわれる
根の近似値の列, 誤差を
以下,
は,
にほぼ等しいということをテイラー展開を用いて証明する.
なので, テイラー展開の公式から,
を満たす定数
が存在して,
したがって, ニュートン法の誤差は, 2 乗, 4 乗, 8 乗,
と急速に減っていく.
たとえば,
の時,
正の根
は約
である.
このとき
の値は
である.
なら, ほぼ
となる.
実際, 精度の高い浮動小数点計算を遂行する次のプログラムで, 誤差の様子をプリントすると 以下のようになる.
setprec(20); def sqrtByNewton(A) { Epsilon = (1/10)^20; P = 0; Q = A; while (!( (Q-P > -Epsilon) && (Q-P < Epsilon))) { Q=eval(Q*exp(0)); print(Q); P = Q; Q = P-(P*P-A)/(2*P); } return(Q); }
|
出力結果
[431] sqrtByNewton(2)$ 2.00000000000000000000000000000000000000 1.500000000000000000000000000000000000000 1.416666666666666666666666666666666666670 1.414215686274509803921568627450980392162 1.414213562374689910626295578890134910125 1.414213562373095048801689623502530243620 [432] eval(2^(1/2)); 1.414213562373095048801688724204443084493 |
関数 setprec で,
何桁の有効数字で計算をおこなうかを指定できる.
setprec(n) の を大きくすると,
高い精度で近似計算をおこなう.
sin などの関数の近似値を高い精度で計算するには,
次のように関数 eval を用いる.
eval(sin(@pi/3));
参考:
1 変数代数方程式の解法については,
Frisco の 最終報告書の Algorithmic Research/Task 3.4: Numerical solving/
3.4.1 Univariate solving
http://www.nag.co.uk/Projects/Frisco/frisco/node7.htm
にいろいろな方法の比較がある.
pari(root, 多項式) は Schoenhage アルゴリズムを用いて根を求めている.