Fortran Example – Newton’s method to find a root

As a warm-up, study the method.

RootFinder.F90: A driver routine.

  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
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
!!------------------------------------------------------------------
!! A Fortran example code for finding a root of a user-defined 
!! function f(x) = 0.
!! 
!! This code is written by Prof. Dongwook Lee for AMS 209.
!!
!!
!! * Two methods of iteration:
!!   1. Newton's method
!!   2. Modified Newton's method
!!
!! * This routine is a driver routine which calls subroutines:
!!   
!!   RootFinder:
!!       |- setup_init (from setup_module)
!!       |
!!       |   / newton_method (from findRootMethod_module)
!!       |- |
!!       |   \ modified_newton_method (from findRootMethod_module)
!!       | 
!!       |- output_write (from output_module)
!!
!!------------------------------------------------------------------

   
program RootFinder

  !! include a C-type header file:
  !! this is why the file extensions are .F90, instead of .f90
!#include "definition.h"

  !! Begining of the real implementation of the driver
  !! define usages of module variables and subroutines 
  use setup_module,          only : setup_init, threshold, maxIter, methodType,xInit
  use findRootMethod_module, only : newton_method, modified_newton_method
  use output_module,         only : output_write


  !! Always start with implicit none after use module section
  implicit none

  !! These are local variables whose scopes are within RootFinder.F90
  integer :: nIter
  real    :: x, xNew, residual, f

  !! Define allocatable arrays
  real, allocatable, dimension(:) :: x_history,f_history,res_history
!!$  real, dimension(:,:), allocatable :: dummyVar

  !! Initial values for residual and number of iteration
  residual = 1.e10
  nIter = 1
  
  !! Call setup_init to initialize the runtime parameters
  call setup_init()

  !! allocate arrays with size of a user-defined integer value, maxIter
  allocate(  x_history(maxIter))
  allocate(  f_history(maxIter))
  allocate(res_history(maxIter))


  !! The first initial guess for Newton's search.
  !! This is a user-defined value in rootFinder.init
  x = xInit

  !! Start root searches using either 
  !! (a) Newton's method, or
  !! (b) modified Newton's method.

  if (methodType == 'newton') then
print*,'first if'
     !! Keep search iteration until 
     !! (a) residual is bigger then a user-defined threshold value, and
     !! (b) iteration number is less than a user-defined maximum iteration number.

     do while ((residual > threshold) .and. (nIter < maxIter))

        !! Search using conventional Newton's method
        call newton_method(x,xNew,f,residual)

        !! Store a newly updated root into an array 
        x_history(nIter) = xNew

        !! Save for the next search iteration
        x = xNew

        !! Store a newly updated residual into an array
        res_history(nIter) = residual

        !! Update iteration number
        nIter = nIter + 1
     end do
     
  elseif (methodType == 'modified_newton') then
print*,'second if'
     do while ((residual > threshold) .and. (nIter < maxIter))
        !! Search using a modified Newton's method
        call modified_newton_method(x,xNew,f,residual)

        !! Store a newly updated root into an array
        x_history(nIter) = x

        !! Save for the next search iteration
        x = xNew

        !! Store a newly updated residual into an array
        res_history(nIter) = residual


        !! Update iteration number
        nIter = nIter + 1
     end do
  else
     !! Exit the program if the search method is unknown.
     print *, "Unknown method type. Please choose either Newton or Modified Newton"
     print *, "Aborting the program."
     stop
  end if


  !! Prepare to report outputs:
  !! (a) Short summary onto the screen
  if ((nIter < maxIter) .and. (residual < threshold)) then
     print *, '+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++'
     print *, '             Solution Convergence Summary                          '
     print *, '+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++'
     print *, 'Your converged solution x = ', xNew
     print *, 'Solution converged in Nstep=', nIter-1
     print *, 'Threshold value = ', threshold
     print *, '+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++'
  end if

  !! (b) File output
  call output_write(nIter-1,x_history,f_history,res_history)

  !! Make sure allocatable arrays are deallocated before exiting the program
  deallocate(  x_history)
  deallocate(  f_history)
  deallocate(res_history)
  
