Linear Algebra and the C Language/a0bd


Install and compile this file in your working directory. 

/* ------------------------------------ */
/*  Save as :   c00d.c                  */
/* ------------------------------------ */
#include "v_a.h"
/* ------------------------------------ */
/* ------------------------------------ */
#define   RA R5
#define   CA C5
#define   Cb C1 
/* ------------------------------------ */
/* ------------------------------------ */
int main(void)
{
double xy[6] ={
   10,     10,
   -5,      1,
    7,    -10    };
   
double ab[RA*(CA+Cb)]={
/* x**2    y**2    x       y            =  0   */
  +1,     +0,     +0,     +0,     +0,     +1,   
  +0,     +1,     +0,     +0,     +0,     +1,   
+100,   +100,    +10,    +10,     +1,     +0,   
 +25,     +1,     -5,     +1,     +1,     +0,   
 +49,   +100,     +7,    -10,     +1,     +0,            
}; 

double **XY = ca_A_mR(xy,i_mR(R3,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, d  of a circle   \n\n");
  printf("     ax**2 + ay**2 + bx + cy + d  = 0            \n\n");
  printf(" that passes through these three XY.             \n\n");
  printf("    x     y");
  p_mR(XY,S5,P0,C6);
  printf(" Using the given points, we obtain this matrix.\n");
  printf("  (a = 1. This is my choice)\n\n");  
  printf(" Ab :\n");
  printf("   x**2    y**2    x       y            =  0     ");
  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(" The coefficients a, b, c, d  of the curve are : \n\n"
         "  %+.2fx**2 %+.2fy**2 %+.2fx %+.2fy %+.2f = 0\n\n"
            ,x[R1][C1],x[R2][C1],x[R3][C1],x[R4][C1],x[R5][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;
}
/* ------------------------------------ */
/* ------------------------------------ */



Screen output example:

 Find the coefficients a, b, c, d  of a circle   

     ax**2 + ay**2 + bx + cy + d  = 0            

 that passes through these three XY.             

    x     y
  +10   +10 
   -5    +1 
   +7   -10 

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

 Ab :
   x**2    y**2    x       y            =  0     
  +1.00   +0.00   +0.00   +0.00   +0.00   +1.00 
  +0.00   +1.00   +0.00   +0.00   +0.00   +1.00 
+100.00 +100.00  +10.00  +10.00   +1.00   +0.00 
 +25.00   +1.00   -5.00   +1.00   +1.00   +0.00 
 +49.00 +100.00   +7.00  -10.00   +1.00   +0.00 

 Press return to continue. 


 Q :
   +0.0088    -0.0213    +0.0046    -0.6732    -0.7391 
   +0.0000    +0.0186    -0.0191    +0.7387    -0.6735 
   +0.8761    -0.2704    +0.3985    +0.0196    +0.0028 
   +0.2190    -0.5131    -0.8298    +0.0014    +0.0109 
   +0.4293    +0.8141    -0.3901    -0.0270    +0.0039 

 R :
 +114.1359  +130.7652   +10.6715    +4.6874    +1.5245 
   -0.0000   +53.8745    +5.5600   -11.3588    +0.0306 
   -0.0000    -0.0000    +5.4043    +7.0561    -0.8214 
   -0.0000    -0.0000    -0.0000    +0.4678    -0.0060 
   -0.0000    +0.0000    +0.0000    +0.0000    +0.0176 

 Press return to continue. 


 Q_T :
 +8.7615e-03  +0.0000e+00  +8.7615e-01  +2.1904e-01  +4.2931e-01 
 -2.1266e-02  +1.8562e-02  -2.7044e-01  -5.1309e-01  +8.1413e-01 
 +4.5779e-03  -1.9096e-02  +3.9854e-01  -8.2984e-01  -3.9005e-01 
 -6.7319e-01  +7.3872e-01  +1.9620e-02  +1.4100e-03  -2.7021e-02 
 -7.3910e-01  -6.7349e-01  +2.7660e-03  +1.0935e-02  +3.8595e-03 

 invR :
 +8.7615e-03  -2.1266e-02  +4.5779e-03  -6.7319e-01  -7.3910e-01 
 -0.0000e+00  +1.8562e-02  -1.9096e-02  +7.3872e-01  -6.7349e-01 
 -0.0000e+00  +0.0000e+00  +1.8504e-01  -2.7909e+00  +7.7024e+00 
 -0.0000e+00  +0.0000e+00  -0.0000e+00  +2.1376e+00  +7.2929e-01 
 -0.0000e+00  +0.0000e+00  -0.0000e+00  +0.0000e+00  +5.6945e+01 

 Press return to continue. 


 Solving this system yields a unique
 least squares solution, namely   

 x = invR * Q_T * b :
     +1.00 
     +1.00 
    -11.07 
     -0.89 
    -80.44 

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

  +1.00x**2 +1.00y**2 -11.07x -0.89y -80.44 = 0

 Press return to continue.



Copy and paste in Octave:

function xy = f (x,y)
  xy = +1.000000000*x^2 +1.000000000*y^2 -11.065934068*x -0.890109889*y -80.439560439;
endfunction

f (+10,+10)
f (-5,+1)
f (+7,-10)
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.