# 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`?