end program RootFinder
  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
 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
!!------------------------------------------------------------------
!! A Fortran example code for finding a root of a user-defined 
!! function f(x) = 0.
!! 
!! This code is written by Prof. Dongwook Lee for AMS 209.
!!
!! This module has one subroutine which initialize your setup
!! by reading in runtime parameters from 'rootFinder.init' file.
!! The setup_init subroutine calls read_initFile*** subroutines
!! that are implemented as subroutines in read_initFile_module. 
!!
!!------------------------------------------------------------------


module setup_module

#include "definition.h"
  
  use read_initFile_module
  
  implicit none

  real, save :: xBeg,xEnd
  real, save :: threshold
  real, save :: xInit
  integer, save :: maxIter
  integer, save :: ftnType
  character(len=MAX_STRING_LENGTH), save :: runName, methodType
  integer, save :: multiplicity
  
contains

  subroutine setup_init()

    implicit none

    call read_initFileChar('rootFinder.init','run_name',runName)
    call read_initFileChar('rootFinder.init','method_type',methodType)
    call read_initFileInt ('rootFinder.init','multiplicity', multiplicity)
    call read_initFileInt ('rootFinder.init', 'ftn_type', ftnType)
    if (ftnType > 2) then
       print*,'Not a valid function input type.'
       print*,'Aborting now!'
       stop
    end if
    call read_initFileInt ('rootFinder.init', 'max_iter', maxIter)
    call read_initFileReal('rootFinder.init', 'x_beg', xBeg)
    call read_initFileReal('rootFinder.init', 'x_end', xEnd)
    call read_initFileReal('rootFinder.init', 'threshold', threshold)
    call read_initFileReal('rootFinder.init', 'init_guess', xInit)
    if ((xInit > xEnd) .or. (xInit < xBeg)) then
       print*,'Not a valid initial guess -- out of domain.'
       print*,'Choose your initial guess between x_min and x_max. Aborting now!'
       stop
    end if
    
  end subroutine setup_init

  
end module setup_module
  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
  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
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
!!------------------------------------------------------------------
!! A Fortran example code for finding a root of a user-defined 
!! function f(x) = 0.
!! 
!! This code is written by Prof. Dongwook Lee for AMS 209.
!!
!! This module has four subroutines which can read in various
!! types of runtime parameter values (real, integer, logical, and character)
!! from a user-defined runtime parameter file, 'rootFinder.init'
!!
!!------------------------------------------------------------------

module read_initFile_module

#include "definition.h"
  
  implicit none
  integer, parameter :: fileLen=50

contains

subroutine read_initFileReal(fileName,varName,varValue)

  implicit none
  character(len=*),intent(IN) :: fileName,varName
  real, intent(OUT) :: varValue

  integer :: i,openStatus,inputStatus
  real :: simInitVars
  character(len=MAX_STRING_LENGTH) :: simCharVars
  integer :: pos1,pos2

  open(unit = 10, file=fileName, status='old',IOSTAT=openStatus,FORM='FORMATTED',ACTION='READ')

  do i=1,fileLen
     read(10, FMT = 100, IOSTAT=inputStatus) simCharVars
     pos1 = index(simCharVars,varName)
     pos2 = pos1+len_trim(varName)
     if (pos2 > len_trim(varName)) then
        read(simCharVars(pos2+1:),*)simInitVars
        !print*,varName,len_trim(varName)
        !print*,simCharVars
        !print*,pos1,pos2,simCharVars(pos2+1:),simInitVars;stop
        varValue = simInitVars
     endif
  end do

  close(10)

100 FORMAT(A, 1X, F3.1)

end subroutine read_initFileReal



