DECLARE FUNCTION func# (x#) DECLARE FUNCTION IClenCurtStGew% (n%, StStelle#(), Gewicht#()) DECLARE FUNCTION IClenCurt% (m%, t#(), n%, Tk#(), Ak#(), Res#) '************************************************************************ '* This program computes the value of the integral of func(x) over the * '* interval (a,b) by using the summed Clenshaw-Curtis formula * '* -------------------------------------------------------------------- * '* SAMPLE RUN: * '* (Integrate function func(x) = cos(x) - x*sin(x) from x=0 to x=PI) * '* * '* Here parameters m=4 (number of subintervals in [0, PI]), * '* n=6 (number of nodes for Clenshaw-Curtis quadrature) * '* * '* Integral = -3.141592653571974 * '* Error code = 0 * '* * '* (Exact value is -PI). * '* -------------------------------------------------------------------- * '* REF.: "Numerical Algorithms with C, by Gisela Engeln-Muellges * '* and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11]. * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '************************************************************************ 'Program Test_Clencurt DEFDBL A-H, O-Z DEFINT I-N ISIZE = 10 DIM t(ISIZE), Tk(ISIZE), Ak(ISIZE) PI = 4# * ATN(1#) m = 4: n = 6 t(0) = 0#: t(1) = PI / 4#: t(2) = PI / 2#: t(3) = 3# * PI / 4#: t(4) = PI ierror = IClenCurtStGew(n, Tk(), Ak()) ierror = IClenCurt(m, t(), n, Tk(), Ak(), Res) CLS PRINT PRINT " Integral = "; Res PRINT " Error code = "; ierror PRINT END 'of main program 'Function to integrate FUNCTION func (x) func = COS(x) - x * SIN(x) END FUNCTION FUNCTION IClenCurt (m, t(), n, Tk(), Ak(), Res) '************************************************************************ '* Computes the value of the integral of func (x) over the interval * '* (a,b) with the partition * '* t: a = t[0] < t[1] < .. < t[m] = b * '* by using the summed Clenshaw-Curtis formula. * '* This program uses precomputed [0..n] weight vectors and Chebyshev * '* nodes. * '* * '* Parameter: * '* integer m number of subintervals * '* double t(0:n) partition * '* integer n n + 1 = number of nodes, n > 1, n even * '* (n + 2 = global error order) * '* double Tk(0:n) Chebyshev nodes * '* double Ak(0:n) Weights * '* Tk and Ak must be made available before * '* calling this function by the procedure * '* ClenCurtStGew for example * '* double Res Compute integral value * '* * '* Return value : * '* 0: o.k. * '* 1: improper number of nodes * '* 2: improper number of sub intervals * '* * '* Author: Uli Eggermann, 10.3.1991 * '************************************************************************ IF n < 2 OR n MOD 2 <> 0 THEN ' n positive ? n even ? IClenCurt = 1 EXIT FUNCTION END IF IF m < 1 THEN ' partition IClenCurt = 2 EXIT FUNCTION END IF Res = 0# FOR j = 0 TO m - 1 ' loop over intervals v = .5# * (t(j + 1) - t(j)) ' half the interval size h = .5# * (t(j + 1) + t(j)) ' Interval center sum = 0# FOR k = 0 TO n ' Chebyshev loop sum = sum + Ak(k) * func(v * Tk(k) + h) NEXT k Res = Res + v * sum NEXT j IClenCurt = 0 END FUNCTION FUNCTION IClenCurtStGew (n, StStelle(), Gewicht()) '************************************************************************ '* Computes the nodes and weights of a Clenshaw-Curtis quadrature * '* formula of local error order n+3 for the reference interval [-1,1]. * '* -------------------------------------------------------------------- * '* Parameters: * '* integer n n+1 is the number of nodes and weights: * '* n > 1, n odd * '* double StStelle() computed nodes * '* double Gewicht() computed weights * '* * '* Return value : * '* 0: all ok * '* 1: n too small or odd * '* * '* Author: Uli Eggermann, 9.24.1990 * '************************************************************************ DIM i AS DOUBLE PI = 4# * ATN(1#) i = 1# IF n < 2 OR n MOD 2 <> 0 THEN ' n - Test IClenCurtStGew = 1 EXIT FUNCTION END IF m = n / 2 d = 1# * (n * n - 1) h = 2# / (n * d) f = 4# / n p = PI / n Gewicht(0) = 1# / d: Gewicht(n) = 1# / d ' left end point StStelle(0) = -1#: StStelle(n) = 1# ' right end point FOR k = 1 TO m ' Interval [-1,0] g = 0# FOR j = 1 TO m - 1 g = g + COS(2# * j * k * p) / (4# * j * j - 1#) NEXT j Gewicht(k) = h * (d + i) - f * g StStelle(k) = -COS(k * p) i = -i NEXT k StStelle(m) = 0# ' center of interval FOR k = m + 1 TO n - 1 ' Interval [0, 1] Gewicht(k) = Gewicht(n - k) StStelle(k) = -StStelle(n - k) NEXT k IClenCurtStGew = 0 END FUNCTION 'IClenCurtStGew