Linear Algebra and the C Language/a04w

__NOTOC__ Install and compile this file in your working directory. /* ------------------------------------ */ /* Save as: c00a.c */ /* ------------------------------------ */ #include "v_a.h" /* ------------------------------------ */ /* ------------------------------------ */ #define RA R5 #define CA C5 #define Cb C1 /* ------------------------------------ */ /* ------------------------------------ */ int main(void) { double ab[RA*(CA+Cb)]={ /* x**4 x**3 x**2 x**1 x**0 */ 1, 1., 1., 1., 1., -5., 16., 8., 4., 2., 1, 8., 81., 27., 9., 3., 1, -7., 256., 64., 16., 4., 1, 1., 625., 125., 25., 5., 1, -4., }; 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 **A_T = transpose_mR(A, i_mR(CA,RA)); double **A_TA = mul_mR(A_T,A, i_mR(CA,CA)); double **invA_TA = inv_mR(A_TA, i_mR(CA,CA)); double **invA_TAA_T = mul_mR(invA_TA,A_T, i_mR(CA,RA)); double **x = mul_mR(invA_TAA_T,b, i_mR(CA,Cb)) ; clrscrn(); printf(" Fitting a Quartic equation Curve to Data:\n\n"); printf(" A:"); p_mR(A, S10,P2,C7); printf(" b:"); p_mR(b, S10,P2,C7); printf(" Ab:"); p_mR(Ab, S10,P2,C7); stop(); clrscrn(); printf(" A_T:"); p_mR(A_T, S10,P2,C7); printf(" A_TA:"); p_mR(A_TA, S10,P2,C7); printf(" inv(A_TA):"); p_mR(invA_TA, S10,P4,C7); stop(); clrscrn(); printf(" inv(A_TA)A_T:"); p_mR(invA_TAA_T, S10,P4,C7); printf(" x = inv(A_TA)A_T b:"); p_mR(x, S10,P4,C7); stop(); clrscrn(); printf(" x = inv(A_TA)A_T b:"); p_mR(x, S10,P4,C7); printf(" The Quartic equation Curve to Data: \n\n" " y = %+.3f*x^4 %+.3f*x^3 %+.3f*x^2 %+.3f*x %+.3f\n\n" ,x[R1][C1],x[R2][C1],x[R3][C1],x[R4][C1],x[R5][C1]); stop(); f_mR(A); f_mR(b); f_mR(Ab); f_mR(A_T); f_mR(A_TA); f_mR(invA_TA); f_mR(invA_TAA_T); f_mR(x); 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 Pinv_Rn_mR() function to solve the system that will give us the coefficients a, b, c, d, e Screen output example: Fitting a Quartic equation Curve to Data : A : +1.00 +1.00 +1.00 +1.00 +1.00 +16.00 +8.00 +4.00 +2.00 +1.00 +81.00 +27.00 +9.00 +3.00 +1.00 +256.00 +64.00 +16.00 +4.00 +1.00 +625.00 +125.00 +25.00 +5.00 +1.00 b : -5.00 +8.00 -7.00 +1.00 -4.00 Ab : +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. A_T : +1.00 +16.00 +81.00 +256.00 +625.00 +1.00 +8.00 +27.00 +64.00 +125.00 +1.00 +4.00 +9.00 +16.00 +25.00 +1.00 +2.00 +3.00 +4.00 +5.00 +1.00 +1.00 +1.00 +1.00 +1.00 A_TA : +462979.00 +96825.00 +20515.00 +4425.00 +979.00 +96825.00 +20515.00 +4425.00 +979.00 +225.00 +20515.00 +4425.00 +979.00 +225.00 +55.00 +4425.00 +979.00 +225.00 +55.00 +15.00 +979.00 +225.00 +55.00 +15.00 +5.00 inv(A_TA) : +0.1215 -1.4583 +6.0243 -9.8958 +5.2500 -1.4583 +17.5694 -72.9167 +120.3889 -64.1667 +6.0243 -72.9167 +304.3299 -505.7292 +271.2500 -9.8958 +120.3889 -505.7292 +847.1528 -458.3333 +5.2500 -64.1667 +271.2500 -458.3333 +251.0000 Press return to continue. inv(A_TA)*A_T : +0.0417 -0.1667 +0.2500 -0.1667 +0.0417 -0.5833 +2.1667 -3.0000 +1.8333 -0.4167 +2.9583 -9.8333 +12.2500 -6.8333 +1.4583 -6.4167 +17.8333 -19.5000 +10.1667 -2.0833 +5.0000 -10.0000 +10.0000 -5.0000 +1.0000 x = inv(A_TA)*A_T*b : -3.6250 +44.7500 -191.8750 +329.7500 -184.0000 Press return to continue. x = inv(A_TA)*A_T*b : -3.6250 +44.7500 -191.8750 +329.7500 -184.0000 The Quartic equation Curve to Data : y = -3.625*x^4 +44.750*x^3 -191.875*x^2 +329.750*x -184.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)