!       ニュートン法
!       f(x)=tanh(x)+0.2x+0.3
!       f(x)=0の解
!       伊理テキストP64の表9-1(左・中央)の作成
!       作成者    福士 邦彦
!       2011/06/06 新規作成
!       2011/06/09 最終更新


!---------------メインプログラム-----------------------

program newton1

    implicit none
    integer :: loop  ! ループ回数のカウンター
    real :: xs       ! xの初期値
    real :: xe       ! xの計算値
    real :: lim      ! 計算終了判定の基準値
    real :: hantei   ! 計算終了判定の判定値
    real :: ff       ! 関数f(x)
    real :: dfdx     ! f'(x)の値

    ! xの初期値を代入
    write(*,*) '最初のxの値を入力してください >'
    read(*,*) xs


    ! 計算終了判定の基準値を設定
    lim=epsilon(0.0e0)

    
    ! ループカウンターの初期値を設定
    loop=0
    
    !xeの初期値を設定
    xe=0.0e0

    ! ループ内に入るために最初のhanteiを設定
    hantei=abs((xs-xe)/xs)

    ! abs((xs-xe)/xs)<limとなるまで計算を続けるループ
    do while (hantei>lim)
        
        ! ループカウンターの更新
        loop=loop+1

        ! xeを計算
        xe=xs-ff(xs)/dfdx(xs)

        ! 判定値を更新
        hantei=abs((xs-xe)/xs)

        ! xsを更新
        xs=xe
        
        ! ループ回数とxsとf(xs)の値を表示
        write(*,*) loop,xs,ff(xs)

        !ループ数が10000に到達したら強制的にループを脱出
        if (loop==10000) then
           exit
        endif

    end do
    

end program newton1



!-------------function　f(x)の値の計算-----関数名ff--------------

function ff(x)

    implicit none
    real,intent(in) :: x        ! xの値
    real :: ff       ! f(x)の値

    ! f(x)の計算
    ff=tanh(x)+0.20e0*x+0.30e0

    return

end function ff



!-------------function　f'(x)の値の計算-----関数名dfdx------------

function dfdx(x)

    implicit none
    real,intent(in) :: x        ! xの値
    real :: dfdx     ! f'(x)の値

    ! f'(x)の計算
    dfdx=1.0e0/((cosh(x))**2)+0.20e0

    return

end function dfdx
