Best Linear Interpolator

bli.wpj
TBLI.FOR
bli.for
bli2.for
cbli.wpj
tbli.c
bli.c

The code bli2.for is the same as bli.for. Bli saves temporary results. These will not be unique if a function called by BLI also calls BLI. When this happens, make the second one BLI2 and all will be well. The parameters passed in common by tBli are given in definitions\LennardJones\Lennard Jones Potential.htm

The routine starts with the 4 red points. Points 1 and 3 are interpolated to find a value for point 2. Points 2 and 4 are interpolated to find a value for point 3. The region with the largest value of d times delta x has a point placed in the middle of it. This continues until the maximum number of points is reached or until the maximum interpolation error falls below some tolerance.

The object is to use exactly N function evaluations to get the best estimate of a function possible with N function calls. - The assumption is that each call is relatively time consuming.

Figure 1 Region of a curve showing linear interpolations for points 3 and 6.

The BLI code

The code begins with a small number of points as shown in figure 1. Each point is assigned a relative error based on the difference between its value and a linear interpolation between the neighboring points. The error at the beginning and ending points are taken to be those at the second and second to last points. Each interval is assigned an error based on the absolute value of the sum of the end points of the interval times the size of the interval to a power. The interval with the largest error has a value of the underlying function evaluated at its mid point. In figure 1, this would involve placing a new point in the interval from 2 to 3. The new point changes the error estimates for points 1, 2, new, and 3. All others remain the same.

Ide’s

bli.wpj

cbli.wpj

The FORTRAN and C ide’s superficially differ in that the Lennard-Jones potential is accessed directly in FORTRAN, but is coded into tbli.c.

Set the options to those shown below.

Compilation and running can check for everything.

Full debugging is an art of its own. The line number of the error almost always is sufficient to correct the problem.

TBLI.FOR

IMPLICIT REAL*8 (A-H,O-Z) ß single precision is almost never adequate for physics

DIMENSION XI(:),YI(:),ERR(:) ß watcom has extended f77 to allocatable arrays, f90

code requires allocatable:: XI,YI,ERR – produces error in watcom

CHARACTER*1 ANS ß will be used to finally stop the code

EXTERNAL ALJ ß tells the compiler that ALJ is a function outside this routine

COMMON/PLJ/EPS,SIGMA ß named common is used here to pass parameters to ALJ

RMIN=1.4D0 ß these three lines are described in

SIGMA=RMIN/2**(1D0/6) ß ..\definitions\LennardJones\Lennard Jones Potential.htm

EPS=0.174 ß

PRINT*,'ENTER NP,NP2'

READ(*,*)NP,NP2

IF(NP2.LT.NP)NP2=NP

ALLOCATE(XI(NP2),YI(NP2),ERR(NP2))

PRINT*,' ENTER B,E'

READ(*,*)XI(1),XI(2)

NDO=2 ß This uses half of the first Np values in a uniform grid

CALL BLI(XI,YI,ERR,NP,NDO,ALJ)

IF(NP2.GT.NP)CALL BLI(XI,YI,ERR,NP2,NP,ALJ) ß NP2 can fill in many more values

OPEN(1,FILE='BLI.OUT')

WRITE(1,'(2G14.6)')(XI(I),YI(I),I=1,NP2) ßThis is the results file.

CLOSE(1)

CALL FSYSTEM('GPLOT BLI.OUT') ß Always make it easy to look at the data

READ(*,'(A)')ANS

END

tbli.c

#include <stdio.h>

#include <math.h>

double sigma; ßnames outside the function are the C version of common

double rmin=1.4,eps=0.174;

void *calloc( size_t n, size_t size );

int system(const char *command);

int bli(double xi[],double fi[],double err[],

int np,int ndo,double (*fext)(double));

double alj(double xt) ßthe Cversion of the Lennard-Jones potential

{double xi,xi6,xi12,retval;

extern double sigma,eps; ß tells the code that these are as defined in the file

if(xt < 0.01)xt=0.01;

xi=sigma/xt;

xi=xi*xi*xi;

xi6=xi*xi;

xi12=xi6*xi6;

retval=4*eps*(xi12-xi6);

return retval;}

void main(void)

