Introduction
This code is designed to be used with the codes given in Bracketing.doc .htm.
Newton
The black is the log(tan(x)) near x = p/2. The blue squares are (x1,f1) and (x2,f2) for the points found by Bracketing.doc .htm such that f1´f2 < 0.
The function f can be evaluated at any location c. The problem is to find c such that
Expand the function about c0, the last postion at which it was evaluated.
(1.2) (1.3)
Solve for the c which makes f(c) equal to zero in this approximation.
Then set c0 equal to c and repeat. The convergence of (1.4) depends on eventually finding a c0 for which f(c0)=0. If this happens (1.1) has been solved. The use of the exact first derivative in (1.3) can speed convergence, but is not necessary.
A derivative that is too large will in general merely slow convergence, while one that is too small may destroy convergence. Note that the use of the actual first derivative at the blue point to the left in figure 1 will make the first estimate of x completely outside the figure. One does somewhat better to estimate the first derivative as
(1.5)
Then in solving
fun(x)-log(1270) = 0, the first red circle is the first
The actual problem frequently is to find the x for which tan(x)-1270 = 0. In this case the value of 1013 at the large x value completely dominates the numerical derivative estimate and all of the red circles are nearly equal to the blue square to the left.
Iteratively
(1.6)
Or
(1.7)
Aitkins.htm .doc. The convergence of
Newton
The complication in this routine is the necessity of keeping f(xlower)*f(xupper)<0
SUBROUTINE NEWAIT(X1,F1,X2,F2,FUN1,NCALL)
C THIS ROUTINE USES (FUNCP-FUNCM)/(XP-XM) TO ESTIMATE THE
C DERIVATIVE IN THE VICINITY OF FUNC = 0. FUNCP*FUNCM < 0.
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION XP(4),XM(4)
NP=0
NM=0
NCALL=0
5 CONTINUE
TEST=ABS((X2-X1)/MAX(MIN(ABS(X2),ABS(X1)),1D-9))
PRINT*,' TEST = ',TEST
IF(TEST.LT.1D-14)RETURN ß normal return
IF(F1.EQ.0)THEN
X2=X1
F2=0
RETURN
ENDIF
IF(F2.EQ.0)THEN
X1=X2
F1=0
RETURN
ENDIF
FP=(F2-F1)/(X2-X1)
IF(ABS(F1).GT.ABS(F2))THEN ß these should both produce the same result
X=X1-F1/FP
ELSE
X=X2-F2/FP
ENDIF
7 CONTINUE
XLOWER=MIN(X1,X2) ß derivatives too small can go outside the region “known” to have a solution
XUPPER=MAX(X1,X2)
IF(X.LE.XLOWER.OR.X.GE.XUPPER)X=(X1+X2)/2 ß defaults to halving
FT=FUN1(X)
print'(a,4E15.6)','X1 x X2, ft',X1,x,X2,ft
read(*,*)itest ß remove in pactice
NCALL=NCALL+1
IF(NCALL.GT.500)STOP' NCALL LIMIT IN NEWAIT ' ß There needs to be some limit
IF(F1*FT.GT.0)THEN
NP=NP+1
NM=-1
IF(NP.GE.1)XP(NP)=X
F1=FT
X1=X
IF(NP.EQ.3)THEN
C AITKINS'S EXTRAPOLATION ß three iterative steps deserve extrapolation
PRINT*,' AITKINS NP '
XN=XP(3)*XP(1)-XP(2)*XP(2)
XD=XP(3)+XP(1)-2*XP(2)
X=XP(3)
IF(XD.NE.0)X=XN/XD
GOTO 7
ENDIF
IF(NP.EQ.4)THEN ß sometimes just halving is best
PRINT*,' NP = 4'
X=(X1+X2)/2
NP=-1
GOTO 7
ENDIF
GOTO 5
ELSE
NM=NM+1
NP=-1
IF(NM.GE.1)XM(NM)=X
F2=FT
X2=X
IF(NM.EQ.3)THEN
C AITKINS'S EXTRAPOLATION ß negative version of the steps above
PRINT*,' AITKINS NM '
XN=XM(3)*XM(1)-XM(2)*XM(2)
XD=XM(3)+XM(1)-2*XM(2)
X=XM(3)
IF(XD.NE.0)X=XN/XD
GOTO 7
ENDIF
IF(NM.EQ.4)THEN
PRINT*,' NM = 4'
X=(X1+X2)/2
NM=-1
GOTO 7
ENDIF
GOTO 5
ENDIF
END
cpp/anewai.c.
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
double anewai(double *x1,double *f1,double *x2, double *f2,double (*fun)(double))
{double xp[3],xm[3],x,xn,xd,xlower,xupper,ft,test,fp;
int np=0,nm=0,ncall=0,itest;
ncall=0;
loop5:
test=fabs((*x2-*x1)/max(min(fabs(*x2),fabs(*x1)),1e-9));
printf(" test = %lg \n",test);
if(test<1e-14)return (*x1+*x2)/2;
if(*f1==0)return *x1;
if(*f2==0)return *x2;
fp=(*f2-*f1)/(*x2-*x1);
if(fabs(*f1)>fabs(*f2))x=*x1-*f1/fp;
else x=*x2-*f2/fp;
loop7:
xlower=min(*x1,*x2);
xupper=max(*x1,*x2);
if(x<=xlower||x>=xupper)x=(*x1+*x2)/2;
ft = fun(x);
ncall=ncall+1;
if(*f1*ft > 0){
np = np + 1;
nm=0;
*f1=ft;
*x1=x;
if(np>0 && np<4 )xp[np]=x;
if(np == 3){
// Aitkin's extrapolation positive
xn = xp[3]*xp[1]-xp[2]*xp[2];
xd = xp[3]+xp[1]-2*xp[2];
x=xp[3];
if(xd != 0)x=xn/xd;
goto loop7;}
if(np == 4){
x=(*x1+*x2)/2;
np=-1;
goto loop7;}
goto loop5;}
else{
nm =nm + 1;
*f2=ft;
*x2=x;
np=0;
if(nm>0 && nm<4 )xp[nm]=x;
if(nm == 3){
// Aitkin's extrapolation negative
xn = xm[3]*xm[1]-xm[2]*xm[2];
xd = xm[3]+xm[1]-2*xm[2];
x=xm[3];
if(xd != 0)x=xn/xd;
goto loop7;}
if(nm==4){
x=(*x1+*x2)/2;
nm=-1;
goto loop7;}
goto loop5;}}
Newton-Raphson (Press 9.6)
In multi-dimensions there are M functions to be solved for M coordinate values. (These frequently arise from the derivatives of a c2(c) with respect to cm ..\optimization\Extremal.doc .htm)
Expand about c0
Solve for the cm’s for which all fm are zero in this approximation
(2.3)
Or
Direct solution of (2.4) for the vector (cm-cm,0) is possible. If the array has special properties, such as Tridiagonal-abc.doc or ConstCoeffs.doc, this may be faster, easier, and more reliable than inverting the matrix.
After a first section on solving (2.4) by Gauss elimination, matrix inversion is considered in GaussJ.doc .htm. More details on matrix inversion including FFT inversion for diagonal matrices are in Welcome.htm#MatrixInversion.
Multiply 2.3 by the inverse of the matrix and sum over m.
(2.5)
So that
(2.6)
Or
(2.7)
More
1. The matrix may not be invertible because there is no set of c values for which the equations in (2.1) are satisfied. – This does not occur when the equations are those that minimize c2(c).
2. The expansion in (2.2) is truncated too soon to give a reasonable estimate of the expansion. When the derivatives are too small, the inverse is too large and c values shoot everywhere. This is “solved” by adding a diagonal element to the partial derivative array. As with the one dimensional example, this slows, but does not stop convergence. ..\stanfit3.doc
3.
The use of