Reference: Schaum's Outline of Numerical Analysis by Francis Scheid
This is a function with a multiplicity of maxima and minima.
c optimization see Schaum's Numerical analysis Francis Scheid Chap. 25
c f(x,y)= x1**2 +x2**2
c f(x1,x2)=x1**3 + x2**3-2.*x1**2+3.*x2**2-8.
eq1(x1,x2)=x1-sin(x1+x2)
eq2(x1,x2)=x2-cos(x1-x2)
f(x1,x2)=eq1(x1,x2)**2 +eq2(x1,x2)**2
dfdx1(x1,x2)=(f(x1+epsi,x2)-f(x1,x2))/epsi
dfdx2(x1,x2)=(f(x1,x2+epsi)-f(x1,x2))/epsi
dfdt(t)=(F(x1-(t+epsi)*dfdx1(x1,x2),x2-(t+epsi)*dfdx2(x1,x2) )-
$ F(x1-t*dfdx1(x1,x2),x2-t*dfdx2(x1,x2)))/epsi
dfdt2(t)=(dfdt(t+epsi)-dfdt(t) )/epsi
tol=1.e-5
x1=1.
x2=1.
epsi=1.e-4
print 100,x1,x2,f(x1,x2)
s=f(x1,x2)
print*,'dfdx1,dfdx2=',dfdx1(x1,x2),dfdx2(x1,x2)
do 10 i=1,10
t=.2
c the roots of the first derivative is calculated next it involves both
c first derivative dfdt(t) and the seciond dfdt2(t)
do 20 j=1,10
t=t-dfdt(t)/dfdt2(t)
20 continue
z1=x1-t*dfdx1(x1,x2)
z2=x2-t*dfdx2(x1,x2)
print 100,z1,z2,f(z1,z2)
c compares succesive values of f against the tolerance (tol)
dif=abs(f(z1,z2)-s)
print*,'s,dif,tol=',s,dif,tol
if(dif.le.tol)goto 300
x1=z1-t*dfdx1(z1,z2)
x2=z2-t*dfdx2(z1,z2)
s=f(x1,x2)
print 110 ,dfdx1(z1,z2),dfdx2(z1,z2)
print*, ' '
10 continue
100 format(1x,'z1,z2,f(z1,z2)=',3(3x,e10.3))
110 format(1x,'dfdx1,dfdx2=',2(3x,e10.3))
300 stop
end
RUN
z1,z2,f(z1,z2)= 0.100E+01 0.100E+01 0.823E-02
dfdx1,dfdx2= 0.257105529 0.0755907521
z1,z2,f(z1,z2)= 0.941E+00 0.983E+00 0.278E-03
s,dif,tol= 0.0082269609 0.00794874318 9.99999975E-006
dfdx1,dfdx2= 0.965E-02 -0.320E-01
z1,z2,f(z1,z2)= 0.937E+00 0.994E+00 0.218E-04
s,dif,tol= 8.17394975E-005 5.99335908E-005 9.99999975E-006
dfdx1,dfdx2= 0.435E-02 -0.837E-02
z1,z2,f(z1,z2)= 0.936E+00 0.997E+00 0.196E-05
s,dif,tol= 5.70062184E-006 3.74106708E-006 9.99999975E-006
A relative maxima is obtained in this other RUN with
starting point -.25 , +.25
z1,z2,f(z1,z2)= -0.250E+00 0.250E+00 0.456E+00
dfdx1,dfdx2= 0.601472437 -1.3565402
z1,z2,f(z1,z2)= -0.345E+00 0.464E+00 0.266E+00
s,dif,tol= 0.456359863 0.190110832 9.99999975E-006
dfdx1,dfdx2= 0.321E+00 0.142E+00
z1,z2,f(z1,z2)= -0.205E+00 0.492E+00 0.314E+00
s,dif,tol= 0.246967688 0.0668245926 9.99999975E-006
dfdx1,dfdx2= 0.313E+00 0.334E-01
z1,z2,f(z1,z2)= 0.108E+01 -0.303E+00 0.382E+00
s,dif,tol= 0.363167167 0.0192831401 9.99999975E-006
dfdx1,dfdx2= -0.742E+00 -0.557E+00
z1,z2,f(z1,z2)= -0.651E-02 -0.252E+01 0.325E+01
s,dif,tol= 23.8541565 20.6003609 9.99999975E-006
dfdx1,dfdx2= 0.628E-01 -0.486E+00
z1,z2,f(z1,z2)= -0.149E-02 -0.238E+01 0.322E+01
s,dif,tol= 3.2197032 0.000512188766 9.99999975E-006
dfdx1,dfdx2= 0.770E-01 -0.184E-01
z1,z2,f(z1,z2)= -0.702E-02 -0.237E+01 0.322E+01
s,dif,tol= 3.21895671 0.000107471336 9.99999975E-006
dfdx1,dfdx2= 0.352E-01 -0.532E-02
z1,z2,f(z1,z2)= -0.892E-02 -0.237E+01 0.322E+01
s,dif,tol= 3.2188344 3.30099538E-005 9.99999975E-006
dfdx1,dfdx2= 0.154E-01 -0.607E-02
z1,z2,f(z1,z2)= -0.993E-02 -0.237E+01 0.322E+01
s,dif,tol= 3.21879482 5.23129984E-006 9.99999975E-006