home *** CD-ROM | disk | FTP | other *** search
/ InfoMagic Source Code 1993 July / THE_SOURCE_CODE_CD_ROM.iso / gnu / nihcl-3.0 / ex / ex5-11.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-05-15  |  2.7 KB  |  97 lines
  1. // ex5-11.c -- Using automatic derivatives
  2. //           to solve systems of equations
  3.  
  4. // $Header: /afs/alw.nih.gov/unix/sun4_40c/usr/local/src/nihcl-3.0/share/ex/RCS/ex5-11.c,v 3.0 90/05/15 22:45:14 kgorlen Rel $
  5.  
  6. #include <iostream.h>
  7. #include <iomanip.h>
  8. #include "ArrayPartial.h"
  9. #include "Matrix.h"
  10.  
  11. void print_iteration(int count,const ArrayPartial& X)
  12. {
  13.     if ( count==0 )
  14.       cout << "INITIAL VALUES" << endl;
  15.     else
  16.       cout << "VALUES FROM ITERATION #" << count << endl;
  17.     cout << setiosflags(ios::showpos|ios::uppercase|
  18.                       ios::scientific) 
  19.          << setprecision(11);
  20.     cout << "X = " << X[1][0] << endl
  21.          << "Y = " << X[2][0] << endl
  22.          << "Z = " << X[3][0] << endl
  23.          << endl;
  24. }
  25. void print_function(const ArrayPartial& Y)
  26. {
  27.     cout << setiosflags(ios::showpos|ios::uppercase|
  28.                       ios::scientific) 
  29.          << setprecision(11);
  30.     cout << "FUNCTION VALUES" << endl
  31.          << "F(X,Y,Z) = " << Y[1][0] << endl
  32.          << "G(X,Y,Z) = " << Y[2][0] << endl
  33.          << "H(X,Y,Z) = " << Y[3][0] << endl
  34.          << endl;
  35. }
  36.  
  37. void print_correction(const Matrix& v)
  38. {
  39.     cout << setiosflags(ios::showpos|ios::uppercase|
  40.                       ios::scientific) 
  41.          << setprecision(11);
  42.     cout << "CORRECTION VECTOR" << endl
  43.          << "X = " << v(0,0) << endl
  44.          << "Y = " << v(1,0) << endl
  45.          << "Z = " << v(2,0) << endl
  46.          << endl;
  47. }
  48.  
  49. // NVAR = number of independent variables
  50. static const int NVAR = 3;
  51. // NEQN = number of equations
  52. static const int NEQN = 3;
  53. // NITER = number of iterations
  54. static const int NITER = 16;
  55. // DELTA = convergence threshhold
  56. static const double DELTA =1.0E-12;
  57. // x0 = initial values of independent variables
  58. static double x0[] = { 1, 1, 1 };
  59.  
  60. main()
  61. {
  62.   // X[i] is the i-th variable
  63.   // (X[1][0],X[2][0],X[3][0]) is the solution vector
  64.   // initial value is (1.0,1.0,1.0)
  65.   ArrayPartial X(NVAR, x0);
  66.  
  67.   // correction vector
  68.   Matrix del(NVAR/*rows*/,1/*columns*/);
  69.  
  70.   int count =0;
  71.   while(count <=NITER) {
  72.     print_iteration(count++,X);
  73.  
  74.     // define 3 equations in 3 variables
  75.     // (f[1][0],f[2][0],f[3][0]) is the function value vector
  76.     // ( f[i][j] ) is the jacobian matrix
  77.     ArrayPartial f(NEQN);
  78.     f[1]= 16.0*X[1].pow(4) + 16.*X[2].pow(4) + X[3].pow(4) - 16.0;
  79.     f[2]=      X[1].pow(2) +     X[2].pow(2) + X[3].pow(2) -  3.0;
  80.     f[3]=      X[1].pow(3) -     X[2];
  81.  
  82.     print_function(f);
  83.  
  84.     // recompute jacobian and new correction vector
  85.     Matrix J = f.jacobian();
  86.     Matrix v = f.value();
  87.     del = -(~J*v);
  88.  
  89.     print_correction(del);
  90.     if ( norm(del) < DELTA ) break;
  91.  
  92.     // compute next estimate for solution
  93.     // (X[1][0],X[2][0],X[3][0])
  94.     X += del;
  95.     }
  96. }
  97.