{double *xi,*fi,*err;

char c1;

extern double sigma,rmin;

int np,np2,ndo=2,i,ibli,itest;

FILE *fp; ßfopen vs open fp points to a structure with information about the file

sigma=rmin/pow(2.0,1.0/6.0);

fp=fopen("blic.out","w"); ßthe File fp contains information about the file (fopen has open inside it)

if(fp ==NULL){printf(" cannot open bli.out file\n");

scanf(" %d ",&itest);

return;}

printf("np np2 \n");

scanf("%d %d",&np,&np2);

if(np2<np)np2=np;

xi = (double *) calloc(np2, sizeof(double)); ßallocating arrays

fi = (double *) calloc(np2, sizeof(double));

err = (double *) calloc(np2, sizeof(double));

printf("enter ndo \n");

scanf("%d",&ndo);

if(ndo < 2){printf(" ndo must be at least 2 %d",ndo);

scanf(" %d ",&itest);

return;}

else if(ndo == 2){printf(" enter the beginning and ending x values \n");

scanf(" %lg %lg",&xi[0],&xi[1]);}

ibli=bli(xi,fi,err,np,ndo,alj);

if(np2>np){ibli=bli(xi,fi,err,np2,np,alj);}

for(i=0;i<ibli;i++)fprintf(fp,"%14.6lg %14.6lg\n",xi[i],fi[i]);

fclose(fp);

system("gplot blic.out"); ß look at the data

c1=getc(stdin);

return;}

bli.for

FUNCTION FERRC(I,XI,FI,NDO) ß the interpolation routine

IMPLICIT REAL*8 (A-H,O-Z)

DIMENSION XI(NDO),FI(NDO)

C INTERPOLATES TO X(I)

IF(I.EQ.1)THEN ß the spacing and endif statement make this programmable

IF(XI(3).EQ.XI(2))THEN

FERRC=ABS(FI(3)-FI(2)) ß might be better set to 0

ELSE

FP=(FI(3)-FI(2))/(XI(3)-XI(2))

FINT=FI(2)+(XI(1)-XI(2))*FP

FERRC=ABS(FI(1)-FINT)

ENDIF

ELSEIF(I.EQ.NDO)THEN

IF(XI(NDO-1).EQ.XI(NDO-2))THEN

FERRC=ABS(FI(NDO)-FI(NDO-1))

ELSE

FP=(FI(NDO-1)-FI(NDO-2))/(XI(NDO-1)-XI(NDO-2))

FINT=FI(NDO-1)+(XI(I)-XI(NDO-1))*FP

FERRC=ABS(FI(NDO)-FINT)

ENDIF

ELSE

IF(XI(I+1).EQ.XI(I-1))THEN

FERRC=ABS(FI(I+1)-FI(I))

ELSE

FP=(FI(I+1)-FI(I-1))/(XI(I+1)-XI(I-1))

FINT=FI(I-1)+(XI(I)-XI(I-1))*FP

FERRC=ABS(FI(I)-FINT)

ENDIF

ENDIF

RETURN

END

SUBROUTINE BLI(XI,FI,ERR,NP,NDO,FEXT) ß keeps the terms in order

IMPLICIT REAL*8 (A-H,O-Z)

DIMENSION ERR(NP),XI(NP),FI(NP)

C *** NDO.EQ.2 makehalf of the point evenly spaced

C *** FI values need to found externally for NDO.NE.2

IF(NDO.EQ.2)THEN ß if NDO != to 2, the FI’s must be defined in the calling routine

NDO=NP/2

H=(XI(2)-XI(1))/(NDO-1)

DO I=1,NDO

XI(I)=XI(1)+(I-1)*H

FI(I)=FEXT(XI(I))

ENDDO

ENDIF

C *** ERR(I) is the interpolation error at X(I)

DO I=1,NDO

ERR(I)=FERRC(I,XI,FI,NDO)

ENDDO

C *** finding the region with the largest error

50 CONTINUE

ERRM=(ERR(1)+ERR(2))*ABS(XI(2)-XI(1))

ILE=1

DO I=2,NDO-1

ERRT=(ERR(I)+ERR(I+1))*ABS(XI(I+1)-XI(I))

IF(ERRT.GT.ERRM)THEN

ERRM=ERRT

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 ß the loop steps down, values are moved up

XI(I+1)=XI(I)

FI(I+1)=FI(I)

ERR(I+1)=ERR(I)

ENDDO

NDO=NDO+1

IADD=ILE+1

