Linear Algebra and the C Language/a04w


Install and compile this file in your working directory. 

/* ------------------------------------ */
/*  Save as:   c00a.c                   */
/* ------------------------------------ */
#include "v_a.h"
/* ------------------------------------ */
/* ------------------------------------ */
#define   RA R5
#define   CA C5
#define   Cb C1 
/* ------------------------------------ */
/* ------------------------------------ */
int main(void)
{
double ab[RA*(CA+Cb)]={
/*  x**4   x**3   x**2    x**1    x**0  */
    1,     1.,    1.,     1.,     1., -5.,  
  16.,     8.,    4.,     2.,     1,   8.,                      
  81.,    27.,    9.,     3.,     1,  -7.,                     
 256.,    64.,   16.,     4.,     1,   1.,                   
 625.,   125.,   25.,     5.,     1,  -4.,                 
};

double **Ab         =      ca_A_mR(ab, i_Abr_Ac_bc_mR(RA,CA,Cb));
double **A          =    c_Ab_A_mR(Ab,           i_mR(RA,CA));
double **b          =    c_Ab_b_mR(Ab,           i_mR(RA,Cb));

double **A_T        = transpose_mR(A,            i_mR(CA,RA));
double **A_TA       =       mul_mR(A_T,A,        i_mR(CA,CA));  
double **invA_TA    =       inv_mR(A_TA,         i_mR(CA,CA));  
double **invA_TAA_T =       mul_mR(invA_TA,A_T,  i_mR(CA,RA));  

double **x          =       mul_mR(invA_TAA_T,b, i_mR(CA,Cb)) ;  

  clrscrn();
  printf(" Fitting a Quartic equation Curve to Data:\n\n");
  printf(" A:");
  p_mR(A, S10,P2,C7);
  
  printf(" b:");
  p_mR(b, S10,P2,C7);
  
  printf(" Ab:");
  p_mR(Ab, S10,P2,C7);
  stop();
  
  clrscrn();
  printf(" A_T:");
  p_mR(A_T, S10,P2,C7);
  
  printf(" A_TA:");
  p_mR(A_TA, S10,P2,C7);
  
  printf(" inv(A_TA):");
  p_mR(invA_TA, S10,P4,C7);
  stop();
  
  clrscrn();  
  printf(" inv(A_TA)A_T:");
  p_mR(invA_TAA_T, S10,P4,C7);
  
  printf(" x = inv(A_TA)A_T b:");
  p_mR(x, S10,P4,C7);
  stop();
  
  clrscrn();
  printf(" x = inv(A_TA)A_T b:");
  p_mR(x, S10,P4,C7); 
  
  printf(" The Quartic equation Curve to Data: \n\n"
         "  y = %+.3f*x^4 %+.3f*x^3 %+.3f*x^2 %+.3f*x %+.3f\n\n"
            ,x[R1][C1],x[R2][C1],x[R3][C1],x[R4][C1],x[R5][C1]);      
  stop();  
  
  f_mR(A);
  f_mR(b);
  f_mR(Ab);

  f_mR(A_T);
  f_mR(A_TA);       
  f_mR(invA_TA);     
  f_mR(invA_TAA_T);  
    
  f_mR(x); 

  return 0;
}
/* ------------------------------------ */
/* ------------------------------------ */

Presentation :

  Let's calculate the coefficients of a polynomial.
 
              y =  ax**4 + bx**3 + cx**2 + dx + e        
  
         Which passes through these five points.    
          
       x[1],  y[1] 
       x[2],  y[2] 
       x[3],  y[3] 
       x[4],  y[4] 
       x[5],  y[5] 
       
  Using the points we obtain the matrix:

      x**4         x**3        x**2      x**1      x**0     y

     x[1]**4     x[1]**3     x[1]**2   x[1]**1   x[1]**0   y[1]
     x[2]**4     x[2]**3     x[2]**2   x[2]**1   x[2]**0   y[2]
     x[3]**4     x[3]**3     x[3]**2   x[3]**1   x[3]**0   y[3]
     x[4]**4     x[4]**3     x[4]**2   x[4]**1   x[4]**0   y[4]
     x[5]**4     x[5]**3     x[5]**2   x[5]**1   x[5]**0   y[5]

  That we can write:

       x**4          x**3        x**2     x     1    y

      x[1]**4       x[1]**3    x[1]**2   x[1]   1   y[1]
      x[2]**4       x[2]**3    x[2]**2   x[2]   1   y[2]
      x[3]**4       x[3]**3    x[3]**2   x[3]   1   y[3]
      x[4]**4       x[4]**3    x[4]**2   x[4]   1   y[4]
      x[5]**4       x[5]**3    x[5]**2   x[5]   1   y[5]   

     Let's use the Pinv_Rn_mR() function to solve
     the system that will give us the coefficients a, b, c, d, e

