'************************************************************************ '* * '* 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]. * '* * '* Basic Release 1.0 By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '************************************************************************ 'Program TrigPoly DEFDBL A-H, O-Z DEFINT I-N PI = 3.1415926535# CLS 10 PRINT INPUT " Number of points: ", N PRINT INPUT " Number of harmonics: ", NH PRINT IF N < NH THEN PRINT " Error, N must be >= NH." GOTO 10 END IF NN = 2 * N + 1 ' TF stores function values DIM TF(NN) FOR K = 0 TO NN TF(K) = K * 2# * PI / NN NEXT K TF(0) = (TF(0) + TF(NN)) / 2# TF(NN) = TF(0) CF = 2# / NN F\$ = "####.######" PRINT PRINT " Function to interpolate:" PRINT " X Y" FOR K = 0 TO NN X = K * 2# * PI / NN PRINT USING F\$; X; PRINT USING F\$; TF(K) NEXT K CT = PI * CF 'Initialize S1 = SIN(CT): C1 = COS(CT) CP = 1#: SP = 0#: IP = 0: F0 = TF(0) DIM A(NH + 1), B(NH + 1)'to store Fourier coeff. 'Calculate Fourier coefficients A, B, for each harmonic 90 U2 = 0#: U1 = 0#: M = NN - 1 100 U0 = TF(M) + 2# * CP * U1 - U2 U2 = U1: U1 = U0: M = M - 1 IF M > 0 THEN GOTO 100 A(IP) = CF * (F0 + CP * U1 - U2) B(IP) = CF * SP * U1 IF IP > NH THEN GOTO 200 Q = C1 * CP - S1 * SP SP = C1 * SP + S1 * CP CP = Q IP = IP + 1 GOTO 90 200 A(0) = .5# * A(0) PRINT PRINT " Fourier coefficients:" FOR K = 0 TO NH PRINT " A("; K; ") = "; PRINT USING F\$; A(K); PRINT " B("; K; ") = "; PRINT USING F\$; B(K) NEXT K 'Interpolate function PRINT PRINT " Interpolation:" PRINT " X Y" FOR K = 0 TO NN * 2 X = K * 2# * PI / (2# * NN) S = A(0) FOR L = 1 TO NH S = S + A(L) * COS(X * L) + B(L) * SIN(X * L) NEXT L PRINT USING F\$; X; PRINT USING F\$; S NEXT K END 'of main program