High Performance Scientific Computing   Coursera Edition

#### Previous topic

Fortran debugging

OpenMP

# Fortran example for Newton’s method¶

This example shows one way to implement Newton’s method for solving an equation $$f(x)=0$$, i.e. for a zero or root of the function f(x).

See Newton’s method for the square root for a description of how Newton’s method works.

These codes can be found in $UWHPSC/codes/fortran/newton. Here is the module that implements Newton’s method in the subroutine solve:   1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 !$UWHPSC/codes/fortran/newton/newton.f90 module newton ! module parameters: implicit none integer, parameter :: maxiter = 20 real(kind=8), parameter :: tol = 1.d-14 contains subroutine solve(f, fp, x0, x, iters, debug) ! Estimate the zero of f(x) using Newton's method. ! Input: ! f: the function to find a root of ! fp: function returning the derivative f' ! x0: the initial guess ! debug: logical, prints iterations if debug=.true. ! Returns: ! the estimate x satisfying f(x)=0 (assumes Newton converged!) ! the number of iterations iters implicit none real(kind=8), intent(in) :: x0 real(kind=8), external :: f, fp logical, intent(in) :: debug real(kind=8), intent(out) :: x integer, intent(out) :: iters ! Declare any local variables: real(kind=8) :: deltax, fx, fxprime integer :: k ! initial guess x = x0 if (debug) then print 11, x 11 format('Initial guess: x = ', es22.15) endif ! Newton iteration to find a zero of f(x) do k=1,maxiter ! evaluate function and its derivative: fx = f(x) fxprime = fp(x) if (abs(fx) < tol) then exit ! jump out of do loop endif ! compute Newton increment x: deltax = fx/fxprime ! update x: x = x - deltax if (debug) then print 12, k,x 12 format('After', i3, ' iterations, x = ', es22.15) endif enddo if (k > maxiter) then ! might not have converged fx = f(x) if (abs(fx) > tol) then print *, '*** Warning: has not yet converged' endif endif ! number of iterations taken: iters = k-1 end subroutine solve end module newton 

The solve routine takes two functions f and fp as arguments. These functions must return the value $$f(x)$$ and $$f'(x)$$ respectively for any input x.

A sample test code shows how solve is called:

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 ! $UWHPSC/codes/fortran/newton/test1.f90 program test1 use newton, only: solve use functions, only: f_sqrt, fprime_sqrt implicit none real(kind=8) :: x, x0, fx real(kind=8) :: x0vals(3) integer :: iters, itest logical :: debug ! set to .true. or .false. print *, "Test routine for computing zero of f" debug = .true. ! values to test as x0: x0vals = (/1.d0, 2.d0, 100.d0 /) do itest=1,3 x0 = x0vals(itest) print *, ' ' ! blank line call solve(f_sqrt, fprime_sqrt, x0, x, iters, debug) print 11, x, iters 11 format('solver returns x = ', es22.15, ' after', i3, ' iterations') fx = f_sqrt(x) print 12, fx 12 format('the value of f(x) is ', es22.15) if (abs(x-2.d0) > 1d-14) then print 13, x 13 format('*** Unexpected result: x = ', es22.15) endif enddo end program test1  This makes use of a module functions.f90 that includes the definition of the functions used here. Since $$f(x) = x^2 - 4$$, we are attempting to compute the square root of 4.   1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 !$UWHPSC/codes/fortran/newton/functions.f90 module functions contains real(kind=8) function f_sqrt(x) implicit none real(kind=8), intent(in) :: x f_sqrt = x**2 - 4.d0 end function f_sqrt real(kind=8) function fprime_sqrt(x) implicit none real(kind=8), intent(in) :: x fprime_sqrt = 2.d0 * x end function fprime_sqrt end module functions 

This test can be run via:

$make test1 which uses the Makefile in the same directory:   1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 #$UWHPSC/codes/fortran/newton/Makefile OBJECTS = functions.o newton.o test1.o MODULES = functions.mod newton.mod FFLAGS = -g .PHONY: test1 clean test1: test1.exe ./test1.exe test1.exe: $(MODULES)$(OBJECTS) gfortran $(FFLAGS)$(OBJECTS) -o test1.exe %.o : %.f90 gfortran $(FFLAGS) -c$< %.mod: %.f90 gfortran $(FFLAGS) -c$< clean: rm -f *.o *.exe *.mod