/**********************************************************
* Copyright (c) 1991, John Forkosh. All rights reserved.
* Released to the public domain only for use that is both
* (a) by an individual, and (b) not for profit.
*********************************************************/
/* --- standard headers --- */
#include <stdio.h> /* need NULL ptr value */
#include <stdlib.h> /*for malloc() prototype*/
#define msglevel 1 /* to printf debug info */
/* --- set testdrive to compile (or not) test main() --- */
#define TESTDRIVE 1 /* 1=compile it (0=not) */
/* --------------------------------------------------------
Entry point
---------------------------------------------------------*/
int lint1d ( x0,dx,n, f,y ) /* returns 0=ok, 1=error */
double x0,dx; /*1st point and interval*/
int n; /* number of points */
double (*f)(); /*func to be interpolated*/
double *y; /*table for interpolation*/
{
/* --- local allocations and declarations --- */
int iserror = 0; /* no error detected yet */
int i,j; /* row,col indexes */
double x; /* arg for f(x) */
/* --- need temporary array for coefficient matrix --- */
double *a = (double *)malloc((n*n+1)*sizeof(double));
/* --------------------------------------------------------
Set up coefficient matrix as per discussion in article
(since the matrix is symmetric, the columnwise requirement
for simq() input is irrelevant).
---------------------------------------------------------*/
/* --- first check that malloc() was successful --- */
if ( a == NULL ) /*failed to malloc matrix*/
return ( iserror=1 ); /*return error to caller*/
for (i=0;i<n*n;i++) a[i]=0.0; /* init array to zeros */
/* --- 4's on diag (2's at extremes) and 1's offdiag --- */
for ( i=1; i<n; i++ ) /* skip 0,0 element */
{
j = n*i + i; /* index of diag element */
a[j] = 4.0; /* set diagonal element */
a[j-1] = a[j+l] = 1.0; /* and off-diag elements */
} /* --- end-of-for(i=1...n-1) --- */
a[0] = a[n*n-1) = 2.0; /* set 1st,last diag */
a[1]= a[n*n-2] = 1.0; /*and remaining off-diags*/
/*---------------------------------------------------------
Set up vector of constants as per discussion in article
---------------------------------------------------------*/
/* --- interior points --- */
for ( i=1; i<n; i++ ) /* loop skips 1st point */
{
x = x0 + dx*i; /* next arg to be tabled */
y[i] = 2.0*(f(x-0.5*dx)+ /*const vector calculated*/
f(x)+f(x+0.5*dx)); /* as per article */
} /* --- end-of-for(i=1...n-1) --- */
/* --- boundary points --- */
y[0] = f(x0) + 2.0*f(x0+0.5*dx); /* first x in domain */
y[n-1] = f(x) + 2.0*f(x -0.5*dx); /*use last x from loop*/
/* --------------------------------------------------------
Display input to simq (for debugging purposes if necessary)
---------------------------------------------------------*/
if ( msglevel >= 9 /* want debugging output */
&& n < 10 ) /* have room to fit it */
{
printf("lintld> input to simq...\n"); /* stub info */
for ( i=0; i<n; i++ )
{ /* --- display row (really col) of the matrix --- */
for ( j=0; j<n; j++ ) printf("%6.2lf",a[n*i+j]);
/* ------ followed by literal for y* and const y ------ */
printf(" y*[%d] %8.2lf\n", i,y[i]); }
} /* ------ end-of-if(msglevel>=9) --- */
/* --------------------------------------------------------
Solve the simultaneous equations
---------------------------------------------------------*/
iserror = simq(a,y,n); /* y[] returns solution */
free ( (void *)a ); /*don't need coeff matrix*/
return ( iserror ); /*back to caller with y[]*/
} /* --- end-of-function lint1d() --- */
#if TESTDRIVE
/***********************************************************
*
* Purpose: Test driver for lint1d().
* Prepares a table of values for
* interpolation of the function f(x)=x*x,
* and compares the resulting accuracy
* with usual linear interpolation.
*
* Command--Line Arguments:
* x0 (I) Double containing the first point
* in the function's domain
*
* (default=-10.0).
* dx (I) Double containing x-interval
* between points in the table
* (default=l.0).
* n (I) Int containing the number of
* points in the table
* (default=21).
* expon (I) Int containing the exponent for
* the test function f(x)=x**expon
* (default=2).
*
**********************************************************/
/* ------ program defaults (x0,dx,n,expon from cmdline) --- */
#define X0 (-10.0) /* 1st point in table */
#define DX 1.0 /* table interval */
#define N 21 /*number pts in table*/
static int expon=2; /* test f(x)=x**expon */
#define NMAX 99 /* easier than malloc */
#define NRMS 101 /*pts per dx for rms calc*/
/* --------------------------------------------------------
Entry point
-------------------------------------------------------- */
main ( argc, argv )
int argc;
char *argv[];
{
/* ------ local allocations and declarations ------ */
double x0=X0, dx=DX; int n=N; /* defaults */
int i,j; /* loop indexes */
double x; /* arg for f() */
double f(); /* test function */
double y[NMAX]; /* table from lint1d() */
int iserror=0; /*return flag from lint1d*/
double frms=0.0,yrms=0.0; /* mean square errors */
double xa,xb; /*interval bounds for rms*/
double sqrarg; /*dummy arg for sqr macro*/
/* ------ Linear interpolation (u(x=xa)=ya, u(x=xb)=yb) ------ */
#define u(x,xa,ya,xb,yb) (ya*(xb-x)+yb*(x-xa))/(xb-xa)
/* ------ square a double ------ */
#define sqr(x) (sqrarg=(x),sqrarg*sqrarg)
/* --------------------------------------------------------
Pick up command-line arguments
-------------------------------------------------------- */
if ( *(argv[1]) == '?' ) /*request for usage info*/
{ printf("Usage: lint1d x0 dx n expon\n");
exit(0); } /*only a little help here*/
if ( argc > 1 ) /* user supplied x0 */
x0 = atof(argv[1]); /* so replace default x0 */
if ( argc > 2 ) /* user supplied dx */
dx = atof(argv[2]); /* so replace default dx */
if ( argc > 3 ) /* user supplied n */
n = atoi(argv[3]); /* so replace default n */
if ( argc > 4 ) /* user supplied expon */
expon = atoi(argv[4]); /* replace default expon */
if ( n<3 || n>NMAX ) exit(0); /* sorry, Charlie */
/* --------------------------------------------------------
Call lint1d() to create interpolation table. Note: After
this call, a normal application (what mainframe IBM
likes to call a "problem program") would simply go about
its business, using y [i=0...n-1] as a table for linear
interpolation of f(x) in the usual way. The special
algorithm by which our table was generated is completely
transparent to further processing.
-------------------------------------------------------- */
iserror = lintld(x0,dx,n,f,y); /*very easy to use lint1d*/
if ( iserror ) /* unless it fails */
{ printf("lint1d() failed.\n"); /* then just print msg */
exit(0); } /* and quit */
/* --------------------------------------------------------
Print the exact and tabled values of f(x) at each x point
-------------------------------------------------------- */
/* ------ stub header ------ */
print(" Table Returned By Lint1d()\n");
printf(" i x[i] f(x[i]) y[i]\n");
printf(" ------ ---------------- ------------------------ ------------------------\n");
/* ------ table generated by lint1d() in right column ------ */
for ( i=0; i<n; i++ /*for each value in table*/
{ /* ------ f(x[i]) is what a normal table woud contain
y[i] is our table to minimize errors ------ */
x = x0 + dx*i; /* need x to print f(x) */
printf("%3d %.8.2lf %12.3lf %12.6lf\n",
i+l,x,f(x),y[i]); /* f(x) and y[i] */
} /* ------ end-of-for(i=0...n-1) ------ */
/* --------------------------------------------------------
Determine mean square error for both f(x[i]) and for y[i]
-------------------------------------------------------- */
for ( i=1; i<n; i++ ) /*each interval in table*/
{
double fxa, fxb, fx; /* f(xa), f(xb), f(x) */
xa = x0+dx*(i-1); fxa=f(xa); /*lower bound of interval*/
xb = x0+dx*i; fxb=f(xb); /* upper bound */
for ( j=0; j<NRMS; j++ ) /*accum err for NRMS pts*/
{
X = xa + (j*dx)/(NRMS-1); /* xa<=x<=xb */
fx = f(x); /* actual value at x */
frms += sqr(u(x,xa,fxa,xb,fxb)-fx); /* usual error */
yrms += sqr(u(x,xa,y[i-1],xb,y[i])-fx); /* our error */
} /* --- end-of-for(j=0...NRMS-1) --- */
} /* --- end-of-for(i=1...n-1) --- */
/* --- print usual error and our (smaller) error --- */
printf("Mean sq err for\n%d pts/intvl %12.8lf %12.8lf\n",
NRMS,frms/(NRMS*(N-1)),yrms/(NRMS*(N-1)));
} /* --- end-of-main() --- */
/* --------------------------------------------------------
Entry point (dummy function f(x)=x**expon can be replaced
by any function of the user's choice)
-------------------------------------------------------- */
double f(x)
double x;
{
/* --- replace this with any function of your choice --- */
int i=expon; /* countdown index */
double y=x; /* start with x */
while ( --i > 0 ) y *= x; /*additional factors of x*/
return ( y ); /* f(x)=x**expon */
} /* --- end-of-function f() --- */
#endif