Introduction
The group of
points shown to the left becomes a function through the process of
interpolation. In general the region in
which the interpolating function is most valid is that between points 12 and
13. The fact that the polynomial goes through points 12 and 13 give rise to the
hope that it is reasonably accurate between them. The fact that it also goes through points 11
and 14 implies that the first derivatives should be accurate at each side of
the (12,13) region. The set of points
will in general be spaced such that there are more points in the regions where
the function is rapidly changing.
Practical
Handbook of Mathematical
Functions, M. Abramowitz and I.
Stegun, editors
Figure 1 Top of page from Abramawitz and Stegun
Figure 2 Bottom of the same page
Linear interpolation on j0(x) in this table has an error of 4 x 10-4. Interpolation with 6 Lagrange polynomials, gives the full table accuracy.
Linear interpolation
A linear polynomial that goes through the beginning and ending points of the data is
Note that the first derivative in the region is approximately given by
(2)
This makes
(3)
This is mostly of value as a convenient way to remember equation (1)
Locate Fortran #Tlocate in C
Thus the interpolation problem begins with finding the J of equation (1).This starts the test routine for locate. The best way to find what a routine does is to test it.
Test program for locate
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (MVALS=129) ß This is the f77 way to have variable dimensions. Only MVALS needs to be changed.
DIMENSION XI(MVALS),FI(MVALS)
OPEN(1,FILE='LOREN.OUT')
N=1
5 READ(1,*,ERR=10,END=10)XI(N),FI(N)
N=N+1
IF(N.GT.MVALS)THEN
PRINT*,' N => MVALS IN READ
READ(*,*)ITEST
GOTO 10
ENDIF
GOTO 5
10 N=N-1
15 PRINT*,' ENTER A VALUE FOR X '
READ(*,*)X
CALL LOCATE(X,XI,N,J)
IF(J.EQ.0)THEN
PRINT*,' X IS LESS THAN ',XI(1)
ELSEIF(J.EQ.N)THEN
PRINT*,' x IS GREATER THAN last point ',XI(N)
ELSE
PRINT*,' J= ',J
PRINT*,XI(J),X,XI(J+1)
ENDIF
GOTO 15
END
C$INCLUDE LOCATE
The Fortran locate routine begins
The routine TLOCATE.ZIP (Tlocate.for, LOCATE.FOR, loren.out)
searches an ordered vector to find the largest member less than the argument x. The following is “essentially” that given in Press (Fortran 1986) p 90. The method is the old 20 questions game. Is the result larger than x(j) halving the size of j each time. 2N values are search with N such questions.
SUBROUTINE LOCATE(X,R,NMAX,J)
C RETURNS A VALUE J SUCH THAT R(J)<=X<=R(J+1)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION R(*)
JL=0
JU=NMAX+1
10 IF(JU-JL.LE.1)GOTO 20
JM=(JU+JL)/2
IF(X.GT.R(JM))THEN
JL=JM
ELSE
JU=JM
ENDIF
GOTO 10
20 J=JL
RETURN
END
run
*WRN* ST-09 line 27, this statement will never be executed due to the preceding
branch
ENTER A VALUE FOR X
37
J= 72
31.746000000000000 37.000000000000000 39.682500000000000
ENTER A VALUE FOR X
25
J= 71
23.809500000000000 25.000000000000000 31.746000000000000
ENTER A VALUE FOR X
.05
J= 63
2.842170000000000D-014 0.050000000000000 1.984130000000000
ENTER A VALUE FOR X
-2000
X IS LESS THAN -1000.000000000000000
ENTER A VALUE FOR X
+1000
J= 127
968.254000000000000 1000.000000000000000 1000.000000000000000
ENTER A VALUE FOR X
s
*ERR* IO-08 bad character in input field
TraceBack: Executing in MAIN PROGRAM, statement 16 in file tlocate.for
<hold>
Enter a character such as ‘s’ to end. This ending is not elegant,
but is simple.
One last comment, in the ide, the c$include makes the
locate code appear both as a part of tlocate and separately, but we want them both so that we can edit from the ide and so that … I always change the c$ to a c#. for\TLOCATE.ZIP.
Tlocate in C
The C code to do this same thing follows. It is more elegant than the Fortran (and somewhat related to Press C 1992 p 117)
The tlocate source in C
#include <stdio.h>
#include <io.h>
#define MAX 129 /* equivalent to a parameter statement */
#include "locate.c" /* the guilty line */
void main(void)
{
double x,xi[MAX],fi[MAX];
int n,i;
FILE *fp;
fp = fopen("loren.out","r");
if(fp == NULL){printf(" cannot open file loren.out \n");
return;}
n=1;
redbeg:
if(fscanf(fp, "%lg %lg", &xi[n],&fi[n])!= EOF)
{printf("The line read was: %lg %g\n", xi[n],fi[n]);
n = n+1;
goto redbeg;}
fclose(fp);
if(n > MAX )
{printf(" from tlocate or read n > MAX, enter to continue ");
scanf(" \n");}
n=n-1;
while (x != -10000)
{printf("enter a value for x or -10000 to quit \n"); /* a more contrived ending than the s above */
scanf("%lf",&x);
i=locate(x,xi,n);
if(i >= 0 && i != n-1)
{printf(" i = %d \n",i);
printf(" %f %f %f \n",xi[i],x,xi[i+1]);}
else if (i == -1 )printf(" %f %f \n",x,xi[i+1]);
else if (i == n-1 )printf(" %f %f \n",xi[i],x); }
return; }
The locate source in C
cpp\TLOCATE.ZIP (cpp\tlocate.wpj cpp\TLOCATE.C cpp\LOCATE.C)
int locate(double x,double r[],int n)
{ int jl = -1;
int ju = n;
int jm;
while (ju-jl>1)
{jm=(ju+jl) >> 1;
if(x > r[jm])jl=jm;
else ju=jm; }
return jl;}
Note that locate is a function rather than a subroutine and that the division by 2 is simply a right shift.