Fortran Example – Newton’s method to find a root

As a warm-up, study the method.

RootFinder.F90: A driver routine.

  1!!------------------------------------------------------------------
  2!! A Fortran example code for finding a root of a user-defined 
  3!! function f(x) = 0.
  4!! 
  5!! This code is written by Prof. Dongwook Lee for AMS 209.
  6!!
  7!!
  8!! * Two methods of iteration:
  9!!   1. Newton's method
 10!!   2. Modified Newton's method
 11!!
 12!! * This routine is a driver routine which calls subroutines:
 13!!   
 14!!   RootFinder:
 15!!       |- setup_init (from setup_module)
 16!!       |
 17!!       |   / newton_method (from findRootMethod_module)
 18!!       |- |
 19!!       |   \ modified_newton_method (from findRootMethod_module)
 20!!       | 
 21!!       |- output_write (from output_module)
 22!!
 23!!------------------------------------------------------------------
 24
 25
 26program RootFinder
 27
 28
 29!! Begining of the real implementation of the driver
 30!! define usages of module variables and subroutines 
 31  use setup_module,          only : setup_init, threshold, maxIter, methodType,xInit
 32  use findRootMethod_module, only : newton_method, modified_newton_method
 33  use output_module,         only : output_write
 34
 35
 36!! Always start with implicit none after use module section
 37  implicit none
 38
 39!! These are local variables whose scopes are within RootFinder.F90
 40  integer :: nIter
 41  real    :: x, xNew, residual, f
 42
 43!! Define allocatable arrays
 44  real, allocatable, dimension(:) :: x_history,f_history,res_history
 45!!$  real, dimension(:,:), allocatable :: dummyVar
 46
 47!! Initial values for residual and number of iteration
 48  residual = 1.e10
 49  nIter = 1
 50
 51!! Call setup_init to initialize the runtime parameters
 52  call setup_init()
 53
 54!! allocate arrays with size of a user-defined integer value, maxIter
 55  allocate(  x_history(maxIter))
 56  allocate(  f_history(maxIter))
 57  allocate(res_history(maxIter))
 58
 59
 60!! The first initial guess for Newton's search.
 61!! This is a user-defined value in rootFinder.init
 62  x = xInit
 63
 64!! Start root searches using either 
 65!! (a) Newton's method, or
 66!! (b) modified Newton's method.
 67
 68  if (methodType == 'newton') then
 69    !print*,'first if'
 70!! Keep search iteration until 
 71!! (a) residual is bigger then a user-defined threshold value, and
 72!! (b) iteration number is less than a user-defined maximum iteration number.
 73
 74    do while ((residual > threshold) .and. (nIter < maxIter))
 75
 76!! Search using conventional Newton's method
 77      call newton_method(x,xNew,f,residual)
 78
 79!! Store a newly updated root into an array 
 80      x_history(nIter) = xNew
 81
 82!! Save for the next search iteration
 83      x = xNew
 84
 85!! Store a newly updated funtion value into an array
 86      f_history(nIter) = f
 87
 88!! Store a newly updated residual into an array
 89      res_history(nIter) = residual
 90
 91!! Update iteration number
 92      nIter = nIter + 1
 93    end do
 94
 95  elseif (methodType == 'modified_newton') then
 96    !print*,'second if'
 97    do while ((residual > threshold) .and. (nIter < maxIter))
 98!! Search using a modified Newton's method
 99      call modified_newton_method(x,xNew,f,residual)
