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
|
- 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
|
- 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
|
- 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.
- 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
|
- 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
|
- 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? - Exercise: Try to delete
newton
andmodified_newton
fromdefinition.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
#######################################################################################
|
- Exercise: What do you expect if you type
make debug
?