I am adapting a version of the Newton-Raphson code found in Numerical Recipes in Fortran90 (page 275/572) for solving f(x) = 0. The code (a fortran function) together with a subroutine plus a driver (main) program I have made appear to work with some function f(x). But a different f(x) gives something very different from what I get when I do a paper and pencil calculation.
My ultimate aim is to incorporate this N-R piece of code in a code for solving a stiff ODE. I need to verify it's proper functioning before I can do that.
So, I would highly appreciate help and corrections on my code or rather an example piece of code showing how the numerical recipes version of the code can be implemented.
Here are my pieces of code:
[program rtnewt_driver
! Driver program for testing the N-R iteration code
use nrtype; use nrutil, only : nrerror
use real_type_mod
implicit none
real(kind=wp) :: rtnewt, y, ans
real(kind=wp) :: x, fval, fderiv
real(kind=wp) :: dt, x0, TOL
interface
subroutine funcd(x,fval,fderiv)
use nrtype
use real_type_mod
implicit none
real(kind=wp) :: x, y, dt
real(kind=wp) :: fval,fderiv
end subroutine funcd
end interface
! Specify needed initial values and the tolerance TOL
y = 0.5_wp
dt = 0.01_wp
x0 = y ! initial guess
TOL = 5.0e-10
! Compute and Print result to screen. The function below
! is called with the argument y so that this y-value is
! passed on to x0, the initial guess of the N-R iteration
ans = rtnewt(funcd,x0,TOL)
print*, 'result is:', ans
end program rtnewt_driver]
[function rtnewt(funcd,x0,TOL)
! Function for doing Newton-Raphson iterations
use nrtype; use nrutil, only : nrerror
use real_type_mod
implicit none
real(kind=wp) :: x0, TOL
real(kind=wp) :: rtnewt, prev_rtnewt
integer(I4B), parameter :: maxit = 50
integer(I4B) :: j
real(kind=wp) :: fderiv, dx, fval
interface
subroutine funcd(x,fval,fderiv)
use nrtype
use real_type_mod
implicit none
real(kind=wp) :: x, y, dt
real(kind=wp) :: fval,fderiv
end subroutine funcd
end interface
rtnewt = x0 ! initialising the iterations
do j = 1, maxit
prev_rtnewt = rtnewt ! store previous iterate (x_n)
call funcd(rtnewt, fval, fderiv)
dx = fval/fderiv ! fractional term f(x)/df(x)
rtnewt = rtnewt - dx ! full iterative step
write(*,'(4x, I2, 4x, f22.18, 4x, ES16.8E3)'), j, &
rtnewt, abs(rtnewt - prev_rtnewt)
if ( abs(rtnewt - prev_rtnewt) < TOL) then
print*, 'value of j is:',j
return
end if
end do
call nrerror('rtnewt exceeded maximum iterations')
end function
!-----------------------------------------------------------------
! subroutine "funcd" takes input "x" and returns both the function
! value and the first derivative of the function 'f'.
!-----------------------------------------------------------------
subroutine funcd(x,fval,fderiv)
use nrtype
use real_type_mod
implicit none
real(kind=wp) :: x, y, dt
real(kind=wp) :: fval,fderiv
!real(kind=wp) :: y, dt
fval = x*sin(x) - 1.0
fderiv = x*cos(x) + sin(x)
! fval = 9*dt*(x**3) - 8*dt*(x**2) + (1 - dt)*x - y
! fderiv = 27*dt*(x**2) - 16*dt*x + (1 - dt)
end subroutine funcd]
My ultimate aim is to incorporate this N-R piece of code in a code for solving a stiff ODE. I need to verify it's proper functioning before I can do that.
So, I would highly appreciate help and corrections on my code or rather an example piece of code showing how the numerical recipes version of the code can be implemented.
Here are my pieces of code:
[program rtnewt_driver
! Driver program for testing the N-R iteration code
use nrtype; use nrutil, only : nrerror
use real_type_mod
implicit none
real(kind=wp) :: rtnewt, y, ans
real(kind=wp) :: x, fval, fderiv
real(kind=wp) :: dt, x0, TOL
interface
subroutine funcd(x,fval,fderiv)
use nrtype
use real_type_mod
implicit none
real(kind=wp) :: x, y, dt
real(kind=wp) :: fval,fderiv
end subroutine funcd
end interface
! Specify needed initial values and the tolerance TOL
y = 0.5_wp
dt = 0.01_wp
x0 = y ! initial guess
TOL = 5.0e-10
! Compute and Print result to screen. The function below
! is called with the argument y so that this y-value is
! passed on to x0, the initial guess of the N-R iteration
ans = rtnewt(funcd,x0,TOL)
print*, 'result is:', ans
end program rtnewt_driver]
[function rtnewt(funcd,x0,TOL)
! Function for doing Newton-Raphson iterations
use nrtype; use nrutil, only : nrerror
use real_type_mod
implicit none
real(kind=wp) :: x0, TOL
real(kind=wp) :: rtnewt, prev_rtnewt
integer(I4B), parameter :: maxit = 50
integer(I4B) :: j
real(kind=wp) :: fderiv, dx, fval
interface
subroutine funcd(x,fval,fderiv)
use nrtype
use real_type_mod
implicit none
real(kind=wp) :: x, y, dt
real(kind=wp) :: fval,fderiv
end subroutine funcd
end interface
rtnewt = x0 ! initialising the iterations
do j = 1, maxit
prev_rtnewt = rtnewt ! store previous iterate (x_n)
call funcd(rtnewt, fval, fderiv)
dx = fval/fderiv ! fractional term f(x)/df(x)
rtnewt = rtnewt - dx ! full iterative step
write(*,'(4x, I2, 4x, f22.18, 4x, ES16.8E3)'), j, &
rtnewt, abs(rtnewt - prev_rtnewt)
if ( abs(rtnewt - prev_rtnewt) < TOL) then
print*, 'value of j is:',j
return
end if
end do
call nrerror('rtnewt exceeded maximum iterations')
end function
!-----------------------------------------------------------------
! subroutine "funcd" takes input "x" and returns both the function
! value and the first derivative of the function 'f'.
!-----------------------------------------------------------------
subroutine funcd(x,fval,fderiv)
use nrtype
use real_type_mod
implicit none
real(kind=wp) :: x, y, dt
real(kind=wp) :: fval,fderiv
!real(kind=wp) :: y, dt
fval = x*sin(x) - 1.0
fderiv = x*cos(x) + sin(x)
! fval = 9*dt*(x**3) - 8*dt*(x**2) + (1 - dt)*x - y
! fderiv = 27*dt*(x**2) - 16*dt*x + (1 - dt)
end subroutine funcd]