100
101!! Store a newly updated root into an array
102      x_history(nIter) = x
103
104!! Save for the next search iteration
105      x = xNew
106
107!! Store a newly updated funtion value into an array
108      f_history(nIter) = f
109
110!! Store a newly updated residual into an array
111      res_history(nIter) = residual
112
113
114!! Update iteration number
115      nIter = nIter + 1
116    end do
117  else
118!! Exit the program if the search method is unknown.
119    print *, "Unknown method type. Please choose either Newton or Modified Newton"
120    print *, "Aborting the program."
121    stop
122  end if
123
124
125!! Prepare to report outputs:
126!! (a) Short summary onto the screen
127  if ((nIter < maxIter) .and. (residual < threshold)) then
128    print *, '+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++'
129    print *, '             Solution Convergence Summary                          '
130    print *, '+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++'
131    print *, 'Your converged solution x = ', xNew
132    print *, 'Solution converged in Nstep=', nIter-1
133    print *, 'Threshold value = ', threshold
134    print *, '+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++'
135  end if
136
137!! (b) File output
138!! You may want to modify the code to include the information on the initial condition too.
139  call output_write(nIter-1,x_history,f_history,res_history)
140
141!! Make sure allocatable arrays are deallocated before exiting the program
142  deallocate(  x_history)
143  deallocate(  f_history)
144  deallocate(res_history)
145
146end program RootFinder

Download this code

  1. Exercise: Add a couple of lines in RootFinder.F90 to show on the screen what the target function and the initial search value are.

setup_module.F90: A setup routine to initialize several runtime parameters which are read in from the parameter file, rootFinder.init.

 1!!------------------------------------------------------------------
 2!! A Fortran example code for finding a root of a user-defined 
 3!! function f(x) = 0.
 4!! 
 5!! This code is written by Prof. Dongwook Lee for AMS 209.
 6!!
 7!! This module has one subroutine which initialize your setup
 8!! by reading in runtime parameters from 'rootFinder.init' file.
 9!! The setup_init subroutine calls read_initFile*** subroutines
10!! that are implemented as subroutines in read_initFile_module. 
11!!
12!!------------------------------------------------------------------
13
14
15module setup_module
16
17!! include a C-type header file:
18!! this is why the file extensions are .F90, instead of .f90
19#include "definition.h"
20
21  use read_initFile_module
22
23  implicit none
24
25  real, save :: xBeg,xEnd
26  real, save :: threshold
27  real, save :: xInit
28  integer, save :: maxIter
29  integer, save :: ftnType
30  character(len=MAX_STRING_LENGTH), save :: runName, methodType
31  integer, save :: multiplicity
32
33contains
34
35  subroutine setup_init()
36
37    implicit none
38
39    call read_initFileChar('rootFinder.init','run_name',runName)
40    call read_initFileChar('rootFinder.init','method_type',methodType)
41    call read_initFileInt ('rootFinder.init','multiplicity', multiplicity)
42    call read_initFileInt ('rootFinder.init', 'ftn_type', ftnType)
43    if (ftnType > 2) then
44      print*,'Not a valid function input type.'
45      print*,'Aborting now!'
46      stop
47    end if
48    call read_initFileInt ('rootFinder.init', 'max_iter', maxIter)
49    call read_initFileReal('rootFinder.init', 'x_beg', xBeg)
50    call read_initFileReal('rootFinder.init', 'x_end', xEnd)
51    call read_initFileReal('rootFinder.init', 'threshold', threshold)
52    call read_initFileReal('rootFinder.init', 'init_guess', xInit)
53    if ((xInit > xEnd) .or. (xInit < xBeg)) then
54      print*,'Not a valid initial guess -- out of domain.'
55      print*,'Choose your initial guess between x_min and x_max. Aborting now!'
56      stop
57    end if
58
59  end subroutine setup_init
60
61
62end module setup_module

Download this code

  1. Exercise: Modify setup_module.F90 to show the runtime parameters on the screen.

