Linear Algebra and the C Language/a0bl


Install and compile this file in your working directory. 

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

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 \n\n"
            ,x[R1][C1],x[R3][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    -8 
   -2    +8 
   +2    -8 
   +5    +8 

 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   -8.00 
  -8.00   +4.00   -2.00   +1.00   +8.00 
  +8.00   +4.00   +2.00   +1.00   -8.00 
+125.00  +25.00   +5.00   +1.00   +8.00 

 Press return to continue. 


 Q :
   -0.7057    +0.6982    +0.0452    -0.1117 
   -0.0452    +0.1117    -0.7057    +0.6982 
   +0.0452    +0.1117    +0.7057    +0.6982 
   +0.7057    +0.6982    -0.0452    -0.1117 

 R :
 +177.1384    +0.0000    +7.2373    +0.0000 
   +0.0000   +35.8050    +0.0000    +1.6199 
   -0.0000    +0.0000    +2.3710    +0.0000 
   -0.0000    +0.0000    +0.0000    +1.1730 

 Press return to continue. 


 Q_T :
 -7.0566e-01  -4.5162e-02  +4.5162e-02  +7.0566e-01 
 +6.9823e-01  +1.1172e-01  +1.1172e-01  +6.9823e-01 
 +4.5162e-02  -7.0566e-01  +7.0566e-01  -4.5162e-02 
 -1.1172e-01  +6.9823e-01  +6.9823e-01  -1.1172e-01 

 invR :
 +5.6453e-03  +0.0000e+00  -1.7232e-02  +0.0000e+00 
 +0.0000e+00  +2.7929e-02  +0.0000e+00  -3.8569e-02 
 -0.0000e+00  -0.0000e+00  +4.2176e-01  -0.0000e+00 
 -0.0000e+00  -0.0000e+00  -0.0000e+00  +8.5250e-01 

 Press return to continue. 


 Solving this system yields a unique
 least squares solution, namely   

 x = invR * Q_T * b :
     +0.27 
     +0.00 
     -5.07 
     +0.00 


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

 y = +0.27x**3   -5.07x 

 Press return to continue.



Copy and paste in Octave:

function y = f (x)
  y = +0.266666667*x^3  -5.066666667*x;
endfunction

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