« 平温な一日 | トップページ | ゴキブリ:1 »

2015年2月27日 (金)

√2 の話:その28 割線法(Secant Method:セカント法)

★忘れ物を思い出しました。
この方法、ある意味で「はさみうち法」とペアでお話すべきでした。
Hasamiutisiki
これがはさみうち法の考え方。
aとbが真値を挟むように設定します。
f(a)とf(b)を結ぶ直線が軸と交わる点を新たなaとして真値に接近していきます。
Kassen
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回あたりで収束。
なかなかいいですね。

★微分を学んでいなくても、この方法なら理解できそう。
実際には「差分」という考え方を実行しているのですけどね。
平方根の計算アルゴリズムとしてはかなりよい方法のようです。

« 平温な一日 | トップページ | ゴキブリ:1 »

理科おじさん」カテゴリの記事

コメント

コメントを書く

(ウェブ上には掲載しません)

« 平温な一日 | トップページ | ゴキブリ:1 »

2017年11月
      1 2 3 4
5 6 7 8 9 10 11
12 13 14 15 16 17 18
19 20 21 22 23 24 25
26 27 28 29 30    
サイト内検索
ココログ最強検索 by 暴想
無料ブログはココログ