read_initFile_module.F90: A file that reads in user defined values in a user provided parameter file, which is rootFinder.init in this example.

  1!!------------------------------------------------------------------
  2!! A Fortran example code for finding a root of a user-defined 
  3!! function f(x) = 0.
  4!! 
  5!! This code is written by Prof. Dongwook Lee for AMS 209.
  6!!
  7!! This module has four subroutines which can read in various
  8!! types of runtime parameter values (real, integer, logical, and character)
  9!! from a user-defined runtime parameter file, 'rootFinder.init'
 10!!
 11!!------------------------------------------------------------------
 12
 13module read_initFile_module
 14
 15#include "definition.h"
 16
 17  implicit none
 18  integer, parameter :: fileLen=50
 19
 20contains
 21
 22  subroutine read_initFileReal(fileName,varName,varValue)
 23
 24    implicit none
 25    character(len=*),intent(IN) :: fileName,varName
 26    real, intent(OUT) :: varValue
 27
 28    integer :: i,openStatus,inputStatus
 29    real :: simInitVars
 30    character(len=MAX_STRING_LENGTH) :: simCharVars
 31    integer :: pos1,pos2
 32
 33    open(unit = 10, file=fileName, status='old',IOSTAT=openStatus,FORM='FORMATTED',ACTION='READ')
 34
 35    do i=1,fileLen
 36      read(10, FMT = 100, IOSTAT=inputStatus) simCharVars
 37      pos1 = index(simCharVars,varName)
 38      pos2 = pos1+len_trim(varName)
 39      if (pos2 > len_trim(varName)) then
 40        read(simCharVars(pos2+1:),*)simInitVars
 41        !print*,varName,len_trim(varName)
 42        !print*,simCharVars
 43        !print*,pos1,pos2,simCharVars(pos2+1:),simInitVars;stop
 44        varValue = simInitVars
 45      endif
 46    end do
 47
 48    close(10)
 49
 50100 FORMAT(A, 1X, F3.1)
 51
 52  end subroutine read_initFileReal
 53
 54
 55
 56  subroutine read_initFileInt(fileName,varName,varValue)
 57
 58    implicit none
 59    character(len=*),intent(IN) :: fileName,varName
 60    integer, intent(OUT) :: varValue
 61
 62    integer :: i,openStatus,inputStatus
 63    integer :: simInitVars
 64    character(len=MAX_STRING_LENGTH) :: simCharVars
 65    integer :: pos1,pos2
 66
 67    open(unit = 11, file=fileName, status='old',IOSTAT=openStatus,FORM='FORMATTED',ACTION='READ')
 68
 69    do i=1,fileLen
 70      read(11, FMT = 101, IOSTAT=inputStatus) simCharVars
 71      pos1 = index(simCharVars,varName)
 72      pos2 = pos1+len_trim(varName)
 73      if (pos2 > len_trim(varName)) then
 74        read(simCharVars(pos2+1:),*)simInitVars
 75        varValue = simInitVars
 76      endif
 77    end do
 78
 79    close(11)
 80
 81101 FORMAT(A, 1X, I5)
 82
 83  end subroutine read_initFileInt
 84
 85
 86  subroutine read_initFileBool(fileName,varName,varValue)
 87
 88    implicit none
 89    character(len=*),intent(IN) :: fileName,varName
 90    logical, intent(OUT) :: varValue
 91
 92    integer :: i,openStatus,inputStatus
 93    logical :: simInitVars
 94    character(len=MAX_STRING_LENGTH) :: simCharVars
 95    integer :: pos1,pos2
 96
 97    open(unit = 12, file=fileName, status='old',IOSTAT=openStatus,FORM='FORMATTED',ACTION='READ')
 98
 99    do i=1,fileLen
