√2 の話:その28 割線法(Secant Method:セカント法)
★忘れ物を思い出しました。
この方法、ある意味で「はさみうち法」とペアでお話すべきでした。
これがはさみうち法の考え方。
aとbが真値を挟むように設定します。
f(a)とf(b)を結ぶ直線が軸と交わる点を新たなaとして真値に接近していきます。
aとbを真値より大きい側に設定します。
ニュートン法なら接線を使うわけですが、割線法では上図のように
f(a)とf(b)を結ぶ直線が軸と交わる点を新たなaとし、以前のaをbとして繰り返しながら接近していきます。
「接線」の代わりに、割線を用いるので「割線法」なのです。
かっ‐せん【割線】
〔数〕曲線と二つ以上の点で交わる直線。
広辞苑第六版より引用
そうなると、f(a)とf(b)を結ぶ直線が軸と交わる点を求める計算は、はさみうち法でも割線法でも同じですね。
出発点の与え方が違うだけのようです。
で、はさみうち法のプログラムを改定して割線法にしてみました。
!******************************
!割線法(secant method)
DECLARE EXTERNAL FUNCTION f
INPUT PROMPT "n = ": n
LET a = n
LET b = n + 1
FOR i=1 TO 50
IF (f(b, n) - f(a, n)) <> 0 THEN
LET c = (a * f(b, n) - b * f(a, n)) / (f(b, n) - f(a, n))
LET b = a
LET a = c
PRINT i;
PRINT c;
PRINT c*c
ELSE
EXIT FOR
END IF
NEXT i
END
!******************************
EXTERNAL FUNCTION f(x, n)
LET f = x*x - n
END FUNCTION
!******************************
参考のためにはさみうち法を議論した記事へのリンクをここに置きます。プログラムを比較してください。
http://yamada-kuebiko.cocolog-nifty.com/blog/2015/01/post-c048.html
2015年1月29日 (木)√2 の話:その16:√2の求め方。★ちょっと凝った考え方:その2「はさみうち法」
★今回実は、最初「ゼロ除算エラー」を起こしました。
どうも、「(f(b, n) - f(a, n)) 」←ここがアンダーフローしてしまうようです。
計算機が扱える限界を超えて数値が小さくなってしまうのですね。
で、ゼロ除算が発生する条件になったらループを脱出しておしまい、というように書き換えました。
↓実行結果です。7,8回もループを回ると計算機の限界に来ます。かなり収束が速いようです。
n = 2
1 1.6 2.56
2 1.44444444444444 2.08641975308641
3 1.41605839416058 2.00522137567264
4 1.41423305925716 2.00005514589587
5 1.41421357508149 2.00000003594477
6 1.41421356237318 2.00000000000024
7 1.4142135623731 2.00000000000001
8 1.4142135623731 2.00000000000001
n = 3
1 2.14285714285714 4.59183673469387
2 1.83333333333333 3.3611111111111
3 1.74251497005988 3.03635842088278
4 1.73234719508792 3.00102680432898
5 1.73205170010706 3.00000309184376
6 1.73205080764524 3.00000000026453
7 1.73205080756888 3.00000000000001
8 1.73205080756888 3.00000000000001
n = 10
1 5.71428571428571 32.6530612244897
2 4.27272727272727 18.2561983471074
3 3.44603381014304 11.875149020649
4 3.20309987288071 10.2598487956484
5 3.16401977492572 10.011021136121
6 3.16228882957984 10.0000706416854
7 3.16227766324417 10.000000019453
8 3.16227766016838 10
9 3.16227766016838 10
大きく出て1000でやったら
n = 1000
1 500.749625187406 250750.187125328
2 334.332667332667 111778.332445776
3 201.677079448659 40673.6443749407
4 127.660432152269 16297.1859373041
5 81.2120762910019 6595.40133549551
6 54.4234798538829 2961.915159406
7 39.9588717882051 1596.71143458621
8 33.6365941144189 1131.42046361816
9 31.850879980996 1014.47855556381
10 31.6297910511845 1000.44368194159
11 31.6228018065156 1000.00159409417
12 31.6227766044789 1000.00000017678
13 31.6227766016838 1000
14 31.6227766016838 1000
13回あたりで収束。
なかなかいいですね。
★微分を学んでいなくても、この方法なら理解できそう。
実際には「差分」という考え方を実行しているのですけどね。
平方根の計算アルゴリズムとしてはかなりよい方法のようです。