! robust_newton.f90 -- ! ! Example belonging to "Modern Fortran in Practice" by Arjen Markus ! ! This work is licensed under the Creative Commons Attribution 3.0 Unported License. ! To view a copy of this license, visit http://creativecommons.org/licenses/by/3.0/ ! or send a letter to: ! Creative Commons, 444 Castro Street, Suite 900, Mountain View, California, 94041, USA. ! ! Robust version of the Newton-Raphson method ! module robust_newton implicit none integer, parameter :: bracketed_root = 1 integer, parameter :: small_value_solution = 2 integer, parameter :: convergence_reached = 3 integer, parameter :: invalid_start_value = -1 integer, parameter :: no_convergence = -2 integer, parameter :: invalid_arguments = -3 contains subroutine find_root( f, xinit, scalex, tolx, smallf, maxevals, result, success ) interface subroutine f( x, indomain, value ) real, intent(in) :: x logical, intent(out) :: indomain real, intent(out) :: value end subroutine f end interface real, intent(in) :: xinit real, intent(in) :: scalex real, intent(in) :: tolx real, intent(in) :: smallf integer, intent(in) :: maxevals real, intent(out) :: result integer, intent(out) :: success real :: eps real :: epsinit real :: epsf real :: fx1 real :: fx2 real :: fxnew real :: fprime real :: x real :: xnew integer :: i integer :: evals logical :: indomain result = 0.0 success = no_convergence ! ! Sanity check ! if ( scalex <= 0.0 .or. tolx <= 0.0 .or. smallf <= 0.0 ) then success = invalid_arguments return endif ! ! Starting value for the stepsize ! epsinit = scalex * sqrt( epsilon(xinit) ) epsf = 100.0 * epsilon(fx1) ! ! Check the initial value ! x = xinit call f( x, indomain, fx1 ) evals = 1 if ( .not. indomain ) then success = invalid_start_value return endif outerloop: & do ! ! Is the function value small enough? ! Then stop ! if ( abs(fx1) <= smallf ) then success = small_value_solution result = x exit endif ! ! Determine the derivative - be careful about the domain ! eps = epsinit epsloop: & do while ( evals < maxevals ) call f( x+eps, indomain, fx2 ) write(*,*) 'Eps: ', eps evals = evals + 1 if ( evals >= maxevals ) exit if ( .not. indomain ) then eps = -eps call f( x+eps, indomain, fx2 ) evals = evals + 1 if ( evals >= maxevals ) exit outerloop if ( .not. indomain ) exit outerloop endif if ( abs(fx2-fx1) < epsf * (abs(fx2)+abs(fx1))/2.0 + smallf ) then eps = 2.0 * eps else exit epsloop endif enddo & epsloop write(*,*) evals, fx1, fx2, eps fprime = (fx2 - fx1) / eps ! ! Determine the next estimate ! newxloop: & do while ( evals < maxevals ) xnew = x - fx1 / fprime call f( xnew, indomain, fxnew ) evals = evals + 1 if ( .not. indomain ) then fx1 = fx1 / 2.0 else exit newxloop endif enddo & newxloop fx1 = fxnew ! ! Have we reached convergence? ! write(*,*) 'Difference: ', abs(xnew-x), ' -- ', scalex * tolx if ( evals < maxevals ) then if ( abs(xnew-x) <= scalex * tolx ) then success = convergence_reached if ( abs(fx1) < smallf ) then success = small_value_solution endif result = xnew write(*,*) 'Result: ',result exit outerloop endif else exit outerloop endif x = xnew write(*,*) evals, x enddo & outerloop ! ! Simply a small value or a bracketed root? ! call f( x - scalex*tolx, indomain, fx1 ) evals = evals + 1 if ( indomain ) then call f( x + scalex*tolx, indomain, fx2 ) evals = evals + 1 if ( indomain ) then if ( fx1 * fx2 <= 0.0 ) then success = bracketed_root endif endif endif write(*,*) 'Number of evaluations: ', evals end subroutine find_root end module module functions implicit none contains subroutine f1( x, indomain, result ) real, intent(in) :: x logical, intent(out) :: indomain real, intent(out) :: result indomain = .true. result = exp(x) + x ! exp(x)-x is also interesting end subroutine subroutine f2( x, indomain, result ) real, intent(in) :: x logical, intent(out) :: indomain real, intent(out) :: result indomain = .true. result = x**2 + 1.0 end subroutine subroutine f3( x, indomain, result ) real, intent(in) :: x logical, intent(out) :: indomain real, intent(out) :: result indomain = .true. ! f3 = sqrt(abs(x)) + 1.0 - cycle: 12.4 -- -130.08 ! f3 = sqrt(abs(x)) + 0.1 - slow divergence ! f3 = sqrt(abs(x)) + 0.001 - cycle: 1.003 -- 1.004 -- 1.003 -- -1.004 result = sqrt(abs(x)) + 0.001 end subroutine subroutine f4( x, indomain, result ) real, intent(in) :: x logical, intent(out) :: indomain real, intent(out) :: result indomain = .true. result = abs(x) ** 0.3333 if ( x > 0.0 ) then result = result - 1.0 else result = -result - 1.0 endif end subroutine subroutine fln( x, indomain, result ) real, intent(in) :: x logical, intent(out) :: indomain real, intent(out) :: result indomain = .true. if ( x <= 0.0 ) then indomain = .false. result = 0.0 else result = log(x) - 1.0 endif end subroutine subroutine fparabola( x, indomain, result ) real, intent(in) :: x logical, intent(out) :: indomain real, intent(out) :: result indomain = .true. result = x ** 2 end subroutine subroutine fparabola2( x, indomain, result ) real, intent(in) :: x logical, intent(out) :: indomain real, intent(out) :: result indomain = .true. result = x ** 2 + 1.0 end subroutine subroutine fsqrt( x, indomain, result ) real, intent(in) :: x logical, intent(out) :: indomain real, intent(out) :: result indomain = .true. result = sqrt(abs(x)) end subroutine subroutine fcos( x, indomain, result ) real, intent(in) :: x logical, intent(out) :: indomain real, intent(out) :: result indomain = .true. result = cos(1.0e4*x) end subroutine end module program test_newton use robust_newton use functions implicit none real :: x real :: root real :: froot integer :: maxiter = 20 integer :: conclusion logical :: indomain write(*,*) 'fln: 0.1' x = 0.1 call find_root( fln, x, 1.0, 1.0e-5, 1.0e-5, maxiter, root, conclusion ) call fln( root, indomain, froot ) write(*,*) root, froot, conclusion write(*,*) 'fln: 10.0' x = 10.0 call find_root( fln, x, 1.0, 1.0e-5, 1.0e-5, maxiter, root, conclusion ) call fln( root, indomain, froot ) write(*,*) root, froot, conclusion write(*,*) 'fparabola: 10.0' x = 10.0 call find_root( fparabola, x, 1.0, 1.0e-5, 1.0e-5, 2*maxiter, root, conclusion ) call fparabola( root, indomain, froot ) write(*,*) root, froot, conclusion write(*,*) 'fparabola2: 10.0' x = 10.0 call find_root( fparabola2, x, 1.0, 1.0e-5, 1.0e-5, maxiter, root, conclusion ) call fparabola2( root, indomain, froot ) write(*,*) root, froot, conclusion write(*,*) 'fsqrt: 1.0' x = 1.0 call find_root( fsqrt, x, 1.0, 1.0e-5, 1.0e-5, maxiter, root, conclusion ) call fsqrt( root, indomain, froot ) write(*,*) root, froot, conclusion write(*,*) 'fcos: 0.0' x = 0.0 call find_root( fcos, x, 1.0e-4, 1.0e-5, 1.0e-5, 100*maxiter, root, conclusion ) call fsqrt( root, indomain, froot ) write(*,*) root, froot, conclusion write(*,*) 'fcos: 0.75e-4*pi' x = 0.75e-4*3.14 call find_root( fcos, x, 1.0e-4, 1.0e-5, 1.0e-5, 100*maxiter, root, conclusion ) call fsqrt( root, indomain, froot ) write(*,*) root, froot, conclusion ! TODO: function with different scale ! ! x = -2.0 ! x = 2.0 works, -2.0 does not ! call find_root( f4, x, 1.0e-5, maxiter, root, conclusion ) ! ! write(*,*) root, f4(root), conclusion ! end program test_newton