'********************************************************** '* This program calculates the Fourier coefficients * '* of a periodic function F(x) by using the subroutine * '* 1000 (AnalyticFourier). * '* ------------------------------------------------------ * '* SAMPLE RUN: * '* * '* Example: F(x) equals 0 from -pi to 0 and sin(x) * '* from 0 to pi. * '* * '* Calculate the Fourier coefficients of a periodic * '* function F(x): * '* * '* Lowest harmonic: 0 * '* Highest harmonic: 10 * '* * '* a0 = 0.318309886 * '* b0 = 0.000000000 * '* * '* a1 = 0.000000000 * '* b1 = 0.500000000 * '* * '* a2 = -0.212206591 * '* b2 = 0.000000000 * '* ----------------- * '* a9 = -0.000000000 * '* b9 = -0.000000000 * '* * '* a10 = -0.006430503 * '* b10 = -0.000000000 * '* * '* Basic Release By J-P Moreau, Paris * '* (www.jpmoreau.fr) * '********************************************************** 'PROGRAM Analytic_Fourier DEFDBL A-H, O-Z DEFINT I-N 'integer icount,ih1,ih2 'REAL*8 t1,t2, a,b, PI MAXITER = 15 DIM T(MAXITER, MAXITER) CLS PI = 3.1415926535# PRINT PRINT " Calculate the Fourier coefficients of a periodic function F(x):" PRINT t1 = -PI: t2 = PI INPUT " Lowest harmonic : ", ih1 INPUT " Highest harmonic: ", ih2 icount = 1 FOR ii = ih1 TO ih2 PRINT GOSUB 1000 'call AnalyticFourier(t1,t2,i,a,b) PRINT " a"; ii; " = "; : PRINT USING "##.#########"; a PRINT " b"; ii; " = "; : PRINT USING "##.#########"; b IF (icount MOD 5) = 0 THEN icount = 0 PRINT INPUT " Press any key to continue ", R$ END IF icount = icount + 1 NEXT ii PRINT PRINT END ' analytical function(s) to study: ' ------------------------------- ' Note: Exact Fourier coefficients for this periodic function are: ' an=0 if n is even, -1/2n if n=4p+1, 1/2n if n=4p+3 ' bn=0 if n=4p, -1/2n if n=4p+1, 1/n if n=4p+2, -1/2n if n=4p+3 ' ---------------------------------------------------------------- 500 'Function F(x) (example 1) IF (x < -PI) THEN F = 0# ELSEIF (x < -PI / 2#) THEN F = PI / 4# ELSEIF (x <= PI) THEN F = -PI / 4# ELSE F = 0# END IF RETURN ' Note: Exact Fourier coefficients for this periodic function are: ' a0=1/pi, an = 0 if n is odd, an = -2/(pi*(n2-1)) if n is even; ' b1=1/2, bn = 0 for n>=2. '----------------------------------------------------------------- 501 'Function F1(x) (example 2) IF (x < 0#) THEN F1 = 0# ELSEIF (x < PI) THEN F1 = SIN(x) ELSE F1 = 0# END IF RETURN ' Function to integrate by Romberg method 600 'Function FUNC(om,x,flag) IF (iflag = 0) THEN GOSUB 501 FUNC = (F1 * COS(om * x)) 'for an ELSE GOSUB 501 FUNC = (F1 * SIN(om * x)) 'for bn END IF RETURN '**************************************************************** '* Calculate the Fourier harmonic #n of a periodic function F(x)* '* analytically defined. * '* ------------------------------------------------------------ * '* Inputs: * '* t1 : -period/2 * '* t2 : period/2 * '* n : order of harmonic * '* * '* Outputs: * '* a : coefficient an of the Fourier series. * '* b: : coefficient bn of the Fourier series. * '* ------------------------------------------------------------ * '* Reference: "Mathematiques en Turbo-Pascal by Marc Ducamp and * '* Alain Reverchon, 1. Analyse, Editions Eyrolles, * '* Paris, 1991" [3]. * '* * '* Basic version by J-P Moreau, Paris. * '**************************************************************** ' Note: When a periodic function f(x), of period T, can be developped ' in a Fourier series, this one is unique: ' n=inf. ' F(x) = a0 + Sum ( an cos(2 n pi/T x) + bn sin(2 n pi/T x) ' n=1 ' (b0 is always zero).' ' ' The coefficients an, bn are given by: ' T/2 ' a0 = 1/T Sum (f(x) dx) ' -T/2 ' T/2 ' an = 2/T Sum (f(x) cos(2npi/T x) dx) ' -T/2 ' T/2 ' an = 2/T Sum (f(x) sin(2npi/T x) dx) ' -T/2 ' Here, the integrals are calculated by the Romberg method. ' -------------------------------------------------------------------- 1000 'Subroutine AnalyticFourier(t1,t2,n,a,b) 'integer iflag T = ABS(t2 - t1): n = ii om = 2# * PI * n / T prec = 1E-10: itermin = 5: itermax = 15 iflag = 0 'calculate an aa = t1: bb = t2 GOSUB 2000 'xintegral=RombergIntegral(t1,t2,prec,obtprec,niter,5,itermax) a = xintegral * 2# / T iflag = 1 'calculate bn GOSUB 2000 'xintegral=RombergIntegral(t1,t2,prec,obtprec,niter,5,itermax) b = xintegral * 2# / T IF n = 0 THEN a = a / 2# RETURN '******************************************************* '* Integral of a function FUNC(X) by Romberg's method * '* --------------------------------------------------- * '* INPUTS: * '* aa begin value of x variable * '* bb end value of x variable * '* prec desired precision * '* itermin minimum number of iterations * '* itermax maximum number of iterations * '* * '* OUTPUTS: * '* obtprec obtained precision for integral * '* niter number of iterations done * '* xintegral the integral of FUNC(X) from a to b * '* * '******************************************************* 2000 'Function RombergIntegral IF itermax > MAXITER THEN itermax = MAXITER x = aa: GOSUB 600: R = FUNC x = bb: GOSUB 600 ta = (R + FUNC) / 2# niter = 0 pas = bb - aa T(0, 0) = ta * pas 2100 niter = niter + 1 pas = pas / 2# s = ta FOR i = 1 TO 2 ^ niter - 1 x = aa + pas * i: GOSUB 600 s = s + FUNC NEXT i T(0, niter) = s * pas R = 1# FOR i = 1 TO niter R = R * 4# j = niter - i T(i, j) = (R * T(i - 1, j + 1) - T(i - 1, j)) / (R - 1#) NEXT i obtprec = ABS(T(niter, 0) - T(niter - 1, 0)) IF niter > itermax THEN GOTO 2200 IF niter < itermin THEN GOTO 2100 IF obtprec > prec THEN GOTO 2100 2200 xintegral = T(niter, 0) RETURN ' end of file analfour.bas