100      read(12, FMT = 102, IOSTAT=inputStatus) simCharVars
101      pos1 = index(simCharVars,varName)
102      pos2 = pos1+len_trim(varName)
103
104      if (pos2 > len_trim(varName)) then
105        read(simCharVars(pos2+1:),*)simInitVars
106        varValue = simInitVars
107      endif
108    end do
109
110    close(12)
111
112102 FORMAT(A, 1X, L)
113
114  end subroutine read_initFileBool
115
116
117
118
119  subroutine read_initFileChar(fileName,varName,varValue)
120
121    implicit none
122    character(len=*),intent(IN)  :: fileName,varName
123    character(len=*),intent(OUT) :: varValue
124
125    integer :: i,openStatus,inputStatus
126    character(len=MAX_STRING_LENGTH) :: simInitVars
127    character(len=MAX_STRING_LENGTH) :: simCharVars
128    integer :: pos1,pos2
129
130    open(unit = 13, file=fileName, status='old',IOSTAT=openStatus,FORM='FORMATTED',ACTION='READ')
131
132    do i=1,fileLen
133      read(13, FMT = 103, IOSTAT=inputStatus) simCharVars
134      pos1 = index(simCharVars,varName)
135      pos2 = pos1+len_trim(varName)
136
137      if (pos2 > len_trim(varName)) then
138        read(simCharVars(pos2+1:),*)simInitVars
139        varValue = simInitVars
140      endif
141    end do
142
143    close(13)
144
145103 FORMAT(A, 1X, A)
146
147  end subroutine read_initFileChar
148
149end module read_initFile_module

Download this code

Note

Check the intrinsic functions list to see what index and len_trim functions are. To see how they work, you can try:

print*,index("I like science","science")
print*,len_trim("science ")

rootFinder.init: A runtime parameter file, including a set of runtime parameter values that are read in by read_initFile_module.F90.

 1run_name	'newton' # [char] Specify your output file name, e.g., 'rootFinder_[run_name].dat'
 2method_type	'newton' # [char] Choose a search method between 'newton' and 'modified_newton'
 3
 4x_beg		 -10.0    # [real] Setting up the search domain 
 5x_end		 10.0    # [real] Setting up the search domain
 6
 7max_iter	 10000    # [int]  Maximum number of iteration
 8threshold	 1.e-8    # [real] Threshold value for solution convergence
 9
10ftn_type	 1       # [int]  Types of function -- either 1 or 2
11init_guess	 -0.5   # [real] Initial guess for root search. Users are responsible to pick a good one.
12multiplicity	 4       # [int]  The multiplicity of the root when using the modified newton method

Download this code

  1. Exercise: In the Newton’s root finding algorithm, it is important to choose a reasonable initial search value. Use a simple plotting tool (e.g., GNU plot) to verify your initial value is closer to the root of the target function.

  2. Exercise: Do you need to recompile your code everytime you change your input runtime parameters in rootFinder.init?

findRootMethod_module.F90: A module that holds two different methods of root finder. Two methods are Newton’s and modified Newton’s methods.

 1!!------------------------------------------------------------------
 2!! A Fortran example code for finding a root of a user-defined 
 3!! function f(x) = 0.
 4!! 
 5!! This code is written by Prof. Dongwook Lee for AMS 209.
 6!!
 7!! This module has two subroutines each of which implements
 8!!   1. Newton's method,
 9!!   2. Modified Newton's method
