home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsv / ch3_4 / xcoord.c
Encoding:
C/C++ Source or Header  |  1995-04-04  |  8.0 KB  |  273 lines
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3.  
  4. #define SQR(a) ((a)*(a))
  5.  
  6. typedef double MATX[10][10];
  7. typedef double VECT[10];
  8. typedef struct {double x; double y;} Point2;
  9.  
  10. Point2 pt[1023];                 /* From coordinates */
  11. double zeta[1023], eta[1023];    /* To   coordinates */
  12.  
  13. int Gauss(MATX ain, VECT bin, int n, VECT v)
  14. /*  Gaussian elimination by converting to upper triangular system.
  15.     Row interchanges are done via re-indexing sub[]. See Vandergraft:
  16.     Intro. Numerical Methods, 2ed, Academic Press, 1983, Chapter 6. */
  17. {   MATX a;  VECT b;
  18.     int  i, j, k, last, index;
  19.     double big, absv;
  20.     int  sub[21];
  21.  
  22.     for(k=0; k < n; k++) {       /* make local copies */
  23.         for(j=0; j < n; j++) a[k][j] = ain[k][j];
  24.         b[k] = bin[k];
  25.         }
  26.  
  27.     last= n-1;
  28.     for (k= 0; k <= last; k++) sub[k]= k;
  29.     for (k= 0; k <= last-1; k++) {
  30.         big= 0.0;
  31.         for (i= k; i <= last; i++) {
  32.             absv= abs(a[sub[i]][k]);
  33.             if (absv > big)
  34.                { big= absv; index= i; }
  35.             }
  36.         if (big == 0.0) return 0;
  37.         j= sub[k];
  38.         sub[k]= sub[index];
  39.         sub[index]= j;
  40.         big= 1.0/a[sub[k]][k];
  41.         for (i= k+1; i <= last; i++) {
  42.             a[sub[i]][k]= -a[sub[i]][k]*big;
  43.             for (j= k+1; j <= last; j++)
  44.                  a[sub[i]][j] += a[sub[i]][k] * a[sub[k]][j];
  45.             b[sub[i]]    += a[sub[i]][k] * b[sub[k]];
  46.             }
  47.         }
  48.         v[last]= b[sub[last]] / a[sub[last]][last];
  49.         for (k= last-1; k >= 0; k--) {
  50.             v[k]= b[sub[k]];
  51.             for (i= k+1; i <= last; i++)
  52.                  v[k] = v[k] -a[sub[k]][i] * v[i];
  53.             v[k] = v[k] /a[sub[k]][k];
  54.             }
  55.         return 1;
  56. }
  57.  
  58. void PrintMatrix(MATX a, VECT v, int size)
  59. {   int r, c;
  60.     for(r= 0; r < size; r++) {
  61.         for(c= 0; c < size; c++) printf("%14.6lf  ",a[r][c]);
  62.         printf(" %14.6lf\n", v[r]);
  63.         }
  64.      printf("\n");
  65. }
  66.  
  67. void PrintSolution(VECT v, int vectorsize, char which)
  68. /*  Print the solution vector */
  69. {   int k;
  70.     printf("Solution vector %c\n", which);
  71.     for(k = 0; k < vectorsize; k++)
  72.         if (abs(v[k]) < 1.0E6) printf("%14.6f  ",v[k]);
  73.         else printf("%14.6e  ", v[k]);
  74.     printf("\n");
  75. }
  76.  
  77. void FirstOrderExact(VECT xv, VECT yv)
  78. {   int k, ok;  VECT b;  MATX c;
  79.     for(k= 0; k<=2; k++) b[k] = zeta[k];
  80.     for(k= 0; k<=2; k++) {
  81.         c[k][0] = pt[k].x;
  82.         c[k][1] = pt[k].y;
  83.         c[k][2] = 1.0;
  84.         };
  85.  
  86.      printf("Augmented matrix:\n");
  87.      PrintMatrix(c, b, 3);
  88.      ok =Gauss(c, b, 3, xv);
  89.      PrintSolution(xv, 3, 'X');
  90.  
  91.      for(k= 0; k<=2; k++) b[k] = eta[k];
  92.      for(k= 0; k<=2; k++) {
  93.          c[k][0] = pt[k].x;
  94.          c[k][1] = pt[k].y;
  95.          c[k][2] = 1.0;
  96.          };
  97.  
  98.      PrintMatrix(c, b, 3);
  99.      ok =Gauss(c, b, 3, yv);
  100.      PrintSolution(yv, 3, 'Y');
  101. }
  102.  
  103. void SecondOrderExact(VECT xv, VECT yv)
  104. {   int k, ok;  VECT b;  MATX c;
  105.     for(k= 0; k<=5; k++) b[k] = zeta[k];
  106.     for(k= 0; k<=5; k++) {
  107.         c[k][0] = pt[k].x*pt[k].x;
  108.         c[k][1] = pt[k].y*pt[k].y;
  109.         c[k][2] = pt[k].x;
  110.         c[k][3] = pt[k].y;
  111.         c[k][4] = pt[k].x*pt[k].y;
  112.         c[k][5] = 1;
  113.         }
  114.  
  115.      printf("Augmented matrix:\n");
  116.      PrintMatrix(c, b, 6);
  117.      ok =Gauss(c, b, 6, xv);
  118.      printf("x = a5*x^2 + a4*y^2 + a3*x + a2*y + a1*x*y + a0:\n");
  119.      PrintSolution(xv, 6, 'X');
  120.  
  121.      for(k= 0; k<=5; k++) b[k] = eta[k];
  122.      for(k= 0; k<=5; k++) {
  123.          c[k][0] = SQR(pt[k].x);
  124.          c[k][1] = SQR(pt[k].y);
  125.          c[k][2] = pt[k].x;
  126.          c[k][3] = pt[k].y;
  127.          c[k][4] = pt[k].x*pt[k].y;
  128.          c[k][5] = 1;
  129.          }
  130.  
  131.      printf("Augmented matrix:\n");
  132.      PrintMatrix(c, b, 6);
  133.      ok =Gauss(c, b, 6, yv);
  134.      printf("y = b5*x^2 + b4*y^2 + b3*x + b2*y + b1*x*y + b0:\n");
  135.      PrintSolution(yv, 6, 'Y');
  136. }
  137.  
  138. void FirstOrderLeastSquares(int npoints, VECT xv, VECT yv)
  139. {   MATX c;  VECT b;
  140.     double sumx= 0, sumxx= 0, sumy= 0, sumyy= 0, sumxy= 0,
  141.            sumd= 0, sumdx= 0, sumdy = 0;
  142.     int k, ok;
  143.     double xt, yt;
  144.     for(k=0; k < npoints; k++) {
  145.         sumx  += pt[k].x;
  146.         sumxx += SQR(pt[k].x);
  147.         sumy  += pt[k].y;
  148.         sumyy += SQR(pt[k].y);
  149.         sumxy += pt[k].x*pt[k].y;
  150.         sumd  += zeta[k];
  151.         sumdx += pt[k].x*zeta[k];
  152.         sumxy += pt[k].y*zeta[k];
  153.         }
  154.  
  155.     c[0][0] = sumxx;   c[0][1] = sumxy;   c[0][2] = sumx;
  156.     c[1][0] = sumxy;   c[1][1] = sumyy;   c[1][2] = sumy;
  157.     c[2][0] = sumx;    c[2][1] = sumy;    c[2][2] = npoints;
  158.  
  159.     b[0] = sumdx;      b[1] = sumdy;      b[2] = sumd;
  160.     ok = Gauss(c, b, 3, xv);
  161.  
  162.     sumd = sumdx = sumdy = 0;
  163.  
  164.     for(k=0; k < npoints; k++) {
  165.         sumd += eta[k];
  166.         sumdx+= pt[k].x*eta[k];
  167.         sumxy+= pt[k].y*eta[k];
  168.         };
  169.  
  170.      b[0] = sumdx;  b[1] = sumdy; b[2] = sumd;
  171.      ok = Gauss(c, b, 3, yv);
  172.      printf("residuals\n");
  173.      for(k=0; k < npoints; k++) {
  174.         xt = zeta[k] -(pt[k].x*xv[0] + pt[k].y*xv[1] + xv[2]);
  175.         yt =  eta[k] -(pt[k].x*yv[0] + pt[k].y*yv[1] + yv[2]);
  176.         printf("%4d   %12.6   %12.6\n", xt, yt);
  177.         }
  178. }
  179.  
  180. void SecondOrderLeastSquares(MATX c, int npoints, VECT xv, VECT yv)
  181. {   int j, k, ok;
  182.     VECT b;
  183.     double sumd=0, sumdx=0, sumdx2=0, sumdy=0,
  184.                    sumdy2=0, sumdxy =0;
  185.     double px2, py2, xt, yt;
  186.  
  187.     for(j=0; j<= 5; j++)
  188.         for(k=0; k<= 5; k++) c[j][k] = 0;
  189.  
  190.     for(k =0; k < npoints; k++) {
  191.         px2 = SQR(pt[k].x);
  192.         py2 = SQR(pt[k].y);
  193.         c[0][0] +=  px2 *px2;    /* coefficients for normal equations */
  194.         c[0][1] +=  px2 *py2;
  195.         c[0][2] +=  px2 *pt[k].x;
  196.         c[0][3] +=  px2 *pt[k].y;
  197.         c[0][4] +=  px2 *pt[k].x *pt[k].y;
  198.         c[0][5] +=  px2;
  199.         c[1][1] +=  py2 *py2;
  200.         c[1][2] +=  pt[k].x *py2;
  201.         c[1][3] +=  pt[k].y *py2;
  202.         c[1][4] +=  pt[k].x *py2 *pt[k].y;
  203.         c[1][5] +=  py2;
  204.         c[2][2] +=  px2;
  205.         c[2][3] +=  pt[k].x *pt[k].y;
  206.         c[2][4] +=  px2 *pt[k].y;
  207.         c[2][5] +=  pt[k].x;
  208.         c[3][3] +=  py2;
  209.         c[3][4] +=  pt[k].x *py2;
  210.         c[3][5] +=  pt[k].y;
  211.         c[4][4] +=  px2 *py2;
  212.         c[4][5] +=  pt[k].x *pt[k].y;
  213.  
  214.         sumd    += zeta[k];
  215.         sumdx   += pt[k].x *zeta[k];
  216.         sumdx2  += px2 *zeta[k];
  217.         sumdy   += pt[k].y *zeta[k];
  218.         sumdy2  += py2 *zeta[k];
  219.         sumdxy  += pt[k].x *pt[k].y *zeta[k];
  220.         }
  221.  
  222.     c[1][0] =c[0][1];  /* Coefficient matrix is symmetric about diagonal */
  223.     c[2][0] =c[0][2];  c[2][1] =c[1][2];
  224.     c[3][0] =c[0][3];  c[3][1] =c[1][3];  c[3][2] =c[2][3];
  225.     c[4][0] =c[0][4];  c[4][1] =c[1][4];  c[4][2] =c[2][4];
  226.                                           c[4][3] =c[3][4];
  227.     c[5][0] =c[0][5];  c[5][1] =c[1][5];  c[5][2] =c[2][5];
  228.                                           c[5][3] =c[3][5];
  229.     c[5][4] =c[4][5];  c[5][5] =npoints;
  230.  
  231.     b[0] =sumdx2;    b[1] =sumdy2;    b[2] =sumdx;  /* new vector */
  232.     b[3] =sumdy;     b[4] =sumdxy;    b[5] =sumd;
  233.  
  234.     printf("Augmented matrix:\n");
  235.     PrintMatrix(c, b, 6);
  236.     ok =Gauss(c, b, 6, xv);
  237.     printf("x = a5*x^2 + a4*y^2 + a3*x + a2*y + a1*x*y + a0:\n");
  238.     PrintSolution(xv, 6, 'X');
  239.  
  240.     sumd = sumdx = sumdx2 = sumdy = sumdy2 = sumdxy =0;
  241.  
  242.     for(k =0; k < npoints; k++) {
  243.         sumd   += eta[k];
  244.         sumdx  += pt[k].x *eta[k];
  245.         sumdx2 += px2 *eta[k];
  246.         sumdy  += pt[k].y *eta[k];
  247.         sumdy2 += py2 *zeta[k];
  248.         sumdxy += pt[k].x *pt[k].y *eta[k];
  249.         }
  250.  
  251.  /* Coefficient matrix must remain unchanged. */
  252.  
  253.     b[0] =sumdx2;    b[1] =sumdy2;    b[2] =sumdx;  /* New vector */
  254.     b[3] =sumdy;     b[4] =sumdxy;    b[5] =sumd;
  255.  
  256.     ok =Gauss(c, b, 6, yv);
  257.     printf("y = b5*x^2 + b4*y^2 + b3*x + b2*y + b1*x*y + b0:\n");
  258.     PrintSolution(yv, 6, 'Y');
  259.  
  260.     printf("Residuals:\n");
  261.     for(k =0; k < npoints; k++) {
  262.         xt = SQR(pt[k].x)*xv[0] + SQR(pt[k].y)*xv[1] +
  263.                  pt[k].x *xv[2] + pt[k].y*xv[3] +
  264.                  pt[k].x *pt[k].y*xv[4] +xv[5];
  265.         xt = zeta[k] -xt;
  266.         yt = SQR(pt[k].x)*yv[0] + SQR(pt[k].y)*yv[1] +
  267.                  pt[k].x*yv[2] + pt[k].y*yv[3] +
  268.                  pt[k].x *pt[k].y*yv[4] +yv[5];
  269.         yt = eta[k] -yt;
  270.         printf("%4d  %12.6f  %12.6f\n", (k+1), xt, yt);
  271.         }
  272. }
  273.