/*********************************************************************** * * * Interpolate a function by trigonometric polynoms * * * * -------------------------------------------------------------------- * * SAMPLE RUN: * * * * Number of points: 3 * * * * Number of harmonics: 3 * * * * * * Function to interpolate: * * X Y * * 0.000000 3.141593 * * 0.897598 0.897598 * * 1.795196 1.795196 * * 2.692794 2.692794 * * 3.590392 3.590392 * * 4.487990 4.487990 * * 5.385587 5.385587 * * 6.283185 3.141593 * * * * Fourier coefficients: * * A(0) = 3.141593 B(0) = 0.000000 * * A(1) = -0.000000 B(1) = -1.863881 * * A(2) = -0.000000 B(2) = -0.715810 * * A(3) = -0.000000 B(3) = -0.204871 * * * * Interpolation: * * X Y * * 0.000000 3.141593 * * 0.448799 1.573507 * * 0.897598 0.897598 * * 1.346397 1.174039 * * 1.795196 1.795196 * * 2.243995 2.293325 * * 2.692794 2.692794 * * 3.141593 3.141593 * * 3.590392 3.590392 * * 4.039191 3.989860 * * 4.487990 4.487990 * * 4.936788 5.109147 * * 5.385587 5.385587 * * 5.834386 4.709678 * * 6.283185 3.141593 * * * * -------------------------------------------------------------------- * * "Méthodes de calcul numérique, Tome 2 By Claude Nowakowski, * * PSI Edition, 1984" [BIBLI 04]. * * * * C++ Release 1.0 By J-P Moreau, Paris. * * (www.jpmoreau.fr) * ***********************************************************************/ #include #include #define SIZE 50 #define PI 3.1415926535 // Labels: e10, e90, e100, e200 int IP, K, L, M, N, NH, NN; double A[SIZE], B[SIZE], TF[SIZE]; double C1, CF, CP, CT, F0, S1, SP, X; double Q, S, U0, U1, U2; void main() { e10:printf("\n Number of points : "); scanf("%d", &N); printf("\n Number of harmonics: "); scanf("%d", &NH); printf("\n"); if (N < NH) { printf(" Error, N must be >= NH.\n"); goto e10; } NN = 2 * N + 1; // TF stores function values for (K=0; K<=NN; K++) TF[K] = K * 2.0 * PI / NN; TF[0] = (TF[0] + TF[NN]) / 2.0; TF[NN] = TF[0]; CF = 2.0 / NN; printf("\n Function to interpolate:\n"); printf(" X Y\n"); for (K=0; K<=NN; K++) { X = K * 2.0 * PI / NN; printf("%11.6f%11.6f\n", X, TF[K]); } CT = PI * CF; // Initialize S1 = sin(CT); C1 = cos(CT); CP = 1.0; SP = 0.0; IP = 0; F0 = TF[0]; // Calculate Fourier coefficients A, B, for each harmonic e90: U2 = 0.0; U1 = 0.0; M = NN - 1; e100:U0 = TF[M] + 2.0 * CP * U1 - U2; U2 = U1; U1 = U0; M--; if (M > 0) goto e100; A[IP] = CF * (F0 + CP * U1 - U2); B[IP] = CF * SP * U1; if (IP > NH) goto e200; Q = C1 * CP - S1 * SP; SP = C1 * SP + S1 * CP; CP = Q; IP++; goto e90; e200:A[0] = 0.5 * A[0]; printf("\n Fourier coefficients:\n"); for (K=0; K<=NH; K++) printf(" A(%d) = %11.6f B(%d) = %11.6f\n", K, A[K], K, B[K]); // Interpolate function printf("\n Interpolation:\n"); printf(" X Y\n"); for (K=0; K<=NN*2; K++) { X = K * 2.0 * PI / (2 * NN); S = A[0]; for (L=1; L<=NH; L++) S += A[L] * cos(X * L) + B[L] * sin(X * L); printf("%11.6f%11.6f\n", X, S); } printf("\n"); } // end of main program // End of file trigpoly.cpp