/******************************************************** * Approximation of a discrete real function F(x) by * * least squares * * ----------------------------------------------------- * * Ref.: "Méthodes de calcul numérique, Tome 2 by Claude * * Nowakowski, PSI Edition, 1984" [BIBLI 04]. * * ----------------------------------------------------- * * SAMPLE RUN: * * * * Number of points : 11 * * * * Degree of polynomial: 3 * * * * Function to approximate: * * X(1), Y(1) = 0 0 * * X(2), Y(2) = 0.1 0.2 * * X(3), Y(3) = 0.2 0.4 * * X(4), Y(4) = 0.3 0.6 * * X(5), Y(5) = 0.4 0.8 * * X(6), Y(6) = 0.5 1 * * X(7), Y(7) = 0.6 0.8 * * X(8), Y(8) = 0.7 0.6 * * X(9), Y(9) = 0.8 0.4 * * X(10), Y(10) = 0.9 0.2 * * X(11), Y(11) = 1 0 * * * * Polynomial approximation of degree 3 (11 points) * * Coefficients of polynomial: * * A(0) = -0.069930070 * * A(1) = 3.496503497 * * A(2) = -3.496503497 * * A(3) = 0.000000000 * * * * Approximated function: * * X Y * * 0.000000 -0.069930 * * 0.100000 0.244755 * * 0.200000 0.489510 * * 0.300000 0.664336 * * 0.400000 0.769231 * * 0.500000 0.804196 * * 0.600000 0.769231 * * 0.700000 0.664336 * * 0.800000 0.489510 * * 0.900000 0.244755 * * 1.000000 -0.069930 * * * * C++ version by J-P Moreau, Paris. * * (www.jpmoreau.fr) * ********************************************************/ #include #include #define SIZE 25 int i,ij,j,k,n,n1,m,m1,m2; double C[SIZE][SIZE]; double A[SIZE],B[SIZE],X[SIZE],Xc[SIZE],Y[SIZE],Yx[SIZE]; double p,xx,s,yc; double IntPower(double x,int k) { if (x==0) return 0; else return (exp(k*log(x))); } void main() { printf("\n Number of points : "); scanf("%d", &n); n--; printf("\n Degree of polynomial: "); scanf("%d", &m); n1=n+1; m1=m+1; m2=m+2; printf("\n Function to approximate:\n"); for (i=1; i<=n1; i++) { printf(" X(%d), Y(%d) = ", i, i); scanf("%lf %lf", &X[i], &Y[i]); } for (k=1; k<=m2; k++) { Xc[k]=0; for (i=1; i<=n1; i++) Xc[k] += IntPower(X[i],k); } yc=0; for (i=1; i<=n1; i++) yc += Y[i]; for (k=1; k<=m; k++) { Yx[k]=0; for (i=1; i<=n1; i++) Yx[k] += Y[i]*IntPower(X[i],k); } for (i=1; i<=m1; i++) for (j=1; j<=m1; j++) { ij=i+j-2; if (i==1 && j==1) C[1][1]= n1; else C[i][j]=Xc[ij]; } B[1]=yc; for (i=2; i<=m1; i++) B[i]=Yx[i-1]; for (k=1; k<=m; k++) for (i=k+1; i<=m1; i++) { B[i] -= C[i][k]/C[k][k]*B[k]; for (j=k+1; j<=m1; j++) C[i][j] -= C[i][k]/C[k][k]*C[k][j]; } A[m1]=B[m1]/C[m1][m1]; for (i=m; i>0; i--) { s=0; for (k=i+1; k<=m1; k++) s = s + C[i][k]*A[k]; A[i] = (B[i]-s)/C[i][i]; } printf("\n Polynomial approximation of degree %d (%d points)\n",m,n+1); printf(" Coefficients of polynomial:\n"); for (i=1; i<=m1; i++) printf(" A(%d) = %15.9f\n", i-1, A[i]); printf("\n Approximated function:\n"); printf(" X Y\n"); for (i=1; i<=n1; i++) { xx=X[i]; p=0; for (k=1; k<=m1; k++) p = p*xx + A[m1+1-k]; printf(" %11.6f %11.6f\n", xx, p); } printf("\n\n"); } // end of file approx.cpp