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

! x^5 - 7x^4 + 16x^3 - 8x^2 -16x +16
 
program newton 
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 = 2.40d0
delta = tiny(0.0d0)
n = 1

! 十分な数だけループを回し, 
! x1 - x0 < 予想誤差deitaとなればループを抜ける
do while(n < 10000)  ! ループ変数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(func(x1)) < 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 = x**5.0d0 - 7.0d0*x**4.0d0 + 16.0d0*x**3.0d0 - 8.0d0*x**2.0d0 - 16.d0*x + 16.d0

return

endfunction

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

real(8) :: x

dFunc = 5.0d0*x**4.d0 - 28.0d0*x**3.0d0 +48.0d0*x**2.0d0 - 16.0d0*x -16.0d0

return 

endfunction

