Linear Algebra and the C Language/a045


Install and compile this file in your working directory.

/* ------------------------------------ */
/*  Save as:   c00d.c                   */
/* ------------------------------------ */
#include "v_a.h"
/* ------------------------------------ */
/* ------------------------------------ */
#define   RA R5
#define   CA C5
/* ------------------------------------ */
#define FACTOR_E        +1.E-0         
/* ------------------------------------ */
int main(void)
{
double  txy[8] ={
   1,       2,
   2,      -8,
   3,      -8,
   4,      -3  };
   
double tA[RA*CA]={
/* x**2    y**2    x       y       e     */
  +1,     +0,     +0,     +0,     +0,      
  +1,     +4,     +1,     +2,     +1,        
  +4,    +64,     +2,     -8,     +1,        
  +9,    +64,     +3,     -8,     +1,        
 +16,     +9,     +4,     -3,     +1,                    
};

double tb[RA*C1]={
/*    = 0   */
       +1,   
       +0,   
       +0,   
       +0,   
       +0 
};

double **xy      =    ca_A_mR(txy,    i_mR(R4,C2));
double **A       =    ca_A_mR(tA,     i_mR(RA,CA));
double **b       =    ca_A_mR(tb,     i_mR(RA,C1));
double **Pinv    = Pinv_Rn_mR(A,      i_mR(CA,RA),FACTOR_E);            
double **Pinvb   =     mul_mR(Pinv,b, i_mR(CA,C1));         

  clrscrn();
  printf("\n");
  printf(" Find the coefficients a, b, c, d, e, of the curve \n\n");
  printf("     ax**2 + by**2 + cx + dy + e  = 0 \n\n");
  printf(" that passes through these four points.\n\n");
  printf("         x          y");
  p_mR(xy, S10,P0,C6);
  stop();
  
  clrscrn();
  printf(" Using the given points, we obtain this matrix.\n");
  printf("  (a = 1. This is my choice)\n\n");  
  printf(" A:");
  p_mR(A, S10,P2,C7);
  printf(" b:");
  p_mR(b, S10,P2,C7);
   
  printf(" Pinv = V invS_T U_T ");
  pE_mR(Pinv, S12,P4,C10); 
  stop();
  
  clrscrn();
  printf(" Solving this system yields a unique\n"
         " least squares solution, namely   \n\n");  
  printf(" Pinv b ");   
  p_mR(Pinvb, S10,P4,C10);   
  printf(" The coefficients a, b, c, d, e, of the curve are: \n\n"
         "  %+.2f*x^2 %+.2f*y^2 %+.2f*x %+.2f*y %+.2f = 0\n\n"
            ,Pinvb[R1][C1],Pinvb[R2][C1],Pinvb[R3][C1],
             Pinvb[R4][C1],Pinvb[R5][C1]);      
  stop();  

  f_mR(xy);
  f_mR(b);
  f_mR(A);
  f_mR(Pinv);
  f_mR(Pinvb); 

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


Screen output example:
  Find the coefficients a, b, c, d, e, of the curve 

     ax**2 + by**2 + cx + dy + e  = 0 

 that passes through these four points.

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

 Press return to continue. 


 Using the given points, we obtain this matrix.
  (a = 1. This is my choice)

 A :
     +1.00      +0.00      +0.00      +0.00      +0.00 
     +1.00      +4.00      +1.00      +2.00      +1.00 
     +4.00     +64.00      +2.00      -8.00      +1.00 
     +9.00     +64.00      +3.00      -8.00      +1.00 
    +16.00      +9.00      +4.00      -3.00      +1.00 

 b :
     +1.00 
     +0.00 
     +0.00 
     +0.00 
     +0.00 

 Pinv = V * invS_T * U_T 
 +1.0000e+00  +8.4318e-12  +5.5131e-11  -5.7896e-11  +1.5745e-11 
 +4.0000e-02  +2.0000e-02  -8.0000e-02  +1.0000e-01  -4.0000e-02 
 -5.0000e+00  -3.7921e-11  -1.0000e+00  +1.0000e+00  -7.7166e-11 
 +4.0000e-02  +2.2000e-01  -6.8000e-01  +7.0000e-01  -2.4000e-01 
 +3.7600e+00  +4.8000e-01  +2.6800e+00  -2.8000e+00  +6.4000e-01 

 Press return to continue. 


 Solving this system yields a unique
 least squares solution, namely   

 Pinv * b 
   +1.0000 
   +0.0400 
   -5.0000 
   +0.0400 
   +3.7600 

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

  +1.00*x^2 +0.04*y^2 -5.00*x +0.04*y +3.76 = 0

 Press return to continue.


Copy and paste in Octave:
function xy = f (x,y)
  xy = +1.00*x^2 +0.04*y^2 -5.00*x +0.04*y +3.76;
endfunction

f (+1,+2) 
f (+2,-8)
f (+3,-8) 
f (+4,-3)