subroutine read_initFileInt(fileName,varName,varValue)

  implicit none
  character(len=*),intent(IN) :: fileName,varName
  integer, intent(OUT) :: varValue

  integer :: i,openStatus,inputStatus
  integer :: simInitVars
  character(len=MAX_STRING_LENGTH) :: simCharVars
  integer :: pos1,pos2

  open(unit = 11, file=fileName, status='old',IOSTAT=openStatus,FORM='FORMATTED',ACTION='READ')
  
  do i=1,fileLen
     read(11, FMT = 101, IOSTAT=inputStatus) simCharVars
     pos1 = index(simCharVars,varName)
     pos2 = pos1+len_trim(varName)
     if (pos2 > len_trim(varName)) then
        read(simCharVars(pos2+1:),*)simInitVars
        varValue = simInitVars
     endif
  end do
  
  close(11)

101 FORMAT(A, 1X, I5)

end subroutine read_initFileInt


subroutine read_initFileBool(fileName,varName,varValue)

  implicit none
  character(len=*),intent(IN) :: fileName,varName
  logical, intent(OUT) :: varValue

  integer :: i,openStatus,inputStatus
  logical :: simInitVars
  character(len=MAX_STRING_LENGTH) :: simCharVars
  integer :: pos1,pos2

  open(unit = 12, file=fileName, status='old',IOSTAT=openStatus,FORM='FORMATTED',ACTION='READ')
  
  do i=1,fileLen
     read(12, FMT = 102, IOSTAT=inputStatus) simCharVars
     pos1 = index(simCharVars,varName)
     pos2 = pos1+len_trim(varName)

     if (pos2 > len_trim(varName)) then
        read(simCharVars(pos2+1:),*)simInitVars
        varValue = simInitVars
     endif
  end do

  close(12)

102 FORMAT(A, 1X, L)

end subroutine read_initFileBool




subroutine read_initFileChar(fileName,varName,varValue)

  implicit none
  character(len=*),intent(IN)  :: fileName,varName
  character(len=*),intent(OUT) :: varValue

  integer :: i,openStatus,inputStatus
  character(len=MAX_STRING_LENGTH) :: simInitVars
  character(len=MAX_STRING_LENGTH) :: simCharVars
  integer :: pos1,pos2

  open(unit = 13, file=fileName, status='old',IOSTAT=openStatus,FORM='FORMATTED',ACTION='READ')
  
  do i=1,fileLen
     read(13, FMT = 103, IOSTAT=inputStatus) simCharVars
     pos1 = index(simCharVars,varName)
     pos2 = pos1+len_trim(varName)

     if (pos2 > len_trim(varName)) then
        read(simCharVars(pos2+1:),*)simInitVars
        varValue = simInitVars
     endif
  end do

  close(13)

103 FORMAT(A, 1X, A)

end subroutine read_initFileChar

end module read_initFile_module

Note

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

print*,index("I am a boy","boy")
print*,len_trim("boy ")

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
run_name	'newton' # [char] Specify your output file name, e.g., 'rootFinder_[run_name].dat'
method_type	'modified_newton' # [char] Choose a search method between 'newton' and 'modified_newton'

x_beg		 -10.0    # [real] Setting up the search domain 
x_end		 10.0     # [real] Setting up the search domain

max_iter	 1000    # [int]  Maximum number of iteration
threshold	 1.e-8   # [real] Threshold value for solution convergence

ftn_type	 2       # [int]  Types of function -- either 1 or 2
init_guess	 1.5     # [real] Initial guess for root search. Users are responsible to pick a good one.
multiplicity	 4       # [int]  The multiplicity of the root when using the modified newton method
  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
 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
