!  ニュートン法
!  2013/10/10 村橋究理基 新規作成

! x^(n+1) = x^(n) - f(x^(n))/f'(x^(n))
 
program newton1 
implicit none

integer :: n     ! ループ変数
real(8) :: x1    ! 一個後の値 
real(8) :: x0    ! 初期値
real(8) :: delta ! 予想される誤差
real(8) :: error ! x1-x0の値
real(8) :: func  ! 与式
real(8) :: dFunc ! 与式の一次導関数

x0 = huge(0.0d0)
delta = tiny(0.0d0)
n = 1

! 十分な数だけループを回し, 
! x1 - x0 < 予想誤差deltaとなればループを抜ける
do while(n < 100)  ! ループ変数nを十分大きな数まで回す

  ! x1 = x0 - f(x0)/f'(x0)を計算する
  x1 = x0 - func(x0)/dFunc(x0)
  write(*, *) "s0 = ", x0, "x1 = ", x1
  write(*, *) "f = ", func(x0), "f' = ", dFunc(x0)
  ! x1 - x0 を deltaで評価し, 
  ! deltaより小さくなればループから抜ける
  error = x1 - x0
  write(*, *) "error = ", error
  if(abs(error) < delta) then
    exit
  endif 
  ! もしdeltaより小さくならなければn=n+1する
  ! x0 = x1とする
  n = n + 1
  x0 = x1
enddo

! 結果の出力
write(*, *) "解 ", x1
write(*, *) "ループ回数 ", n

endprogram
! xの値に応じてf(x)の値を求めるサブルーチン
real(8) function func(x)
implicit none

real(8) :: x

func = tanh(x) + 0.2d0*x + 0.3d0

return

endfunction

! xの値に応じてf(x)の一回微分f'(x)の値を求めるサブルーチン
real(8) function dFunc(x)
implicit none

real(8) :: x

dFunc = 1.0d0/(cosh(x)) + 0.2d0

return 

endfunction

