3.4 二分法

計算の原理は、 中間値の定理「連続な関数 $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の一部である)。