Programming details
In order to make these routines work it is necessary to pass variables to a C function and have them changed and returned to the calling function. The details needed for this are in panda2 and Passing function names in C and Fortran. The Fortran code described here is in fbnewait.zip. The C code is in cbnewait.zip.
Bracketing
The goal is to solve the equation (1.1)
Bracketing in which, one finds a value of x such that g(x) < 0 and another one such that g(x) > 0 is the beginning of the solution process. Press’s method for this is discussed in the section on Bracketing and Bisection 9.1 in the 1986 edition. The test codes to do this, brackbli.for and brackbli.c are described in #BrackBli below. In these routines the subroutine/function BRACK(FUNC,BEG,ENDC,X1,X2,F1,F2,SUCCES)
tries the beginning and end points and then a succession of intermediate points in the hopes of finding a pair x1, x2, for which func(x1)*func(x2) < 0.
The test function
This is the test function for the fortran zip. The test function for the C zip is the Lennard Jones potential described below. cpp/funNG.c is a narrow Gaussian, while the tangent function in C is cpp\funTan.c. Beware of the sign change for the entry ftan, that occurs in Fortran, but not in C.
The classic transcendental equation is
(1.2)
Or
(1.3)
The tan(x) between -p/2 and p/2 is considered in The tan function.doc .htm. This is the natural range.
Figure 1 The tan(x) for x near p/2.
Near ±p/2 the tan goes to ±¥. Numerically p/2 is accurate to only 16 digits, so that the infinity becomes 1016. In fact there is almost no accuracy for x > 1010.
Brack algorithm
Figure 2 Points between -1 and 1 for tan(x)
The above algorithm is designed to give the closest spacing to the input x1 location, since it is presumed that the routine will be used repeatedly to find values close to existing ones.
BrackBli
The function can look very different from that shown above. The traditional function in Robmin is the predicted value of chi-square divided by the original value, c2p/c20, as a function of the log of the Marquardt parameter, l. As l à ¥, this function has the value one for a very large range. As l à -¥, the function becomes very large owing to the inability to solve the equations for the proper parameters. A figure showing this looks very similar to figure 4 at the end of Bli.doc .htm. The codes cpp\brackbli.c and for\brackbli.for, included in the relevant zip files at the beginning of this document, include the bli in the event that the first 20 function evaluations fail to find f1´f2<0. The error for each point in this section is
This differs from that used in Bli.doc .htmwhere the 0.001 is 0.1 and the denominator is not squared. This difference comes about because of an extreme desire to find the function near 0 with as few points as possible.
Figure 3 The Lennard Jones potential as reported on failure to solve VLJ+.25 = 0
The flat region to the right represents no improvement in c2. The 1024 corresponds to a bad solution to the Newton-Raphson equations for the desired coefficients. The point at 0.24 is a “best” solution which we wish to find. Note that the error used in this version of BLI does not produce as detailed a plot in the region near zero as that seen in Bli.doc .htm
Solution output when a -0.24 value for VLJ is requested
C:\temp\cbnewait>cbrack
enter the value of fun(x)
-.24
succes = y
before anewai
x1,f1 = 1.22622 0.0323691
x2,f2 = 1.1271 -0.00985097
test = 0.0879476
test = 0.0660715
test = 0.0561708
test = 0.0266975
test = 0.0260039
test = 0.0258989
test = 0.0129416
test = 0.0129391
test = 0.0129389
test = 0.00646944
test = 0.00646944
test = 0.00646944
test = 0.00323472
test = 0.00323472
test = 0.00323472
test = 0.00161736
test = 0.00080868
test = 0.000101085
test = 5.05425e-005
test = 2.52712e-005
test = 3.15891e-006
test = 1.57945e-006
test = 7.89726e-007
test = 9.87158e-008
test = 4.93579e-008
test = 2.46789e-008
test = 3.08487e-009
test = 1.54243e-009
test = 7.71217e-010
test = 9.64022e-011
test = 4.82011e-011
test = 2.41007e-011
test = 3.01258e-012
test = 1.50629e-012
test = 7.53241e-013
test = 9.41551e-014
test = 4.70775e-014
test = 2.36341e-014
test = 3.04956e-015
x1,f1 1.164993050750716 1.38778e-015
x2,f2 1.164993050750713 -2.77556e-017
at 1.16499 fun1 is 7.21645e-016
enter a 1 to stop
Plus Minus infinity
Tan(x) between 0 and p goes to + infinity as x approaches p/2 from below, then suddenly goes to – infinity as x crosses p/2. This means that an algorithm in which the x’s for which f1*f2 are used to find the function, can easily find x1 = x2, but f1 and f2 are not zero and in fact are very large. This is not just an anomaly of the tangent, it is an always possible result of the algorithm.
The codes for\bnewai.for and cpp\bnewai.c call the bracketing and Newton-Aitkins routines. These check for possible infinities and if present use the infinities as the beginning and ending points for new regions that are checked for possible solutions. Details are in The tan function2.doc .htm.
Two solutions
The assumption in bracketing is that the function is similar to that shown in #figure_2 which has only one zero crossing. The general solution when one wants to find other zeroes is to form the function where a is the first solution. In g(x) is is important to keep x at least e away from a which requires a tiny bit of extra coding. In #figure_3, however the Lennard Jones potential which can be used as a model for the function solved in Robmin.doc .htm, it can be seen that there are almost always two solutions, and that the solution of interest is the one to the right with a positive first derivative. The easiest way to get the desired solution is to start near it since the algorithm is designed to look most closely at nearby values of f(x).
The solution of interest has f1 < 0, x2 > x1 and f2 > 0. If the first call to #brack, produces f1>0, a second call can be made starting at r1 and extending to endc. This is discussed in more detail in Robmin.doc#Bracketing .htm..
Test code
The routine for f(x)=0
cpp\funTan.c
#include <stdlib.h> ß enable the routine to compile separately with watcom
#include <stdio.h>
#include <math.h>
double ftan; ß defines the global
double fun1(double x)
{double fun;
fun=tan(x)-ftan;
return fun;}
cpp\funNG.c
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
double ftan;
double fun1(double x)
{double fun,arg;
arg=(x-.707)/.004;
fun=exp(-arg*arg)-ftan;
return fun;}
cpp\fun1.c – used in cpp\cbrack.tgt
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
double ftan;
double fun1(double x)
{double fun,x3,x6,x12;
x=max(x,.01);
x3=x*x*x;
x6=x3*x3;
x12=x6*x6;
fun=1/x12-1/x6-ftan;
return fun;}
for\fun1.for
FUNCTION FUN1(X)
IMPLICIT REAL*8 (A-H,O-Z)
COMMON /PASS/FTAN ßnote that ftan needs to be passed through a common since it is not in the argument list
FUN1=TAN(X)-FTAN
RETURN
END
for\funLJ.for - used in
FUNCTION FUNLJ(X)
IMPLICIT REAL*8 (A-H,O-Z)
COMMON/PASS/FTAN
XT=MAX(1D-12,X)
X3=XT*XT*XT
X6=X3*X3
X12=X6*X6
FUN1=X12-X6-FTAN
RETURN
END
for\tbrack.for
IMPLICIT REAL*8 (A-H,O-Z)
COMMON/PASS/TXVAL ß This common is both here and in main
EXTERNAL FUN1
DATA PI/3.141592653589793D0/,eps/1d-10/
c details about eps are in "The tan function.doc"
c the maximum size of the tan for this code is 1d10.
c the test function is FUN1
5 PRINT*,' ENTER A 1 TO EVALUATE TAN(X) OR 2 TO FIND ATAN(X)'
READ(*,*)ID
PRINT*,' ENTER X'
READ(*,*)X
IF(ID.EQ.1)THEN
TXVAL=0
PRINT*,' FUN1=',FUN1(X)
GOTO 5
ENDIF
TXVAL=-X ß Note the sign change that is not done in C nor in
BEG=-PI/2+eps
ENDC=PI/2-eps
XTEST=BNEWAI(BEG,ENDC,FUN1)
PRINT*,' FUN1 HAS A ZERO AT ',XTEST
END
cpp\tbrack.c
#include <stdlib.h>
#include <stdio.h>
double ftan; ß This global variable eliminates the need for a common
double fun1(double x); ß These function prototypes allow the functions to be defined in separate segments. If main is last in the same file they are not needed.
double bnewai(double beg,double endc,double (*fun)(double)); ßNote the C method for passing a function.
void main()
{double beg,endc,pi=3.141592653589793,ftest,eps=1e-12;
int itest;
printf(" enter the value of tan(x) \n");
scanf("%lg",&ftan);ßno #sign_change for Fortran equivalent
// details about eps are in "The tan function.doc"
beg=-pi/2+eps;
endc=pi/2-eps;
ftest=bnewai(beg,endc,fun1);
printf(" the solution is %lg %lg\n",ftest,fun1(ftest));
printf(" enter a 1 to stop \n"); ß Stops code in the run version of ide
scanf("%d",&itest);
return;}
Calling code
for/bnewai.for
FUNCTION BNEWAI(BEG,ENDC,FUN1)
IMPLICIT REAL*8 (A-H,O-Z)
CHARACTER*1 SUCCES
EXTERNAL FUN1
X1=(BEG+ENDC)/2
CALL BRACK(FUN1,BEG,ENDC,X1,X2,F1,F2,SUCCES)
IF(SUCCES.NE.'Y')THEN
IF(SUCCES.EQ.'F')THEN
PRINT*,' BRACK FOUND A ZERO AT ',X1
ELSE
PRINT*,' CANNOT FIND A ZERO CROSSING BETWEEN ',BEG,' AND',
2 ENDC
STOP 'STOPPED IN BNEWAI'
ENDIF
ENDIF
PRINT*,X1,F1
PRINT*,X2,F2
PRINT*,SUCCES
CALL NEWAIT(X1,F1,X2,F2,FUN1,NCALL) ß
PRINT*,' NCALL = ',NCALL
PRINT*,' AFTER NEWAIT X1,F1=',X1,F1
PRINT*,' AFTER NEWAIT X2,F2=',X2,F2
C *** +- INFINITY See The tan function2.doc for details -
IF(ABS(F1*F2).GT.1D5)THEN
C TEST THE UPPER END VALUE
IF(X2.LT.X1)THEN
TEMP=X1
X1=X2
X2=TEMP
TEMP=F1
F1=F2
F2=TEMP
ENDIF
X1T=X1
F1T=F1
X2T=X2
F2T=F2
C *** UPPER END TEST
X1=X2
F1=F2
X2=ENDC
F2=FUN1(ENDC)
PRINT*,' X1,F1 ',X1,F1
PRINT*,' X2,F2 ',X2,F2
IF(F2*F1.LT.0)THEN
PRINT*,' TESTING UPPER END'
CALL NEWAIT(X1,F1,X2,F2,FUN1,NCALL)
PRINT*,' NCALL = ',NCALL
PRINT*,' AFTER NRAF2 X1,F1=',X1,F1
PRINT*,' AFTER NRAF2 X2,F2=',X2,F2
ENDIF
ENDIF
PRINT*,' AFTER UPPER END TEST '
IF(ABS(F1*F2).GT.1D5)THEN
C TRY LOWER SOLUTION
TEND=X1T
X1=(BEG+99*TEND)/100
CALL BRACK(FUN1,BEG,TEND,X1,X2,F1,F2,SUCCES)
PRINT*,' SUCESS=',SUCCES
PRINT*,' X1, F1 ',X1,F1
PRINT*,' X2, F2 ',X2,F2
IF(SUCCES.EQ.'Y')CALL NEWAIT(X1,F1,X2,F2,FUN1,NCALL)
PRINT*,' NCALL = ',NCALL
PRINT*,' AFTER NRAF2 X1,F1=',X1,F1
PRINT*,' AFTER NRAF2 X2,F2=',X2,F2
ENDIF
IF(ABS(F1*F2).GT.1D5)THEN
C TRY UPPER SOLUTION
TBEG=X2T
X1=(99*TBEG+ENDC)/100
CALL BRACK(FUN1,TBEG,ENDC,X1,X2,F1,F2,SUCCES)
PRINT*,' SUCESS=',SUCCES
PRINT*,' X1, F1 ',X1,F1
PRINT*,' X2, F2 ',X2,F2
IF(SUCCES.EQ.'Y')CALL NEWAIT(X1,F1,X2,F2,FUN1,NCALL)
PRINT*,' NCALL = ',NCALL
PRINT*,' AFTER NRAF3 X1,F1=',X1,F1
PRINT*,' AFTER NRAF3 X2,F2=',X2,F2
ENDIF
C *** END OF +- INFINITY
BNEWAI=.5D0*(X1+X2)
RETURN
END
cpp/bnewai.c
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
double fun1(double x);
char brack(double (*fun)(double),double beg,double endc,double *x1,double *x2,
double *f1,double *f2);
double anewai(double *x1,double *f1,double *x2, double *f2,double (*fun)(double));
double bnewai(double beg,double endc,double (*fun)(double))
{double x1,x2,f1,f2,ftest,tend,tbeg,temp,fend;
int itest;
char succes;
x1=(beg+endc)/2;
succes=brack(*fun,beg,endc,&x1,&x2,&f1,&f2);
printf(" succes = %c \n",succes);
if(succes == 'f')return x1;
if(succes == 'n')
{fprintf(stderr,"brack could find no evidence of a solution \n");
fprintf(stderr,"in the range %lg to %lg \n",beg,endc);
printf(" enter a 1 to stop \n");
scanf("%d",&itest);
exit(2);}
printf("before anewai\n x1,f1 = %lg %lg\n x2,f2 = %lg %lg\n",x1,f1,x2,f2);
ftest=anewai(&x1,&f1,&x2,&f2,fun);
printf("x1,f1 %24.16g %lg \n",x1,f1);
printf("x2,f2 %24.16g %lg \n",x2,f2);
// make x1 smaller of x1, x2 and x2 the larger
if(x2<x1)
{temp=x1;
x1=x2;
x2=temp;
temp=f1;
f1=f2;
f2=temp;
tend=x1;
tbeg=x2;}
// cheap upper end test
fend=fun(endc);
if(fend*f2<0)
{x1=x2;
f1=f2;
x2=endc;
f2=fend;
ftest=anewai(&x1,&f1,&x2,&f2,fun);
printf("x1,f1 %24.16g %lg \n",x1,f1);
printf("x2,f2 %24.16g %lg \n",x2,f2);}
// lower region
if(fabs(f1)*fabs(f2) > 1e5){
printf(" beg tend %lg %lg \n",beg,tend);
x1=(beg+99*tend)/100;
succes=brack(*fun,beg,tend,&x1,&x2,&f1,&f2);
printf(" succes = %c \n",succes);
if(succes == 'y'){
ftest=anewai(&x1,&f1,&x2,&f2,fun);
printf("x1,f1 %24.16g %lg \n",x1,f1);
printf("x2,f2 %24.16g %lg \n",x2,f2);}}
// upper region
if(fabs(f1)*fabs(f2) > 1e5){
x1=(99*tbeg+endc)/100;
printf(" tbeg endc %g1 %g1 \n",tbeg,endc);
succes=brack(*fun,tbeg,endc,&x1,&x2,&f1,&f2);
printf(" succes = %c \n",succes);
if(succes == 'y'){
ftest=anewai(&x1,&f1,&x2,&f2,fun);
printf("x1,f1 %24.16lg %lg \n",x1,f1);
printf("x2,f2 %24.16lg %lg \n",x2,f2);}}
if(fabs(f1)*fabs(f2) > 1e5){
printf(" no solution ");
exit(3);}
return ftest;}
Subroutine brackbli
cpp/brackbli.c
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
char brack(double (*fun)(double),double beg,double endc,double *x1,double *x2,
double *f1,double *f2)
/* brack tries to find values x1 and x2 such that fun(x1)*fun(x2)<0
the return is 'y' if this has been is the case and 'n' if not */
{double x0,xb,xe,fbel,fabv,fbelc,fabvc,errl,fint,errc,errm;
double err[512],xi[512],fi[512];
double alpha=0.001,beta=0.01,div;
int i,itest,np=512,iprnt=1,ibel,iabv,ndo,ile,iadd,ib,ie;
FILE *fp;
x0=*x1;
*f1=(*fun)(x0);
if(*f1 == 0 )
{*x2=*x1;
*f2=*f1;
return 'f';}
xi[20]=x0;
fi[20]=*f1;
fabv=exp(log(100*(endc-x0))/20);
fbel=exp(log(100*(x0-beg))/20);
fbelc=.01*fbel;
fabvc=.01*fabv;
ibel=19;
iabv=21;
for(i=1;i< 21;i++){
*x2=x0-fbelc;
*f2=(*fun)(*x2);
if(*f2 == 0 )
{*x1=*x2;
*f1=*f2;
return 'f';}
if(*f1**f2<0)
{*x1=xi[ibel+1];
*f1=fi[ibel+1];
return 'y';}
xi[ibel]=*x2;
fi[ibel]=*f2;
*x2=x0+fabvc;
*f2=(*fun)(*x2);
if(*f2 == 0 )
{*x1=*x2;
*f1=*f2;
return 'f';}
if(*f1**f2<0)
{*x1=xi[iabv-1];
*f1=fi[iabv-1];
return 'y';}
xi[iabv]=*x2;
fi[iabv]=*f2;
ibel=ibel-1;
iabv=iabv+1;
fbelc=fbelc*fbel;
fabvc=fabvc*fbel;}
// from cbli.zip
ndo=41; /* there are 41 points numbered 0 to 40 */
fp=fopen("test.out","w"); ß debugging, this and next two lines should be removed once satisfied that code works
for(i=0;i<ndo;i++)fprintf(fp,"%24.16lg %24.16lg\n",xi[i],fi[i]);
fclose(fp);
for (i=1;i < ndo-1;++i)
{xb=xi[i-1];
xe=xi[i+1];
fint=fi[i-1]+((xi[i]-xb)/(xe-xb))*fi[i+1];
div=beta*fabs(fi[i]);
if(div<alpha)div=alpha;
errc=fabs(fint-fi[i])/(div*div); ß error definition appropriate for this case #BrackBli
if(i==1)errl=errc;
err[i-1]=(xi[i]-xi[i-1])*(errc+errl);
errl=errc;}
err[ndo-2]=2*(xi[ndo-1]-xi[ndo-2])*errl;
fp=fopen("err.out","w"); ß debugging this and next two lines need to be removes once satisfied that code works
for(i=0;i<ndo-1;i++)fprintf(fp,"%24.16lg %24.16lg\n",xi[i],err[i]);
fclose(fp);
// There are 40 intervals numbered 0 to 39
/* *** FINDING THE LARGEST ERROR 50 in Fortran */
f_largest_err:
errm=err[0];
ile=0;
for (i=1;i<ndo-1;++i)
{if(err[i]>errm){errm=err[i];
ile=i;}}
/* printf(" largest err is %d %lg \n",ile,errm);
scanf(" %d",&itest);*/
/* *** POINT TO BE INSERTED WILL BE ILE+1, THUS CURRENT ILE+1 BECOMES ILE+2 */
for (i=ndo;i>=ile+1;--i)
{xi[i+1]=xi[i];
fi[i+1]=fi[i];
err[i+1]=err[i];}
iadd=ile+1;
xi[iadd]=(xi[ile]+xi[ile+1])/2;
fi[iadd]=(*fun)(xi[iadd]);
if(fi[iadd]**f1<0)
{*f2=fi[iadd];
*x2=xi[iadd];
*x1=xi[iadd+1];
*f1=fi[iadd+1];
fp=fopen("nret.out","w"); ß normal return – this and next two lines need to be removed after debugging is finished
for(i=0;i<ndo;i++)fprintf(fp,"%24.16lg %24.16lg\n",xi[i],fi[i]);
fclose(fp);
return 'y';}
/* *** ERRL FOR ILE-1 AND ERR(ILE-1 TO ILE+1) NEED TO BE RECALCULATED */
if(ile>0)
{i=ile;
xb=xi[i-1];
xe=xi[i+1];
fint=fi[i-1]+((xi[i]-xb)/(xe-xb))*fi[i+1];
errl=fabs(fint-fi[i]);}
else errl=0;
ib=ile-1;
if(ib<1)ib=1;
ie=ile+3;
if(ie>ndo-1)ie=ndo-1;
for(i=ib;i<=ie;++i)
{xb=xi[i-1];
xe=xi[i+1];
fint=fi[i-1]+((xi[i]-xb)/(xe-xb))*fi[i+1];
div=beta*fabs(fi[i]);
if(div<alpha)div=alpha;
errc=fabs(fint-fi[i])/(div*div);
if(i>ile-1)err[i-1]=(xi[i]-xi[i-1])*(errc+errl);
errl=errc;}
if(ile+3==ndo-1)err[ndo-1]=2*(xi[ndo-1]-xi[ndo-2])*errl;
ndo=ndo+1;
if(np > ndo)goto f_largest_err;
fp=fopen("fail.out","w"); ß This is written on failure to find a zero crossing – probably always important.
for(i=0;i<ndo;i++)fprintf(fp,"%24.16lg %24.16lg\n",xi[i],fi[i]);
fclose(fp);
printf(" after finding xi and fi ");
scanf(" %d",&itest);
return 'n';}
for\brackbli.for (in zip above)
SUBROUTINE BRACK(FEXT,BEG,ENDC,X1,X2,F1,F2,SUCCES)
C BRACK TRIES TO FIND VALUES X1 AND X2 SUCH THAT FEXT(X1)*FEXT(X2)<0
C SUCCES IS RETURNED AS 'Y' IF THIS IS THE CASE AND 'N' IF NOT
C SUCCES IS RETURNED AS 'F' IF FEXT(X1)=0
c This version extends to bli if IBLI=1
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (NP=512)
DIMENSION ERR(NP),XI(NP),FI(NP)
CHARACTER*1 SUCCES
DATA IPRNT/1/
DATA ALPHA,BETA/1D-3,1D-2/
X0=X1
F1=FEXT(X0)
XI(21)=X0
FI(21)=F1
IF(F1.EQ.0)THEN
IF(IPRNT.EQ.1)PRINT*,' IMMEDIATE SUCCESS'
SUCCES='F'
X2=X1
F2=F1
RETURN
ENDIF
ABOUND=ENDC-X0
FABV=EXP(LOG(ABOUND*100)/20)
ABOUND=X0-BEG
FBEL=EXP(LOG(ABOUND*100)/20)
FBELC=FBEL*.01D0
FABVC=FABV*.01D0
IBEL=20
IABV=22
DO I=1,20
X2=X0-FBELC
F2=FEXT(X2)
IF(F2.EQ.0)THEN
IF(IPRNT.EQ.1)PRINT*,' FULL SUCCESS BELOW I = ',I
SUCCES='F'
X1=X2
F1=0
F2=0
RETURN
ELSEIF(F1*F2.LT.0)THEN
IF(IPRNT.EQ.1)PRINT*,' F1*F2<0 BELOW I = ',I
SUCCES='Y'
IF(F2.GT.0)THEN
X1=XI(IBEL+1) ß point 20 above
F1=FI(IBEL+1)
RETURN
ENDIF
XI(IBEL)=X2
FI(IBEL)=F2
X2=X0+FABVC
F2=FEXT(X2)
IF(F2.EQ.0)THEN
IF(IPRNT.EQ.1)PRINT*,' FULL SUCCESS ABOVE I = ',I
SUCCES='F'
X1=X2
F1=0
F2=0
RETURN
ELSEIF(F1*F2.LT.0)THEN
IF(IPRNT.EQ.1)PRINT*,' F1*F2<0 ABOVE I = ',I
X1=XI(IABV-1)
F1=FI(IABV-1)
SUCCES='Y'
RETURN
ENDIF
XI(IABV)=X2
FI(IABV)=F2
IBEL=IBEL-1
IABV=IABV+1
FBELC=FBELC*FBEL
FABVC=FABVC*FABV
ENDDO
C *** TESTING
OPEN(1,FILE='TEST.OUT')
DO I=1,41
WRITE(1,'(2E15.6)')XI(I),FI(I)
ENDDO
CLOSE(1)
IF(IPRNT.EQ.1)PRINT*,' ENTERING BLI SECTION '
C AT THIS POINT THE ORIGINAL 41 TERMS HAVE NOT CROSSED ZERO
NDO=41
DO I=2,NDO-1
XB=XI(I-1)
XE=XI(I+1)
IF(XB.EQ.XE)THEN
FINT=.5D0*(FI(I-1)+FI(I+1))
ELSE
FINT=FI(I-1)+((XI(I)-XB)/(XE-XB))*(FI(I+1)-FI(I-1))
ENDIF
DIV=BETA*ABS(FI(I))
IF(DIV.LT.ALPHA)DIV=ALPHA
ERRC=ABS(FINT-FI(I))/(DIV*DIV)
IF(I.EQ.2)ERRL=ERRC
ERR(I-1)=(XI(I)-XI(I-1))*(ERRC+ERRL)
ERRL=ERRC
ENDDO
ERR(NDO-1)=2*(XI(NDO)-XI(NDO-1))**2*ERRL
C *** FINDING THE LARGEST ERROR
50 CONTINUE
ERRM=ERR(1)
ILE=1
DO I=2,NDO-1
IF(ERR(I).GT.ERRM)THEN
ERRM=ERR(I)
ILE=I
ENDIF
ENDDO
C *** POINT TO BE INSERTED WILL BE ILE+1, THUS CURRENT ILE+1 BECOMES
C *** ILE+2
XI(NDO+1)=XI(NDO)
FI(NDO+1)=FI(NDO)
DO I=NDO-1,ILE+1,-1
XI(I+1)=XI(I)
FI(I+1)=FI(I)
ERR(I+1)=ERR(I)
ENDDO
IADD=ILE+1
XI(IADD)=(XI(ILE)+XI(ILE+1))/2
FI(IADD)=FEXT(XI(IADD))
IF(F1*FI(IADD).LE.0)THEN
X1=XI(IADD+1)
F1=FI(IADD+1)
X2=XI(IADD)
F2=FI(IADD)
SUCCES='Y'
IF(IPRNT.EQ.1)PRINT*,' SUCCESS IN BLI PART IADD= ',IADD
RETURN
ELSEIF(ABS(FI(IADD)).LT.ABS(F1))THEN
X1=XI(IADD)
F1=FI(IADD)
ENDIF
C *** ERRL FOR ILE-1 AND ERR(ILE-1 TO ILE+1) NEED TO BE RECALCULATED
IF(ILE.GT.1)THEN
I=ILE
XB=XI(I-1)
XE=XI(I+1)
FINT=FI(I-1)+((XI(I)-XB)/(XE-XB))*(FI(I+1)-FI(I-1))
ERRL=ABS(FINT-FI(I))
ELSE
ERRL=0
ENDIF
IB=MAX0(2,ILE-2)
IE=MIN0(NDO,ILE+3)
DO I=IB,IE
XB=XI(I-1)
XE=XI(I+1)
FINT=.5D0*(FI(I+1)+FI(I-1))
IF(XE.GT.XB)THEN
FINT=FI(I-1)+((XI(I)-XB)/(XE-XB))*(FI(I+1)-FI(I-1))
ENDIF
DIV=BETA*ABS(FI(I))
IF(DIV.LT.ALPHA)DIV=ALPHA
ERRC=ABS(FINT-FI(I))/(DIV*DIV)
IF(I.EQ.2)ERRL=ERRC
IF(I.GT.ILE-1)ERR(I-1)=(XI(I)-XI(I-1))*(ERRC+ERRL)
ERRL=ERRC
ENDDO
IF(ILE+3.GE.NDO)THEN
ERR(NDO)=2*(XI(NDO)-XI(NDO-1))**2*ERRL
ENDIF
NDO=NDO+1
IF(NP.GT.NDO)GOTO 50
OPEN(1,FILE='FAIL.OUT')
DO I=1,NP
WRITE(1,'(3G15.6)')XI(I),FI(I),ERR(I)
ENDDO
CLOSE(1)
SUCCES='N' ß never found a zero crossing.
RETURN
END
ENDDO
Newton
At this point there exists an {x1, f1} and an {x2,f2} such that f1*f2 < 0. Thus there should be a solution between x1 and x2. The codes cpp/anewai.c and for/newait.for are discussed in Newton.doc .htm.
The Fortran codes are in fbnewait.zip. The C codes are in cbnewait.zip. In addition to the bracketing codes, these
files contain code for using a combination of
Assignment
Use the brack routine to find T such that for T0 = 1.17, d=0.5, p=1.
Write the code such that the uniform grid in T(FA) can be found for FA = (j-1d0/2)*.01 for j between 1 and 100. Send me the code. Plot T(FA) – not the line.