In this post, the downhill simplex method, for the minimization of functions with more than one independent variable, by Nelder and Mead as implemented in F90 in the book of Numerical Recipes in F90 by Press WH, Teukolsky SA, Vetterling WT and Flannery BP is presented.
There are several othe methods for this kind of task and for maximization and the reader is referred to the first edition of Numerical Recipes presented in F77 which contains interesting comments and key references for a better understanding of these topics.
 |
| Fig. 1. Minimization and maximization are part of the major subject optimization |
The method presented here was selected because its application in multidimensional cases.
Contents
_________________________
Some comments on the method
One the reasons this method was chosen is that it does not requires derivatives to determine the minimum of the function but solely evaluations of it. Reading a little about how this method works you will learn that triangles are at the core otf this method and that three guesses are needed to start the process of minimization.
Please, follow this link to read the paper of Nelder JA and Mead R. Comp J 1965 (7) 4 308-313.
This a method that seems to be suitable for cases for which derivatives are not available. This may have the advantage of being not so complicated.
The code
The driver program presented here was adapted into F90 from the F77 version supplied in Numerical Recipes Example Book [Fortran] 1st edition. Also, the subroutines were also created for F77 and are little different from the F90 versions (some arguments in the routine are not needed anymore).
The amoeba subroutine was originaly presented in single precision (SP) but for this post it was modified to work in double precision (DP). This means that the whole program needs to be modified accordingly. Also, the tools in the MyTools.f90 file were also modified.
For this program 4 diferent files are used (as usual),
- main.f90 which is the driver program. This program is supplied below
- PrecisionConstantsStuff.f90 which contains details on numerics for the program to work with double precision. This module is provided below
module PrecisionConstantsStuff
!
INTEGER, PARAMETER:: DP = KIND(1.0D0)
!
INTEGER, PARAMETER :: I4B = SELECTED_INT_KIND(9)
!
end module PrecisionConstantsStuff
- MySubroutines.f90 which is a module with the amoeba routine in F90 as presented in the Numerical Recipes book. You can build your own module file and include the routine supplied in the Numerical Recipes whith not too much trouble
- MyFunctions.f90 which contains, in this case, the working function Eq. (01). This file can bee built as a module to include the function only
- MyTools.f90 which contains minor functions/subroutines used by the amoeba subroutine. If you examine the routine you will need: nrerror, iminloc, etc. These small routines are also supplied in the Numerical Recipes book as well. At this point, you should be aware that some routines needed to be adapted for the amoeba subroutine.
The main program
The adapted code, into F90, is as follows,
program xamoeba
USE MyFunctions
USE PrecisionConstantsStuff
USE MySubroutines
USE MyTools
!
INTEGER(I4B), PARAMETER :: NP = 3, MP = 4
REAL(DP), PARAMETER :: FTOL = 1.0E-16_DP ! This is the desired presicion
REAL(DP) :: P(MP,NP), X(1:NP), Y(MP)
!
! The guess points are supplied below
DATA P/0.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,1.0/
!
DO I = 1,MP
DO J = 1,NP
X(J) = P(I,J)
END DO
! FAMOEB is the function to be minimized
Y(I) = FAMOEB(X)
END DO
!
CALL AMOEBA(P,Y,FTOL,FAMOEB,ITER)
WRITE(*,'(/1X,A,I3)') 'ITERATIONS: ', ITER
WRITE(*,'(/1X,A)') 'VERTICES OF FINAL 3-D SIMPLEX AND'
WRITE(*,'(1X,A)') 'FUNCTION VALUES AT THE VERTICES:'
WRITE(*,'(/3X,A,T11,A,T23,A,T35,A,T45,A/)') 'I', &
'X(I)', 'Y(I)', 'Z(I)', 'FUNCTION'
DO I = 1,MP
WRITE(*,'(1X,I3,4F12.6)') I,(P(I,J),J=1,NP),Y(I)
END DO
WRITE(*,'(1X,A)') 'TRUE MINIMUM IS AT (0.5, 0.6, 0.7)'
end program
As you can see the problem at hand is to determine the minimum of the function,
$\mathbf{fn}(x,y,z)=0.6 - \text{J}_0\left[ (x-0.5)^2 - (y-0.6)^2 - (z-0.7)^2 \right]$ Eq. (01)
which is at $(x,y,z)=(0.5,0.6,0.7)$ as printed at the end of the program. For the program to work, three point guesses, corresponding to the variables and vertices at the base of the triangle, are supplied as a list P. These are: $(1,0,0)$, $(0,1,0)$ and $(0,0,1)$. A fourth point is al so supplied as $(0,0,0)$ which corresponds to the top vertex of the triangle and which is also the aim the program where, nearly, the minimum is located.
The output
Once the program is compiled and run the results come out promptly. The output is not presented neither in the Numerical Recipes nor in the by example book. However, we can easily check if this is working as expected since we already know where the minimum is located.
In our program the output is,
 |
| Fig. 2. Output of the xamoeba program for minimization from the Numerical Recipes in F90 book. |
The program performed 97 iterations but only shows 4. Notice that the fourth iterations is very closed to the exact minimum. However, the evaluation of the function in each iteerations does not change too much. In fact, the evaluation of Eq. (01) at the exact minimum is 4.0.
It is left to the reader to try the program for other cases. Perhaps, a suitable and challenging situation would be that involving experimental data.
Any question? Write in the comments and I shall try to help.
Perhaps the following post may be of interest for you
- Up and running with Numerical Recipes in Fortran 90
- Bracketing roots (zbrac) - Numerical Recipes in F90
- An improvement for bracketing roots (zbrak) - Numerical recipes in F90
- The bisection method (rtbis) - Numerical Recipes in F90
- The false position method (rtfls) - Numerical recipes in F90
- The 4th-order Runge-Kutta method for systems of ode in F90
- Dealing with higher order ODE when using the 4th-order Runge-Kutta mehtod
- Boundary value problem by a shooting method - Case #1 Linear interpolation
- The 4th-order Runge-Kutta method for systems of ode in F90 - NR routine
=========
Ildebrando
No comments:
Post a Comment