Screen output example:

 Fitting a Quartic equation Curve to Data :

 A :
     +1.00      +1.00      +1.00      +1.00      +1.00 
    +16.00      +8.00      +4.00      +2.00      +1.00 
    +81.00     +27.00      +9.00      +3.00      +1.00 
   +256.00     +64.00     +16.00      +4.00      +1.00 
   +625.00    +125.00     +25.00      +5.00      +1.00 

 b :
     -5.00 
     +8.00 
     -7.00 
     +1.00 
     -4.00 

 Ab :
     +1.00      +1.00      +1.00      +1.00      +1.00      -5.00 
    +16.00      +8.00      +4.00      +2.00      +1.00      +8.00 
    +81.00     +27.00      +9.00      +3.00      +1.00      -7.00 
   +256.00     +64.00     +16.00      +4.00      +1.00      +1.00 
   +625.00    +125.00     +25.00      +5.00      +1.00      -4.00 

 Press return to continue. 


 A_T :
     +1.00     +16.00     +81.00    +256.00    +625.00 
     +1.00      +8.00     +27.00     +64.00    +125.00 
     +1.00      +4.00      +9.00     +16.00     +25.00 
     +1.00      +2.00      +3.00      +4.00      +5.00 
     +1.00      +1.00      +1.00      +1.00      +1.00 

 A_TA :
+462979.00  +96825.00  +20515.00   +4425.00    +979.00 
 +96825.00  +20515.00   +4425.00    +979.00    +225.00 
 +20515.00   +4425.00    +979.00    +225.00     +55.00 
  +4425.00    +979.00    +225.00     +55.00     +15.00 
   +979.00    +225.00     +55.00     +15.00      +5.00 

 inv(A_TA) :
   +0.1215    -1.4583    +6.0243    -9.8958    +5.2500 
   -1.4583   +17.5694   -72.9167  +120.3889   -64.1667 
   +6.0243   -72.9167  +304.3299  -505.7292  +271.2500 
   -9.8958  +120.3889  -505.7292  +847.1528  -458.3333 
   +5.2500   -64.1667  +271.2500  -458.3333  +251.0000 

 Press return to continue. 


 inv(A_TA)*A_T :
   +0.0417    -0.1667    +0.2500    -0.1667    +0.0417 
   -0.5833    +2.1667    -3.0000    +1.8333    -0.4167 
   +2.9583    -9.8333   +12.2500    -6.8333    +1.4583 
   -6.4167   +17.8333   -19.5000   +10.1667    -2.0833 
   +5.0000   -10.0000   +10.0000    -5.0000    +1.0000 


 x = inv(A_TA)*A_T*b :
   -3.6250 
  +44.7500 
 -191.8750 
 +329.7500 
 -184.0000 

 Press return to continue. 


 x = inv(A_TA)*A_T*b :
   -3.6250 
  +44.7500 
 -191.8750 
 +329.7500 
 -184.0000 

 The Quartic equation Curve to Data : 

  y = -3.625*x^4 +44.750*x^3 -191.875*x^2 +329.750*x -184.000

 Press return to continue.

Copy and paste in Octave:

function y = f (x)
  y = -3.625*x^4 +44.750*x^3 -191.875*x^2 +329.750*x -184.000;
endfunction

f (+1) 
f (+2)
f (+3) 
f (+4) 
f (+5)
This article is issued from Wikibooks. The text is available under Creative Commons Attribution-Share Alike 4.0 unless otherwise noted. Additional terms may apply for the media files.