home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / LinAlg 3.1 / LinAlg / vali.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1995-12-20  |  5.0 KB  |  204 lines  |  [TEXT/ttxt]
  1. // This may look like C code, but it is really -*- C++ -*-
  2. //
  3. //        Verification of the Aitken-Lagrange Interpolation
  4. //
  5.  
  6. #include "LinAlg.h"
  7. #include "math_num.h"
  8. #include <iostream.h>
  9.  
  10. /*
  11.  *------------------------------------------------------------------------
  12.  *        Set of the test functions to interpolate
  13.  */
  14.  
  15. static double (*F)(const double x);
  16.  
  17. static double f_exp(const double x)
  18. {
  19.   return exp(-x);
  20. }
  21.  
  22. static double f_sinexp(const double x)
  23. {
  24.   return sin(x) * exp(-x/10);
  25. }
  26.  
  27. /*
  28.  *------------------------------------------------------------------------
  29.  *    Set of tests for the interpolation over the uniform grid
  30.  */
  31.  
  32.             // Verify that interpolation exactly at a node
  33.             // really gives the function value at the node
  34. static void test_u_exact(const double s, const int ia, const int ib)
  35. {
  36.   cout << "\nCheck to see that interpolation at a node gives an exact"
  37.           " result\n";
  38.   cout << "\nUniform grid [" << s*ia << ':' << s*ib << "] with the step " << s << endl;
  39.  
  40.   Vector y(ia,ib);
  41.   Vector y_int(ia,ib);
  42.  
  43.   register int i;
  44.   for(i=y.q_lwb(); i<=y.q_upb(); i++)
  45.     y(i) = (*F)(s*i);
  46.  
  47.   for(i=y_int.q_lwb(); i<=y_int.q_upb(); i++)
  48.     y_int(i) = ali(s*i,s*ia,s,y);
  49.  
  50.   assert(y == y_int);
  51.  
  52.   cout << "\nDone\n";
  53. }
  54.  
  55.             // Check the precision for interpolation
  56.             // at arbitrary points
  57. static void test_u(const double s, const int ia, const int ib)
  58. {
  59.   cout << "\nCheck the precision of the interpolation at arbitrary points\n";
  60.   cout << "\nUniform grid [" << s*ia << ':' << s*ib << "] with the step " << s
  61.        << endl;
  62.  
  63.   Vector y(ia,ib);
  64.   register int i;
  65.   for(i=y.q_lwb(); i<=y.q_upb(); i++)
  66.     y(i) = (*F)(s*i);
  67.  
  68.   cout << "\nPoint    Exact function value     Interpolated    Error\n";
  69.   register double q;
  70.   for(q=(ia-1)*s; q < (ib+2)*s; q += 1.1*s)
  71.   {
  72.     double yexact = (*F)(q);
  73.     double yint = ali(q,ia*s,s,y);
  74.     printf("%4.2g\t%12.6g\t\t%12.6g  %12.6g\n",q,yexact,yint,
  75.          abs(yexact-yint));
  76.   }
  77.  
  78.   cout << "\nDone\n";
  79. }
  80.  
  81. /*
  82.  *------------------------------------------------------------------------
  83.  *    Set of tests for the interpolation over the non-uniform grid
  84.  */
  85.  
  86.             // Verify that interpolation exactly at a node
  87.             // really gives the function value at the node
  88. static void test_n_exact(const Vector& x)
  89. {
  90.   cout << "\nCheck to see that interpolation at a node gives an exact"
  91.           " result\n";
  92.  
  93.   Vector y(x);
  94.   Vector y_int(x);
  95.  
  96.   register int i;
  97.   for(i=y.q_lwb(); i<=y.q_upb(); i++)
  98.     y(i) = (*F)(x(i));
  99.  
  100.   for(i=y_int.q_lwb(); i<=y_int.q_upb(); i++)
  101.     y_int(i) = ali(x(i),x,y);
  102.  
  103.   assert(y == y_int);
  104.  
  105.   cout << "\nDone\n";
  106. }
  107.  
  108.             // Check the precision for interpolation
  109.             // at arbitrary points
  110. static void test_n(const Vector& x,const double s, const int ia, const int ib)
  111. {
  112.   cout << "\nCheck the precision of the interpolation at arbitrary points\n";
  113.  
  114.   Vector y(x);
  115.   register int i;
  116.   for(i=y.q_lwb(); i<=y.q_upb(); i++)
  117.     y(i) = (*F)(x(i));
  118.  
  119.   cout << "\nPoint    Exact function value     Interpolated    Error\n";
  120.   register double q;
  121.   for(q=ia*s; q < ib*s; q += s)
  122.   {
  123.     double yexact = (*F)(q);
  124.     double yint = ali(q,x,y);
  125.     printf("%4.2g\t%12.6g\t\t%12.6g  %12.6g\n",q,yexact,yint,
  126.          abs(yexact-yint));
  127.   }
  128.  
  129.   cout << "\nDone\n";
  130. }
  131.  
  132.             // Check the precision for interpolation
  133.             // on example in the book
  134. static void test_n_book(void)
  135. {
  136.   cout << "\nCheck the precision of the interpolation at arbitrary points\n";
  137.   cout << "\nExample from Fig. 4.11 of the book"
  138.           "\n\tNumerical Methods and Software,"
  139.       "\n\tby D.Kahaner, C.Moler, and S.Nash - Prentice Hall, 1989\n";
  140.  
  141.   REAL xy[2][11] =
  142.   { {0,  2,  3,  5,  6,  8,  9,    11, 12, 14, 15},    // abscissae
  143.     {10, 10, 10, 10, 10, 10, 10.5, 15, 50, 60, 85}    // ordinates
  144.   };
  145.  
  146.   Vector x(0,sizeof(xy[0])/sizeof(xy[0][0])-1);
  147.   Vector y(x);
  148.   register int i;
  149.  
  150.   for(i=0; i<sizeof(xy[0])/sizeof(xy[0][0]); i++)
  151.     x(i) = xy[0][i], y(i) = xy[1][i];
  152.  
  153.   cout << "\n\t\tInterpolation nots";
  154.   cout << "\nx ";
  155.   for(i=x.q_lwb(); i<=x.q_upb(); i++)
  156.     printf("%6.2g ",x(i));
  157.   cout << "\ny ";
  158.   for(i=y.q_lwb(); i<=y.q_upb(); i++)
  159.     printf("%6.2g ",y(i));
  160.  
  161.   cout << "\n\nPoint    Interpolated value\n";
  162.   register double q;
  163.   for(q=0; q <= 16; q += 16/50.)
  164.     printf("%7.4g   %12.6g\n",q,ali(q,x,y));
  165.  
  166.   cout << "\nDone\n";
  167. }
  168.  
  169. /*
  170.  *------------------------------------------------------------------------
  171.  *                Root module
  172.  */
  173.  
  174. main()
  175. {
  176.   cout << "\n\n\t\tVerification of the Aitken-Lagrange interpolation\n";
  177.  
  178.   cout << "\n------------- Interpolation over the table with uniform grid";
  179.  
  180.   cout << "\nFunction to interpolate: exp(-x)\n";
  181.   F = f_exp;
  182.  
  183.   const double s = 0.1;
  184.   test_u_exact(s,1,10);
  185.   test_u(s,1,10);
  186.  
  187.   printf("\n------------- Interpolation over the table with non-uniform grid");
  188.  
  189.   printf("\nFunction to interpolate: sin(x) * exp(-x/10)\n");
  190.   F = f_sinexp;
  191.  
  192.   register int i;
  193.   Vector nodes(2,51);
  194.   for(i=nodes.q_lwb(); i<=11; i++)
  195.     nodes(i) = (i-1.)/10;
  196.   for(i=12; i<=nodes.q_upb(); i++)
  197.     nodes(i) = (i-6.)/5;
  198.   cout << "\nGrids 0.1 .. 1 have the mesh 0.1, and 1..9 have the mesh 0.2\n";
  199.  
  200.   test_n_exact(nodes);
  201.   test_n(nodes,0.27,0,15);
  202.   test_n_book();
  203. }
  204.