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 >
Text File  |  1991-03-21  |  12KB  |  269 lines
  1.  
  2. /**********************************************************
  3.  *
  4.  * Copyright (c) 1991, John Forkosh.  All rights reserved.
  5.  * Released to the public domain only for use that is both
  6.  * (a) by an individual, and (b) not for profit.
  7.  * --------------------------------------------------------
  8.  *
  9.  * Function:    lint1d ( x0,dx,n, f,y )
  10.  *
  11.  * Purpose:     Creates a table containing n values,
  12.  *              starting with y[0]=f(x0), y[1]=f(x0+dx),
  13.  *              through y[n-1]=f(x0+(n-1)*dx).
  14.  *              However, the value in y[i] won't
  15.  *              necessarily exactly equal f(x0+i*dx)
  16.  *              as indicated above.  Instead, the table
  17.  *              is constucted to minimize the mean square
  18.  *              error due to sampling random points in
  19.  *              f(x)'s domain.
  20.  *
  21.  * Arguments:
  22.  *              x0 (I)  Double containing the first point
  23.  *                      in the function's domain.
  24.  *              dx (I)  Double containing x-interval
  25.  *                      between points in the table.
  26.  *              n (I)   Int containing the number of
  27.  *                      points in the table.
  28.  *              f (I)   Address of function returning
  29.  *                      double whose values are to be
  30.  *                      represented in the table.
  31.  *              y (O)   Address of double returning the
  32.  *                      table of n values, as discussed.
  33.  *
  34.  * Returns:     (int)   0 for a normal solution, or
  35.  *                      1 for any error (e.g., a singular
  36.  *                      set of equations).
  37.  *
  38.  * Notes:     o See the accompanying article for a complete
  39.  *              description of lint1d().
  40.  *            o Set #defined value of TESTDRIVE to 1/TRUE
  41.  *              to compile a sample main() (see below).
  42.  *
  43.  * Source:      LINT1D.C
  44.  *
  45.  * --------------------------------------------------------
  46.  * Revision History:
  47.  *
  48.  * 01/07/91     J.Forkosh       Installation
  49.  *
  50.  *********************************************************/
  51.  
  52. /* --- standard headers --- */
  53. #include <stdio.h>              /* need NULL ptr value */
  54. #include <stdlib.h>             /*for malloc() prototype*/
  55. #define msglevel 1              /* to printf debug info */
  56.  
  57. /* --- set testdrive to compile (or not) test main() --- */
  58. #define TESTDRIVE 1             /* 1=compile it (0=not) */
  59.  
  60. /* --------------------------------------------------------
  61. Entry point
  62. -------------------------------------------------------- */
  63. int     lint1d ( x0,dx,n, f,y ) /* returns 0=ok, 1=error */
  64. double  x0,dx;                  /*1st point and interval*/
  65. int     n;                      /* number of points */
  66. double  (*f)();                 /*func to be interpolated*/
  67. double  *y;                     /*table for interpolation*/
  68. {
  69. /* --- local allocations and declarations --- */
  70. int     iserror = 0;            /* no error detected yet */
  71. int     i,j;                    /* row,col indexes */
  72. double  x;                      /* arg for f(x) */
  73. /* --- need temporary array for coefficient matrix --- */
  74. double  *a = (double *)malloc((n*n+1)*sizeof(double));
  75.  
  76. /* --------------------------------------------------------
  77. Set up coefficient matrix as per discussion in article
  78. (since the matrix is symmetric, the columnwise requirement
  79. for simq() input is irrelevant).
  80. -------------------------------------------------------- */
  81. /* --- first check that malloc() was successful --- */
  82. if ( a == NULL )                /*failed to malloc matrix*/
  83.   return ( iserror=1 );         /*return error to caller*/
  84. for (i=0;i<n*n;i++) a[i]=0.0;   /* init array to zeros */
  85. /* --- 4's on diag (2's at extremes) and 1's offdiag --- */
  86. for ( i=1; i<n; i++ )           /* skip 0,0 element */
  87.   {
  88.   j = n*i + i;                  /* index of diag element */
  89.   a[j] = 4.0;                   /* set diagonal element */
  90.   a[j-1] = a[j+1] = 1.0;        /* and off-diag elements */
  91.   } /* --- end-of-for(i=1...n-1) --- */
  92. a[0] = a[n*n-1] = 2.0;          /* set 1st,last diag */
  93. a[1] = a[n*n-2] = 1.0;          /*and remaining off-diags*/
  94.  
  95. /* --------------------------------------------------------
  96. Set up vector of constants as per discussion in article
  97. -------------------------------------------------------- */
  98. /* --- interior points --- */
  99. for ( i=1; i<n; i++ )           /* loop skips 1st point */
  100.   {
  101.   x = x0 + dx*i;                /* next arg to be tabled */
  102.   y[i] = 2.0*(f(x-0.5*dx)+      /*const vector calculated*/
  103.                 f(x)+f(x+0.5*dx)); /* as per article */
  104.   } /* --- end-of-for(i=1...n-1) --- */
  105. /* --- boundary points --- */
  106. y[0]   = f(x0) + 2.0*f(x0+0.5*dx); /* first x in domain */
  107. y[n-1] = f(x)  + 2.0*f(x -0.5*dx); /*use last x from loop*/
  108.  
  109. /* --------------------------------------------------------
  110. Display input to simq (for debugging purposes if necessary)
  111. -------------------------------------------------------- */
  112. if ( msglevel >= 9              /* want debugging output */
  113. &&   n < 10 )                   /* have room to fit it */
  114.   {
  115.   printf("lint1d> input to simq...\n"); /* stub info */
  116.   for ( i=0; i<n; i++ )
  117.     { /* --- display row (really col) of the matrix --- */
  118.       for ( j=0; j<n; j++ ) printf("%6.2lf",a[n*i+j]);
  119.       /* --- followed by literal for y* and const y --- */
  120.       printf("    y*[%d]    %8.2lf\n", i,y[i]); }
  121.   } /* --- end-of-if(msglevel>=9) --- */
  122.  
  123. /* --------------------------------------------------------
  124. Solve the simultaneous equations
  125. -------------------------------------------------------- */
  126. iserror = simq(a,y,n);          /* y[] returns solution */
  127. free ( (void *)a );             /*don't need coeff matrix*/
  128. return ( iserror );             /*back to caller with y[]*/
  129. } /* --- end-of-function lint1d() --- */
  130.  
  131.  
  132. #if TESTDRIVE
  133. /**********************************************************
  134.  *
  135.  * Purpose:     Test driver for lint1d().
  136.  *              Prepares a table of values for
  137.  *              interpolation of the function f(x)=x*x,
  138.  *              and compares the resulting accuracy
  139.  *              with usual linear interpolation.
  140.  *
  141.  * Command-Line Arguments:
  142.  *              x0 (I)  Double containing the first point
  143.  *                      in the function's domain
  144.  *                      (default=-10.0).
  145.  *              dx (I)  Double containing x-interval
  146.  *                      between points in the table
  147.  *                      (default=1.0).
  148.  *              n (I)   Int containing the number of
  149.  *                      points in the table
  150.  *                      (default=21).
  151.  *              expon (I) Int containing the exponent for
  152.  *                      the test function f(x)=x**expon
  153.  *                      (default=2).
  154.  *
  155.  *********************************************************/
  156.  
  157. /* --- program defaults (x0,dx,n,expon from cmdline) --- */
  158. #define X0 (-10.0)              /* 1st point in table */
  159. #define DX 1.0                  /* table interval */
  160. #define N 21                    /*number pts in table*/
  161. static  int expon=2;            /* test f(x)=x**expon */
  162. #define NMAX 99                 /* easier than malloc */
  163. #define NRMS 101                /*pts per dx for rms calc*/
  164.  
  165. /* --------------------------------------------------------
  166. Entry point
  167. -------------------------------------------------------- */
  168. main    ( argc, argv )
  169. int     argc;
  170. char    *argv[];
  171. {
  172. /* --- local allocations and declarations --- */
  173. double  x0=X0, dx=DX; int n=N;  /* defaults */
  174. int     i,j;                    /* loop indexes */
  175. double  x;                      /* arg for f() */
  176. double  f();                    /* test function */
  177. double  y[NMAX];                /* table from lint1d() */
  178. int     iserror=0;              /*return flag from lint1d*/
  179. double  frms=0.0,yrms=0.0;      /* mean square errors */
  180. double  xa,xb;                  /*interval bounds for rms*/
  181. double  sqrarg;                 /*dummy arg for sqr macro*/
  182. /* --- linear interpolation (u(x=xa)=ya, u(x=xb)=yb) --- */
  183. #define u(x,xa,ya,xb,yb) (ya*(xb-x)+yb*(x-xa))/(xb-xa)
  184. /* --- square a double --- */
  185. #define sqr(x) (sqrarg=(x),sqrarg*sqrarg)
  186.  
  187. /* --------------------------------------------------------
  188. Pick up command-line arguments
  189. -------------------------------------------------------- */
  190. if ( *(argv[1]) == '?' )        /*request for usage info*/
  191.   { printf("Usage: lint1d x0 dx n expon\n");
  192.     exit(0); }                  /*only a little help here*/
  193. if ( argc > 1 )                 /* user supplied x0 */
  194.   x0 = atof(argv[1]);           /* so replace default x0 */
  195. if ( argc > 2 )                 /* user supplied dx */
  196.   dx = atof(argv[2]);           /* so replace default dx */
  197. if ( argc > 3 )                 /* user supplied n */
  198.   n  = atoi(argv[3]);           /* so replace default n */
  199. if ( argc > 4 )                 /* user supplied expon */
  200.   expon  = atoi(argv[4]);       /* replace default expon */
  201. if ( n<3 || n>NMAX ) exit(0);   /* sorry, Charlie */
  202.  
  203. /* --------------------------------------------------------
  204. Call lint1d() to create interpolation table.  Note: After
  205. this call, a normal application (what mainframe IBM
  206. likes to call a "problem program") would simply go about
  207. its business, using y[i=0...n-1] as a table for linear
  208. interpolation of f(x) in the usual way.  The special
  209. algorithm by which our table was generated is completely
  210. transparent to further processing.
  211. -------------------------------------------------------- */
  212. iserror = lint1d(x0,dx,n,f,y);  /*very easy to use lint1d*/
  213. if ( iserror )                  /* unless it fails */
  214.   { printf("lint1d() failed.\n"); /* then just print msg */
  215.     exit(0); }                  /* and quit */
  216.  
  217. /* --------------------------------------------------------
  218. Print the exact and tabled values of f(x) at each x point
  219. -------------------------------------------------------- */
  220. /* --- stub header --- */
  221. printf("        Table Returned By Lint1d()\n");
  222. printf(" i     x[i]       f(x[i])        y[i]\n");
  223. printf("---  --------  ------------  ------------\n");
  224. /* --- table generated by lint1d() in right column --- */
  225. for ( i=0; i<n; i++ )           /*for each value in table*/
  226.   { /* --- f(x[i]) is what a normal table would contain
  227.               y[i] is our table to minimize errors --- */
  228.   x = x0 + dx*i;                /* need x to print f(x) */
  229.   printf("%3d  %8.2lf  %12.3lf  %12.6lf\n",
  230.     i+1,x,f(x),y[i]);           /* f(x) and y[i] */
  231.   } /* --- end-of-for(i=0...n-1) --- */
  232.  
  233. /* --------------------------------------------------------
  234. Determine mean square error for both f(x[i]) and for y[i]
  235. -------------------------------------------------------- */
  236. for ( i=1; i<n; i++ )           /*each interval in table*/
  237.   {
  238.   double fxa, fxb, fx;          /* f(xa), f(xb), f(x) */
  239.   xa = x0+dx*(i-1); fxa=f(xa);  /*lower bound of interval*/
  240.   xb = x0+dx*i;     fxb=f(xb);  /* upper bound */
  241.   for ( j=0; j<NRMS; j++ )      /*accum err for NRMS pts*/
  242.     {
  243.     x = xa + (j*dx)/(NRMS-1);   /* xa<=x<=xb */
  244.     fx = f(x);                  /* actual value at x */
  245.     frms += sqr(u(x,xa,fxa,xb,fxb)-fx); /* usual error */
  246.     yrms += sqr(u(x,xa,y[i-1],xb,y[i])-fx); /* our error */
  247.     } /* --- end-of-for(j=0...NRMS-1) --- */
  248.   } /* --- end-of-for(i=1...n-1) --- */
  249. /* --- print usual error and our (smaller) error --- */
  250. printf("Mean sq err for\n%d pts/intvl  %12.8lf  %12.8lf\n",
  251.   NRMS,frms/(NRMS*(N-1)),yrms/(NRMS*(N-1)));
  252. } /* --- end-of-main() --- */
  253.  
  254. /* --------------------------------------------------------
  255. Entry point (dummy function f(x)=x**expon can be replaced
  256. by any function of the user's choice)
  257. -------------------------------------------------------- */
  258. double  f(x)
  259. double  x;
  260. {
  261. /* --- replace this with any function of your choice --- */
  262. int     i=expon;                /* countdown index */
  263. double  y=x;                    /* start with x */
  264. while ( --i > 0 ) y *= x;       /*additional factors of x*/
  265. return ( y );                   /* f(x)=x**expon */
  266. } /* --- end-of-function f() --- */
  267. #endif
  268.  
  269.