86
87
88
89
90
91
92
93
94
95
96
97
98
99
!!------------------------------------------------------------------
!! A Fortran example code for finding a root of a user-defined 
!! function f(x) = 0.
!! 
!! This code is written by Prof. Dongwook Lee for AMS 209.
!!
!! This module has two subroutines each of which implements
!!   1. Newton's method,
!!   2. Modified Newton's method
!!
!!------------------------------------------------------------------


module findRootMethod_module

#include "definition.h"
  
  use ftn_module, only: ftn_eval, ftn_derivative
  use setup_module, only: multiplicity,xBeg,xEnd
  
contains

  !! Conventional Newton's method
  subroutine newton_method(x,xNew,fNew,residual)
    implicit none
    real, intent(IN) :: x
    real, intent(OUT) :: fNew, residual
    real :: xNew, f, fprime
    
    !! compute function value evaluated at x
    call ftn_eval(x,f)

    !! compute function derivative evaluated at x
    call ftn_derivative(x,fprime)

    !! Exit if f' is near or become zero
    if (abs(fprime) < 1.e-12) then
       print *, '[Error: newton_method] Function derivative becomes very close to zero or zero.'
       print *, 'f=',f, 'df/dx =',fprime
       print *, 'Aborting now in order to avoid division by zero.'
       stop
    end if

    !! Algorithm
    xNew = x - f/fprime
    fNew = f
    
    !! Search fails if a newly updated value x is out of the search domain
    if ((xNew < xBeg) .or. (xNew > xEnd)) then
       print *, '[Error: newton_method] xNew',xNew, 'is out of domain.'
       print *, 'Failed in search. Aborting now.'
       stop
    end if
    
    !! Calculate a new residual
    residual = abs(xNew - x)
    
  end subroutine newton_method



  !! Modified Newton's method
  subroutine modified_newton_method(x,xNew,fNew,residual)
    implicit none
    real, intent(IN) :: x
    real, intent(OUT) :: fNew, residual
    real :: xNew, f, fprime
    
    !! compute function value evaluated at x
    call ftn_eval(x,f)

    !! compute function derivative evaluated at x
    call ftn_derivative(x,fprime)

    !! Exit if f' is near or become zero
    if (abs(fprime) < 1.e-12) then
       print *, '[Error: modified_newton_method] Function derivative becomes very close to zero or zero.'
       print *, 'f=',f, 'df/dx =',fprime
       print *, 'Aborting now in order to avoid division by zero.'
       stop
    end if
    
    !! Algorithm
    xNew = x - multiplicity*f/fprime
    fNew = f
    
    !! Search fails if a newly updated value x is out of the search domain
    if ((xNew < xBeg) .or. (xNew > xEnd)) then
       print *, '[Error: modified_newton_method] xNew',xNew, 'is out of domain.'
       print *, 'Failed in search. Aborting now.'
       stop
    end if

    !! Calculate a new residual
    residual = abs(xNew - x)
    
  end subroutine modified_newton_method
  
end module findRootMethod_module

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
 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
!!------------------------------------------------------------------
!! A Fortran example code for finding a root of a user-defined 
!! function f(x) = 0.
!! 
!! This code is written by Prof. Dongwook Lee for AMS 209.
!!
!! This module has two subroutines each of which computes
!!   1. function value evaluated at x,
!!   2. function derivative evaluated at x
!!
!!------------------------------------------------------------------

module ftn_module

  use setup_module, only : ftnType
  implicit none

contains

  subroutine ftn_eval(x,f)
    implicit none
    real, intent(IN)  :: x
    real, intent(OUT) :: f

    !! User specific function description for root finding
    if (ftnType == 1) then
       !! (1) f(x) = x + exp(x) + 10/(1+x^2) - 5
       f = x + exp(x) + 10./(1.+x**2) - 5.
    elseif (ftnType == 2) then
       !! (2) f(x) = (x-1)log_10(x)
       !!f = (x - 1.) * log10(x)
       f = (x - 1.) * (x - 2.)
    end if

  end subroutine ftn_eval



  subroutine ftn_derivative(x,fprime)
    implicit none
    real, intent(IN)  :: x
    real, intent(OUT) :: fprime

    !! User specific function derivative
    if (ftnType == 1) then
       !! (1) derivative of the first function
       fprime = 1. + exp(x) - 20.*x/(1.+x**2)**2
    elseif (ftnType == 2) then
       !! (2) derivative of the second function
       !fprime = log10(x) + (x - 1.)/x
       fprime = (x-1.) + (x-2.)
    end if

  end subroutine ftn_derivative

  