10!!
11!!------------------------------------------------------------------
12
13
14module findRootMethod_module
15
16#include "definition.h"
17
18  use ftn_module, only: ftn_eval, ftn_derivative
19  use setup_module, only: multiplicity,xBeg,xEnd
20
21contains
22
23!! Conventional Newton's method
24  subroutine newton_method(x,xNew,fNew,residual)
25    implicit none
26    real, intent(IN) :: x
27    real, intent(OUT) :: fNew, residual
28    real :: xNew, f, fprime
29
30!! compute function value evaluated at x
31    call ftn_eval(x,f)
32
33!! compute function derivative evaluated at x
34    call ftn_derivative(x,fprime)
35
36!! Exit if f' is near or become zero
37    if (abs(fprime) < 1.e-12) then
38      print *, '[Error: newton_method] Function derivative becomes very close to zero or zero.'
39      print *, 'f=',f, 'df/dx =',fprime
40      print *, 'Aborting now in order to avoid division by zero.'
41      stop
42    end if
43
44!! Algorithm
45    xNew = x - f/fprime
46    fNew = f
47
48!! Search fails if a newly updated value x is out of the search domain
49    if ((xNew < xBeg) .or. (xNew > xEnd)) then
50      print *, '[Error: newton_method] xNew',xNew, 'is out of domain.'
51      print *, 'Failed in search. Aborting now.'
52      stop
53    end if
54
55!! Calculate a new residual
56    residual = abs(xNew - x)
57
58  end subroutine newton_method
59
60
61
62!! Modified Newton's method
63  subroutine modified_newton_method(x,xNew,fNew,residual)
64    implicit none
65    real, intent(IN) :: x
66    real, intent(OUT) :: fNew, residual
67    real :: xNew, f, fprime
68
69!! compute function value evaluated at x
70    call ftn_eval(x,f)
71
72!! compute function derivative evaluated at x
73    call ftn_derivative(x,fprime)
74
75!! Exit if f' is near or become zero
76    if (abs(fprime) < 1.e-12) then
77      print *, '[Error: modified_newton_method] Function derivative becomes very close to zero or zero.'
78      print *, 'f=',f, 'df/dx =',fprime
79      print *, 'Aborting now in order to avoid division by zero.'
80      stop
81    end if
82
83!! Algorithm
84    xNew = x - multiplicity*f/fprime
85    fNew = f
86
87!! Search fails if a newly updated value x is out of the search domain
88    if ((xNew < xBeg) .or. (xNew > xEnd)) then
89      print *, '[Error: modified_newton_method] xNew',xNew, 'is out of domain.'
90      print *, 'Failed in search. Aborting now.'
91      stop
92    end if
93
94!! Calculate a new residual
95    residual = abs(xNew - x)
96
97  end subroutine modified_newton_method
98
99end module findRootMethod_module

Download this code

ftn_module.F90: A module which includes two numerical evaluations. They are the function evaluation and the derivative of function evaluation. These subroutines are used in both Newton’s and modified Newton’s methods.

 1!!------------------------------------------------------------------
 2!! A Fortran example code for finding a root of a user-defined 
 3!! function f(x) = 0.
 4!! 
 5!! This code is written by Prof. Dongwook Lee for AMS 209.
 6!!
 7!! This module has two subroutines each of which computes
 8!!   1. function value evaluated at x,
 9!!   2. function derivative evaluated at x
10!!
11!!------------------------------------------------------------------
12
13module ftn_module
14
15  use setup_module, only : ftnType
16  implicit none
17
18contains
19
20  subroutine ftn_eval(x,f)
21    implicit none
22    real, intent(IN)  :: x
23    real, intent(OUT) :: f
24
25!! User specific function description for root finding
26    if (ftnType == 1) then
27!! (1) f(x) = x + exp(x) + 10/(1+x^2) - 5
28      f = x + exp(x) + 10./(1.+x**2) - 5.
29    elseif (ftnType == 2) then
30!! (2) f(x) = (x-1)log_10(x)
31!!f = (x - 1.) * log10(x)
32      f = (x - 1.) * (x - 2.)
33    end if
34
35  end subroutine ftn_eval
36
37
38
39  subroutine ftn_derivative(x,fprime)
40    implicit none
41    real, intent(IN)  :: x
42    real, intent(OUT) :: fprime
43
44!! User specific function derivative
45    if (ftnType == 1) then
46!! (1) derivative of the first function
47      fprime = 1. + exp(x) - 20.*x/(1.+x**2)**2
48    elseif (ftnType == 2) then
49!! (2) derivative of the second function
50      !fprime = log10(x) + (x - 1.)/x
51      fprime = (x-1.) + (x-2.)
52    end if
53
54  end subroutine ftn_derivative
55
56
57end module ftn_module

Download this code

  1. Exercise: Add a couple of more functions of your choice and test them.

output_module.F90: A routine to write the results to a file.

 1!!------------------------------------------------------------------
 2!! A Fortran example code for finding a root of a user-defined 
 3!! function f(x) = 0.
 4!! 
 5!! This code is written by Prof. Dongwook Lee for AMS 209.
 6!!
 7!! This module has one subroutine which writes an output 
 8!! to a file whose name is also defined by users.
 9!! Users can specify a custom file name, 'runName', in
