! 減速ニュートン法を用いて
! f(x) = tanh(x) + 0.2x +0.3
! の解を求める
! 2015/06/15 梅内紫芳

program newton_1

implicit none
integer n        ! ループ変数
real(8) x1       ! 近似解
real(8) x0       ! 初期近似値
real(8) delta    ! 予想誤差
real(8) error    ! 解と近似値の誤差
real(8) func     ! 関数f，サブルーチンから返される
real(8) dfunc    ! fの1階微分，サブルーチンから返される

x0 = 10.0d0       　　　! 初期近似値はある程度検討を付けて設定
delta = epsilon(0.0d0) ! マシンイプシロンで評価する
n = 1

! error が delta より小さくなるまでループ処理
! n は十分な大きさの数に設定
do while(n < 100)

  ! 減速率 μ = 0.3 で計算
  x1 = x0 - 0.3d0 * (func(x0)/dfunc(x0))
  ! 相対誤差を考える
  error = 1 - x1/x0
  ! 収束状況を観察
  write(*,*) 'x0 = ', x0, 'x1 = ', x1
  write(*,*) 'func = ', func(x0), 'dfunc = ', dfunc(x0) 
  write(*,*) 'error = ', error 
 
  ! error が delta より小さくなったらループを抜ける
  if( abs(error) < delta ) then 
    exit
  endif
  
  ! 今求めた近似解x1の値を初期近似値x0として再計算
  x0 = x1
  ! ループ回数のカウント
  n = n + 1

enddo

!近似解とループ回数を出力
write(*,*) 'solution : ', x1
write(*,*) 'n= ', n

endprogram newton_1


! f(x)を求めるサブルーチン
real(8) function func(x)

implicit none
real(8) x

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

endfunction

! f'(x)を求めるサブルーチン
real(8) function dfunc(x)

implicit none
real(8) x

dfunc = 1.0d0 / (cosh(x))**2.0d0 + 0.2d0
return

endfunction