end module ftn_module
  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
 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
!!------------------------------------------------------------------
!! A Fortran example code for finding a root of a user-defined 
!! function f(x) = 0.
!! 
!! This code is written by Prof. Dongwook Lee for AMS 209.
!!
!! This module has one subroutine which writes an output 
!! to a file whose name is also defined by users.
!! Users can specify a custom file name, 'runName', in
!! rootFinder.init
!!
!!------------------------------------------------------------------

module output_module

  use setup_module, only : runName
  implicit none

contains

  subroutine output_write(length,x,f,residual)
    implicit none
    integer, intent(IN) :: length
    real, dimension(:), intent(IN) :: x,f,residual
    integer :: i
    character(len=35) :: ofile

    !! File name for ascii output
    ofile = 'rootFinder_'//trim(runName)//'.dat'

    !! File open
    open(unit=20,file=ofile,status='unknown')
    
    !! Write into a file:
    !!   iteration number, search result x, function value f(x), and residual
    do i=1,length
       write(20,920)i,x(i),f(i),residual(i)
    end do

    !! Output format specifiers
920 format(1x, i5, f16.8, 1x, f16.8, 1x, f16.12)
    
    !! File close
    close(20)

  end subroutine output_write
  
end module output_module

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

1
2
3
#define MAX_STRING_LENGTH 80
#define newton 1
#define modified_newton 2
  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
 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
########################################################################################
#                       Makefile for RootFinding code
#                Written by Prof. Dongwook Lee, AMS 209, UCSC
########################################################################################

FC	= gfortran

FFLAGS_DEBUG = -Wall -Wextra -Wimplicit-interface -fPIC -fmax-errors=1 -g -fcheck=all -fbacktrace

FFLAGS_OPT = -ggdb -O3 -fdefault-real-8 -fdefault-double-8 -ffree-line-length-none -Wuninitialized

EXE_FILE = rootFinder.exe

OBJS  = RootFinder.o \
	read_initFile_module.o\
	findRootMethod_module.o  \
	ftn_module.o \
	output_module.o \
	setup_module.o

.PHONY: clean
########################################################################################
#COMPILING AND LINKING USING GENERIC SUFFIX RULE FOR F90

$(EXE_FILE) : $(OBJS)
	@$(FC) $(FFLAGS_OPT) $(OBJS) -o $(EXE_FILE)
	@echo "code is now linking..."

#LET'S APPLY GENERIC SUFFIX RULE HERE FOR FORTRAN 90

# method 1 using generic suffix rule
#.SUFFIXES: 
#.SUFFIXES: .F90 .o
#.F90.o:
#	$(FC) $(FFLAGS_OPT) -c $<

# method 2 using inference rule
%.o: %.F90
	$(FC) $(FFLAGS_OPT) -c $<

#######################################################################################
#SOME USEFUL COMMANDS
clean:
	@rm -f *.o *.mod *~ rootFinder.exe *.dat

# debug: clean
debug: FFLAGS_OPT = $(FFLAGS_DEBUG)
debug: $(EXE_FILE)

#######################################################################################
#LET'S DEFINE SOME MODULE DEPENDENCIES!
RootFinder.o: setup_module.o findRootMethod_module.o output_module.o

setup_module.o : read_initFile_module.o

ftn_module.o : setup_module.o

findRootMethod_module.o : ftn_module.o setup_module.o

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