二分法

二分法は方程式 f(x)=0 の解を求める計算法です。f(x0) < 0 、f(x1) >　0 なる値x0,x1がわかっていることが必要になります。

program bisectin(input,output);
const
   xneg0 = 0.0; xpos0 = 2.0; eps = 1.0e-11;
var xneg, xpos, xmid : real;
   i : integer;
begin
   i := 0;
   xneg := xneg0; xpos := xpos0;
   while abs(xneg - xpos) > eps do
   begin
      i := i + 1;
      xmid := (xpos + xneg)/2.0;
      if (xmid * xmid * xmid - 2) > 0 then xpos := xmid
      else xneg := xmid
   end;
   writeln('x = ', xneg);
   writeln('i = ', i)
end.

このプログラムの停止性は、xnegとxposの差の絶対値が必ず小さくなっていることからわかる。

実行結果

sino% ./bisection
x =  1.2599210498884E+0000
i =         38

ニュートン法

ニュートン法は、方程式f(x)=0を解く計算法である。

方程式f(x)=0を解くこととします。このとき関数fの導関数をf'とします。適切な初期値 x = x0 を選択して、次の関係式

                f(xn)
    xn+1 = xn - -----
                f'(xn)

によりx0,x1,x2,x3,x4,...を次々に求めていきます。そうするとこの数列が方程式の解に近づいていきます。

２の立方根を求めるためには

f(x) = x3 - 2

f'(x) = 3x2

とすればよく、

xn+1 = (2/3)(xn + (1/xn2))

によりx_n を求めていけばよい。

program newton(input,output);
const eps = 1.0e-11;
var x, x1 : real;
   i : integer;
begin
   x := 1.0; i := 0;
   repeat
      i := i + 1;
      x1 := x;
      x := (2.0/3.0) * (x1 + 1.0/sqr(x1))
   until abs(x - x1) < eps;
   writeln('x = ', x);
   writeln('i = ', i)
end.

sino% pc newton.p
sino% ./newton
x =  1.2599210498949E+0000
i =          6

このプログラムの停止性は数学的な議論しないとわからない。

ここでrepeat文の説明をする。

repeat <文> ; ... ; <文> until <論理式>

意味は、次の通り。

    1.<文>を左から順番に実行する。

    2.<論理式>を計算する。計算結果がtrueならば終了。そうでなければ、1に戻る。