反復演算の簡単な応用例

反復計算を利用する簡単な応用例として、ニュートン法¹⁾のプログラムを考える。ニュートン法は、関数()の導関数()を利用して、反復計算により()=0の解を求める方法である。()=0の真の解をとするとき、その近傍にある近似解をとする。テイラー展開の1次の項まで用いると、()は次のように近似される。

() ()()()

()=であるので、これにより近似解が得られる。

上記の反復計算の過程は、下図のように表せる。()とういう関係を、
点(())を通る接線() + ()() ()で近似する。この近似式から得られる()=0の解、すなわち接線が軸と交わる点を新しい近似解とする。この演算を繰り返して、真の解に十分に近い近似解を求める。ニュートン法では、()・()となるように初期値を選ぶと収束解が得られやすいとされている。

ニュートン法により正の実数の平方根を求めるプログラムを作成する。

ここでは、()=として、()=0の正の解を求めることになる。例題の近似解を求める関係式は次のようになる。

の初期値をとすれば、反復計算により近似解はに近づく。初期値をから新しい近似解を繰り返し求めるプログラムを作成すればよい。

ととの差が十分に小さくなった時点で反復計算を終了し、近似解や反復回数や誤差を出力する。

• ニュートン法によりを求めるプログラム

list1_8.f90

program newton_sqrt
  impplicit none
  real(8) :: x1, x2, a, er, er0=1.0d-6     !er0は許容誤差
  integer :: k, km=100                     !kmは最大反復回数
  write(*,*) 'input a :',                  !入力を促す表示
  read(*,*) a                              !入力されたデータを変数aに格納
  if ( a<= 0.0d0) stop 'a<=0.0d0'          !a>0の時計算実行
  x1 = a                                   !近似解の初期値をaとする
  do k = 1, km                             !最大km回ループを回す
    x2 = x1 - 0.5d0 * ((x1 * x1) - a) / x1 !新しい近似解x2の計算 
    er = abs(x2 - x1)                      !x1とx2の差を計算
    if ( er < er0) exit                    !許容誤差以下ならば計算終了
    x1 = x2                                !x1にx2を入れ替える
  end do                                   !doループ終了 
  write(*,*) 'kai, k, er = ',x2, k, er     !近似解,反復回数,誤差の出力
end program newton_sqrt