3.2 二分法

計算の原理は中間値の定理「連続な関数 $f\colon[a,b]\to\R$ について、 $f(a)$ と $f(b)$ の符号が異なれば、$(a,b)$ に解が少なくとも1つ存在する」と、 その区間縮小法による証明 (付録に書いておいた) に基づく。

$[a,b]$ に解が存在するならば、 2つに分割した区間 $[a,(a+b)/2]$, $[(a+b)/2, b]$ のどちらかに存在する (両方に存在することもある) が、どちらであるか判断できれば、 繰り返すことで区間の幅を半分半分にしていけて、 解を追い詰めることができる。

(詳しいことは事項で解説する -- 暇な時に読んでね。)

以下にサンプル・プログラムを示す。少々長いが、 心臓部分 (プログラム後半部分) は (区間縮小法を理解していれば) 難しくないと信じる。

bisection.BAS

REM bisection.BAS --- 二分法（区間縮小法）で f(x)=0 の近似解を求める
REM 注意！ 1000桁モードにしても超越関数は17桁
REM 超越関数を使わなければ 1000 桁モードで高精度の近似解を得る
REM OPTION ARITHMETIC DECIMAL_HIGH
OPTION ARITHMETIC NATIVE
REM 解きたい方程式 f(x)=0 の定義
FUNCTION F(x)
   LET F=COS(X)-X
END FUNCTION
REM -----------------------------------------------------
LET FMT$="---%.############### "
LET FMT2$=FMT$&FMT$
INPUT PROMPT "左端、右端=": A,B
LET EPS=(B-A)*1.0e-14
REM ------------- 入力した値のチェック ------------------
LET FA=F(A)
LET FB=F(B)
IF FA=0 THEN
   PRINT A;"は解です。"
   STOP
ELSEIF FB=0 THEN
   PRINT B;"は解です。"
   STOP
ELSEIF (FA > 0 AND FB > 0) OR (FA < 0 AND FB < 0) THEN
   PRINT "f(a),f(b)の符号が同じです"
   STOP
END IF
REM ------------ 二分法（区間縮小法）を実行 -------------
LET MAXITR=100
FOR i=1 TO MAXITR
   LET C=(A+B)/2
   LET FC=F(C)
   IF FC=0 THEN
      PRINT "解が見つかりました。"
      PRINT USING FMT$: C
      STOP
   ELSEIF (FA>0 AND FC<0) OR (FA<0 AND FC>0) THEN
      REM 左側[A,C]に解がある
      LET B=C
      LET FB=FC
   ELSE
      REM 右側[C,B]に解がある
      LET A=C
      LET FA=FC
   END IF
   PRINT USING "###": I;
   PRINT USING FMT2$: A,B;
   PRINT FA;FB
   IF B-A < EPS THEN
      PRINT "区間の幅が十分小さくなりました。中点を表示します。"
      PRINT USING FMT$: (A+B)/2
      STOP
   END IF
NEXT I
PRINT "区間の幅は十分小さくなりませんでした。" 
END

実行すると「区間の左端、右端」を尋ねてくる。 例えば 0,1 と答える。

bisection.TXT

左端、右端=0,1
  1   0.500000000000000    1.000000000000000  .377582561890373 -.45969769413186 
  2   0.500000000000000    0.750000000000000  .377582561890373 -1.83111311261791E-2 
  3   0.625000000000000    0.750000000000000  .185963119505218 -1.83111311261791E-2 
  4   0.687500000000000    0.750000000000000  8.53349461524715E-2 -1.83111311261791E-2 
  5   0.718750000000000    0.750000000000000  3.38793724180665E-2 -1.83111311261791E-2 
  6   0.734375000000000    0.750000000000000  7.87472545850132E-3 -1.83111311261791E-2 
  7   0.734375000000000    0.742187500000000  7.87472545850132E-3 -5.19571174375921E-3 
  8   0.738281250000000    0.742187500000000  1.34514975180511E-3 -5.19571174375921E-3 
(中略)
 42   0.739085133214985    0.739085133215212  2.93876034618279E-13 -8.65973959207622E-14 
 43   0.739085133215099    0.739085133215212  1.03583808197527E-13 -8.65973959207622E-14 
 44   0.739085133215156    0.739085133215212  8.54871728961371E-15 -8.65973959207622E-14 
 45   0.739085133215156    0.739085133215184  8.54871728961371E-15 -3.90798504668055E-14 
 46   0.739085133215156    0.739085133215170  8.54871728961371E-15 -1.53210777398272E-14 
 47   0.739085133215156    0.739085133215163  8.54871728961371E-15 -3.44169137633799E-15 
区間の幅が十分小さくなりました。中点を表示します。
   0.739085133215159 

なお、この計算では要求精度 (区間の幅がどこまで小さくなったら反復を停止するか) を (意味は $1\times 10^{-14}$ という意味) としたが、 これは演習に用いている十進BASICの通常の演算精度が、 10 進法 15 桁であることから決めたものである。

やってみよう     $x^2-2=0$ を解くことで、$\sqrt{2}$ を求めてみよ (これは課題8Bの一部である)。