1.4 二分法

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

$[a,b]$ に解が存在するならば、 2つに分割した区間 $\dsp\left[a,\frac{a+b}{2}\right]$ , $\dsp\left[\frac{a+b}{2}, b\right]$ のどちらかに存在します (両方に存在することもある) が、どちらにあるか判断できれば、 繰り返すことで区間の幅を半分半分にして行けて、 解を追い詰めることが出来る、ということです。 (詳しいことは次の小節で -- 暇な時に読んでね。)

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

bisection.BAS

REM bisection.BAS --- 二分法（区間縮小法）で f(x)=0 の近似解を求める
REM 注意！ 1000桁モードにしても超越関数は17桁
REM 超越関数を使わなければ 1000 桁モードで高精度の近似解を得る
REM 解きたい方程式 f(x)=0 の定義
FUNCTION F(x)
   LET F=COS(X)-X
END FUNCTION
REM -----------------------------------------------------
LET FMT$="---%."&REPEAT$("#",15)&" "
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 ------------ 二分法（区間縮小法）を実行 -------------
PRINT USING "  0"&FMT2$: a,b
LET MAXITR=1000
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
   REM 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
  0   0.000000000000000    1.000000000000000 
  1   0.500000000000000    1.000000000000000 
  2   0.500000000000000    0.750000000000000 
  3   0.625000000000000    0.750000000000000 
  4   0.687500000000000    0.750000000000000 
  5   0.718750000000000    0.750000000000000 
  6   0.734375000000000    0.750000000000000 
  7   0.734375000000000    0.742187500000000 
  8   0.738281250000000    0.742187500000000 
(中略)
 42   0.739085133214989    0.739085133215216 
 43   0.739085133215103    0.739085133215216 
 44   0.739085133215160    0.739085133215216 
 45   0.739085133215160    0.739085133215188 
 46   0.739085133215160    0.739085133215174 
 47   0.739085133215160    0.739085133215167 
区間の幅が十分小さくなりました。中点を表示します。
   0.739085133215164 

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

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