home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
The C Users' Group Library 1994 August
/
wc-cdrom-cusersgrouplibrary-1994-08.iso
/
listings
/
v_09_05
/
9n05085a
< prev
next >
Wrap
Text File
|
1991-03-21
|
12KB
|
269 lines
/**********************************************************
*
* 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.
* --------------------------------------------------------
*
* Function: lint1d ( x0,dx,n, f,y )
*
* Purpose: Creates a table containing n values,
* starting with y[0]=f(x0), y[1]=f(x0+dx),
* through y[n-1]=f(x0+(n-1)*dx).
* However, the value in y[i] won't
* necessarily exactly equal f(x0+i*dx)
* as indicated above. Instead, the table
* is constucted to minimize the mean square
* error due to sampling random points in
* f(x)'s domain.
*
* Arguments:
* x0 (I) Double containing the first point
* in the function's domain.
* dx (I) Double containing x-interval
* between points in the table.
* n (I) Int containing the number of
* points in the table.
* f (I) Address of function returning
* double whose values are to be
* represented in the table.
* y (O) Address of double returning the
* table of n values, as discussed.
*
* Returns: (int) 0 for a normal solution, or
* 1 for any error (e.g., a singular
* set of equations).
*
* Notes: o See the accompanying article for a complete
* description of lint1d().
* o Set #defined value of TESTDRIVE to 1/TRUE
* to compile a sample main() (see below).
*
* Source: LINT1D.C
*
* --------------------------------------------------------
* Revision History:
*
* 01/07/91 J.Forkosh Installation
*
*********************************************************/
/* --- 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+1] = 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("lint1d> 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=1.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 = lint1d(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 --- */
printf(" 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 would 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+1,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