XI(IADD)=(XI(ILE)+XI(ILE+1))/2

FI(IADD)=FEXT(XI(IADD))

C *** The error at IADD-1 depends on FI(IADD)

C *** the error at IADD+1 depends on FI(IADD)

IB=MAX(1,IADD-1)

IE=MIN(NDO,IADD+1)

DO I=IB,IE ß calculatesthe changed err values

ERR(I)=FERRC(I,XI,FI,NDO)

ENDDO

IF(NP.GT.NDO)GOTO 50

RETURN

END

bli.c

#include <stdio.h>

#include <math.h>

#include <stdlib.h>

double fabs( double x ); ß beware abs is an integer in C

double ferrc(int i,double xi[],double fi[],int ndo)

{double ftemp,fp,fint;

if(i==0){ ßthe beginning of a multiline if is {

the triple segment USED IN FORTRAN did not easily work

if(xi[2]==xi[1])ftemp=fabs(fi[2]-fi[1]);

else {fp=(fi[2]-fi[1])/(xi[2]-xi[1]);

fint=fi[1]+(xi[0]-xi[1])*fp;

ftemp=fabs(fi[0]-fint);}} ß }} ends the two if segments

if(i==ndo-1){

if(xi[ndo-2]==xi[ndo-3])ftemp=fabs(fi[ndo-2]-fi[ndo-3]); ßfabs is double, abs is integer

else {fp=(fi[ndo-3]-fi[ndo-2])/(xi[ndo-3]-xi[ndo-2]);

fint=fi[ndo-2]+(xi[i]-xi[ndo-2])*fp;

ftemp=fabs(fi[ndo-1]-fint);}}

if(i!=0 && i!= ndo-1){ftemp=0;

if(xi[i+1]==xi[i-1])ftemp=fabs(fi[i+1]-fi[i-1]);

else {fp=(fi[i+1]-fi[i-1])/(xi[i+1]-xi[i-1]);

fint=fi[i-1]+(xi[i]-xi[i-1])*fp;

ftemp=fabs(fi[i]-fint);}}

return ftemp;}

int bli(double xi[],double fi[],double err[],int np,

int ndo, double (*fext)(double)) ß *fext is the external function

{int i,ib,ie,ile,iadd;

// int itest;

double h,errt,errm;

if(ndo == 2) ß if not the values of fi[0 to ndo-1] need to be set in calling code

{ndo=np/2;

h=(xi[1]-xi[0])/(ndo-1);

for (i=0;i < ndo;++i)

{xi[i]=xi[0]+h*i;

fi[i]=(*fext)(xi[i]);}}

for (i=0;i < ndo;++i)err[i]=ferrc(i,xi,fi,ndo); ß much more compact than in FORTRAN

f_largest_err:

errm=(err[0]+err[1])*fabs(xi[1]-xi[0]);

ile=0;

for (i=1;i<ndo-1;++i)

{errt=(err[i]+err[i+1])*fabs(xi[i+1]-xi[i]);

if(errt>errm){

errm=errt;

ile=i;}}

xi[ndo]=xi[ndo-1];

fi[ndo]=fi[ndo-1];

// move points up to ß //is a watcomextension to ansi C

for (i=ndo-2;i>=ile+1;--i)

{xi[i+1]=xi[i];

fi[i+1]=fi[i];

err[i+1]=err[i];}

ndo=ndo+1;

iadd=ile+1;

xi[iadd]=(xi[ile]+xi[ile+1])/2;

fi[iadd]=(*fext)(xi[iadd]);

ib=iadd-1;

if(ib<0)ib=0;

ie=iadd+1;

if(ie>ndo-1)ie=ndo-1;

// udating err

for(i=ib;i<=ie;++i)err[i]=ferrc(i,xi,fi,ndo);

if(np > ndo)goto f_largest_err;

return ndo;}

ß making the execuatbale required re-compiling both files, options, make, make all files.

Figure 2 This is the screen brought up by clicking on the running figure.

Figure 3 Initial screen from gplot

Enter the command <alt>vf

Enter m for a menu

Figure 4 gplot showing the menu

Enter 9, then 1 P to get the points. The command PrtScr will copy this to the clip board. Esc will shrink it down to a window from which select all, followed by copy will also copy this to the clipboard.

Figure 5 Lennard Jones potential represented by 27 points.

The plots for bli.out and blic.out are of course the same.