Linear Algebra and the C Language/a0bk


Install and compile this file in your working directory. 

/* ------------------------------------ */
/*  Save as :   c00c.c                  */
/* ------------------------------------ */
#include "v_a.h"
/* ------------------------------------ */
/* ------------------------------------ */
#define   RA R4
#define   CA C4
#define   Cb C1 
/* ------------------------------------ */
/* ------------------------------------ */
int main(void)
{
double   xy[8] ={  -5,      -3,
                   -2,       0,
                    2,       3,
                    3,      -2     };
   
double ab[RA*(CA+Cb)]={
/* x**3    x**2     x**1   x**0     y   */
 -125,    +25,     -5,     +1,     -3,   
   -8,     +4,     -2,     +1,     +0,   
   +8,     +4,     +2,     +1,     +3,   
  +27,     +9,     +3,     +1,     -2,        
}; 

double **XY = ca_A_mR(xy,i_mR(R4,C2));

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 **Q    = i_mR(RA,CA);
double **R    = i_mR(CA,CA);

double **invR = i_mR(CA,CA);
double **Q_T  = i_mR(CA,RA);


double **invR_Q_T = i_mR(CA,RA);
double **x        = i_mR(CA,Cb); // x = invR * Q_T * b

  clrscrn();
  printf("\n");
  printf(" Find the coefficients a, b, c  of the curve \n\n");
  printf("      y =  ax**3 + bx**2 + cx + d            \n\n");
  printf(" that passes through the points.             \n\n");
  printf("    x     y");
  p_mR(XY,S5,P0,C6);
  printf(" Using the given points, we obtain this matrix.\n");
  printf("   x**3    x**2    x**1    x**0     y\n");
  p_mR(Ab,S7,P2,C6);
  stop();

    
  clrscrn();
  QR_mR(A,Q,R);    
  printf(" Q :");
  p_mR(Q,S10,P4,C6); 
  printf(" R :");
  p_mR(R,S10,P4,C6);
  stop();

  clrscrn();
  transpose_mR(Q,Q_T);   
  printf(" Q_T :");
  pE_mR(Q_T,S12,P4,C6); 
  inv_mR(R,invR); 
  printf(" invR :");
  pE_mR(invR,S12,P4,C6);
  stop();
  
  clrscrn();
  printf(" Solving this system yields a unique\n"
         " least squares solution, namely   \n\n");
  mul_mR(invR,Q_T,invR_Q_T);
  mul_mR(invR_Q_T,b,x);
  printf(" x = invR * Q_T * b :");
  p_mR(x,S10,P2,C6);
  printf("\n The coefficients a, b, c of the curve are :  \n\n" 
         " y = %+.2fx**3 %+.2fx**2  %+.2fx %+.2f\n\n"
            ,x[R1][C1],x[R2][C1],x[R3][C1],x[R4][C1]);  
  
  stop();

  f_mR(XY);  
  f_mR(A);
  f_mR(b);
  f_mR(Ab);
  f_mR(Q);
  f_mR(Q_T);
  f_mR(R);
  f_mR(invR);  
  f_mR(invR_Q_T); 
  f_mR(x); 

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



Presentation :

         Let's calculate the coefficients of a polynomial.
 
              y =  ax**3 + bx**2 + cx + d        
  
        Which passes through these four points.     
          
       x[1],  y[1] 
       x[2],  y[2] 
       x[3],  y[3] 
       x[4],  y[4] 

   Using the points we obtain the matrix:

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

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

  That we can write:

      x**3       x**2      x      1   y
 
      x[1]**3    x[1]**2   x[1]   1   y[1]
      x[2]**3    x[2]**2   x[2]   1   y[2]
      x[3]**3    x[3]**2   x[3]   1   y[3]
      x[4]**3    x[4]**2   x[4]   1   y[4]
   
     Let's use the QR_mR() function to solve
     the system that will give us the coefficients a, b, c, d



Screen output example:

 Find the coefficients a, b, c  of the curve 

      y =  ax**3 + bx**2 + cx + d            

 that passes through the points.             

    x     y
   -5    -3 
   -2    +0 
   +2    +3 
   +3    -2 

 Using the given points, we obtain this matrix.
   x**3    x**2    x**1    x**0     y

-125.00  +25.00   -5.00   +1.00   -3.00 
  -8.00   +4.00   -2.00   +1.00   +0.00 
  +8.00   +4.00   +2.00   +1.00   +3.00 
 +27.00   +9.00   +3.00   +1.00   -2.00 

 Press return to continue. 


 Q :
   -0.9737    +0.2054    +0.0819    -0.0556 
   -0.0623    +0.1700    -0.9033    +0.3889 
   +0.0623    +0.3529    +0.4209    +0.8333 
   +0.2103    +0.8969    -0.0131    -0.3889 

 R :
 +128.3822   -22.4486    +5.7485    -0.7633 
   +0.0000   +15.2990    +2.0292    +1.6252 
   -0.0000    -0.0000    +2.1995    -0.4136 
   -0.0000    +0.0000    -0.0000    +0.7778 

 Press return to continue. 


 Q_T :
 -9.7365e-01  -6.2314e-02  +6.2314e-02  +2.1031e-01 
 +2.0543e-01  +1.7002e-01  +3.5289e-01  +8.9686e-01 
 +8.1914e-02  -9.0331e-01  +4.2088e-01  -1.3124e-02 
 -5.5556e-02  +3.8889e-01  +8.3333e-01  -3.8889e-01 

 invR :
 +7.7892e-03  +1.1429e-02  -3.0902e-02  -3.2672e-02 
 +0.0000e+00  +6.5364e-02  -6.0304e-02  -1.6865e-01 
 -0.0000e+00  -0.0000e+00  +4.5466e-01  +2.4180e-01 
 -0.0000e+00  -0.0000e+00  -0.0000e+00  +1.2857e+00 

 Press return to continue. 


 Solving this system yields a unique
 least squares solution, namely   

 x = invR * Q_T * b :
     -0.14 
     -0.73 
     +1.31 
     +4.43 


 The coefficients a, b, c of the curve are :  

 y = -0.14x**3 -0.73x**2  +1.31x +4.43

 Press return to continue.



Copy and paste in Octave:

function y = f (x)
  y = -0.139285714*x^3 -0.732142857*x^2  +1.307142857*x +4.428571429;
endfunction

f (-5) 
f (-2)
f (+2) 
f (+3)
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.