Linear Algebra and the C Language/a040

__NOTOC__ Install and compile this file in your working directory. /* ------------------------------------ */ /* Save as: c00f.c */ /* ------------------------------------ */ #include "v_a.h" /* ------------------------------------ */ /* ------------------------------------ */ #define RA R5 #define CA C5 /* ------------------------------------ */ #define FACTOR_E +1.E-2 /* ------------------------------------ */ int main(void) { double tA[RA*CA]={ /* x**4 x**3 x**2 x**1 x**0 */ +1, +1, +1, +1, +1, +16, +8, +4, +2, +1, +81, +27, +9, +3, +1, +256, +64, +16, +4, +1, +625, +125, +25, +5, +1, }; double tb[RA*C1]={ /* y */ -2, -2, 3, -9, 4, }; double **A = ca_A_mR(tA, i_mR(RA,CA)); double **b = ca_A_mR(tb, i_mR(RA,C1)); double **Pinv = i_mR(CA,RA); double **Pinvb = i_mR(CA,C1); clrscrn(); printf(" Fitting a linear Curve to Data:\n\n"); printf(" A:"); p_mR(A, S10,P2,C7); printf(" b:"); p_mR(b, S10,P2,C7); stop(); clrscrn(); printf(" Pinv = V invS_T U_T "); Pinv_Rn_mR(A,Pinv,FACTOR_E); pE_mR(Pinv, S12,P4,C10); printf(" Pinv b "); mul_mR(Pinv,b,Pinvb); p_mR(Pinvb, S10,P4,C10); printf(" The Quartic equation Curve to Data: \n\n" " y = %+.9f*x^4 %+.9f*x^3 %+.9f*x^2 %+.9f*x %+.9f\n\n" ,Pinvb[R1][C1],Pinvb[R2][C1],Pinvb[R3][C1], Pinvb[R4][C1],Pinvb[R5][C1]); stop(); f_mR(b); f_mR(A); f_mR(Pinv); f_mR(Pinvb); 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 linear 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 : -2.00 -2.00 +3.00 -9.00 +4.00 Press return to continue. Pinv = V * invS_T * U_T +4.1667e-02 -1.6667e-01 +2.5000e-01 -1.6667e-01 +4.1667e-02 -5.8333e-01 +2.1667e+00 -3.0000e+00 +1.8333e+00 -4.1667e-01 +2.9583e+00 -9.8333e+00 +1.2250e+01 -6.8333e+00 +1.4583e+00 -6.4167e+00 +1.7833e+01 -1.9500e+01 +1.0167e+01 -2.0833e+00 +5.0000e+00 -1.0000e+01 +1.0000e+01 -5.0000e+00 +1.0000e+00 Pinv * b +2.6667 -30.3333 +117.8333 -181.1666 +89.0000 The Quartic equation Curve to Data : y = +2.666665810*x^4 -30.333323719*x^3 +117.833298608*x^2 -181.166625268*x +88.999995106 Press return to continue. Copy and paste in Octave: function y = f (x) y = +2.666665810*x^4 -30.333323719*x^3 +117.833298608*x^2 -181.166625268*x +88.999995106; endfunction f (+1) f (+2) f (+3) f (+4) f (+5)