Linear Algebra and the C Language/a07l


Install and compile this file in your working directory.

/* ------------------------------------ */
/*  Save as :   c00e.c                  */
/* ------------------------------------ */
#include   "v_a.h"
#include "dpoly.h"
/* -----------------***---------------- */
int main(void)
{
double   xy[10] ={
   1,      -5,
   2,       8,
   3,      -7,
   4,       1,
   5,      -4     };

double **XY =  ca_A_mR(xy,i_mR(R5,C2));
double **A  =             i_mR(R5,C5);
double **b =              i_mR(R5,C1);
double **Ab =   i_Abr_Ac_bc_mR(R5,C5,C1);

  clrscrn();
  printf(" Find the coefficients a, b, c  of the curve \n\n");
  printf("   y =  ax**4 + bx**3 + cx**2 + dx + e       \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\n");
  printf("   x**4    x**3    x**2    x**1    x**0    y");
  i_A_b_with_XY_mR(XY,A,b);
  c_A_b_Ab_mR(A,b,Ab);
  p_mR(Ab,S7,P2,C6);
  stop();

  clrscrn();
  printf(" The Gauss Jordan process will reduce this matrix to : \n");
  gj_TP_mR(Ab);
  p_mR(Ab,S7,P2,C6);
  printf("\n The coefficients a, b, c of the curve are :  \n\n");
  p_eq_poly_mR(Ab);
  stop();

  clrscrn();
  printf("    x     y \n");
  p_mR(XY,S5,P0,C6);
  printf("\n");

  printf(" Verify the result : \n\n");
  verify_X_mR(Ab,XY[R1][C1]);
  verify_X_mR(Ab,XY[R2][C1]);
  verify_X_mR(Ab,XY[R3][C1]);
  verify_X_mR(Ab,XY[R4][C1]);
  verify_X_mR(Ab,XY[R5][C1]);
  printf("\n\n\n");
  stop();

  f_mR(XY);
  f_mR(A);
  f_mR(b);
  f_mR(Ab);

  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 gj_TP_mR() function to solve
     the system that will give us the coefficients a, b, c, d, e
Screen output example:
 Find the coefficients a, b, c  of the curve 

   y =  ax**4 + bx**3 + cx**2 + dx + e       

 that passes through the points.             

    x     y
   +1    -5 
   +2    +8 
   +3    -7 
   +4    +1 
   +5    -4 

 Using the given points, we obtain this matrix

   x**4    x**3    x**2    x**1    x**0    y
  +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. 


 The Gauss Jordan process will reduce this matrix to : 

  +1.00   +0.00   +0.00   +0.00   +0.00   -3.63 
  +0.00   +1.00   +0.00   +0.00   +0.00  +44.75 
  +0.00   +0.00   +1.00   +0.00   +0.00 -191.88 
  -0.00   -0.00   -0.00   +1.00   -0.00 +329.75 
  +0.00   +0.00   +0.00   +0.00   +1.00 -184.00 


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

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


 Press return to continue. 


    x     y 

   +1    -5 
   +2    +8 
   +3    -7 
   +4    +1 
   +5    -4 


 Verify the result : 

 With x =  +1.000,       y = -5.000 
 With x =  +2.000,       y = +8.000 
 With x =  +3.000,       y = -7.000 
 With x =  +4.000,       y = +1.000 
 With x =  +5.000,       y = -4.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)