10!! rootFinder.init
11!!
12!!------------------------------------------------------------------
13
14module output_module
15
16  use setup_module, only : runName
17  implicit none
18
19contains
20
21  subroutine output_write(length,x,f,residual)
22    implicit none
23    integer, intent(IN) :: length
24    real, dimension(:), intent(IN) :: x,f,residual
25    integer :: i
26    character(len=35) :: ofile
27
28!! File name for ascii output
29    ofile = 'rootFinder_'//trim(runName)//'.dat'
30
31!! File open
32    open(unit=20,file=ofile,status='unknown')
33
34!! Write into a file:
35!!   iteration number, search result x, function value f(x), and residual
36    do i=1,length
37      write(20,920)i,x(i),f(i),residual(i)
38    end do
39
40!! Output format specifiers
41920 format(1x, i5, f16.8, 1x, f16.8, 1x, f16.12)
42
43!! File close
44    close(20)
45
46  end subroutine output_write
47
48end module output_module

Download this code

definition.h: A C-style header file which includes some definitions of code parameters.

1#define MAX_STRING_LENGTH 80
2#define newton 1
3#define modified_newton 2

Download this code

  1. Exercise: Do you need to recompile your code everytime you change the values in defninition.h? Without recompiling your code, do you see your code works properly?

  2. Exercise: Try to delete newton and modified_newton from definition.h and see how it works. Do you see any different code behavior?

makefile: A makefile for compilation.

 1########################################################################################
 2#                       Makefile for RootFinding code
 3#                Written by Prof. Dongwook Lee, AMS 209, UCSC
 4########################################################################################
 5
 6FC	= gfortran
 7
 8FFLAGS_DEBUG = -Wall -Wextra -Wimplicit-interface -fPIC -fmax-errors=1 -g -fcheck=all -fbacktrace
 9
10FFLAGS_OPT = -ggdb -O3 -fdefault-real-8 -fdefault-double-8 -ffree-line-length-none -Wuninitialized
11
12EXE_FILE = rootFinder.exe
13
14OBJS  = RootFinder.o \
15	read_initFile_module.o\
16	findRootMethod_module.o  \
17	ftn_module.o \
18	output_module.o \
19	setup_module.o
20
21.PHONY: clean
22########################################################################################
23#COMPILING AND LINKING USING GENERIC SUFFIX RULE FOR F90
24
25$(EXE_FILE) : $(OBJS)
26	@$(FC) $(FFLAGS_OPT) $(OBJS) -o $(EXE_FILE)
27	@echo "code is now linking..."
28
29#LET'S APPLY GENERIC SUFFIX RULE HERE FOR FORTRAN 90
30
31# method 1 using generic suffix rule
32#.SUFFIXES: 
33#.SUFFIXES: .F90 .o
34#.F90.o:
35#	$(FC) $(FFLAGS_OPT) -c $<
36
37# method 2 using inference rule
38%.o: %.F90
39	$(FC) $(FFLAGS_OPT) -c $<
40
41#######################################################################################
42#SOME USEFUL COMMANDS
43clean:
44	@rm -f *.o *.mod *~ rootFinder.exe *.dat
45
46# debug: clean
47debug: FFLAGS_OPT = $(FFLAGS_DEBUG)
48debug: $(EXE_FILE)
49
50#######################################################################################
51#LET'S DEFINE SOME MODULE DEPENDENCIES!
52RootFinder.o: setup_module.o findRootMethod_module.o output_module.o
53
54setup_module.o : read_initFile_module.o
55
56ftn_module.o : setup_module.o
57
58findRootMethod_module.o : ftn_module.o setup_module.o
59
60#######################################################################################

Download this code

  1. Exercise: What do you expect if you type make debug?

  2. Exercise: Modify the code to include information on the intial guess in the output.