/*********************************************************************** * 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.141593 * * Error code = 0 * * * * (Exact value is -PI). * * -------------------------------------------------------------------- * * REF.: "Numerical Algorithms with C, By Gisela Engeln-Muellges * * and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11]. * * * * Visual C++ Release By J-P Moreau, Paris. * * (www.jpmoreau.fr) * ***********************************************************************/ #include #include #define ONE 1.0 #define PI 4.0*atan(1.0) #define REAL double #define ZERO 0.0 #define SIZE 10 //Function to integrate REAL func(REAL x) { return (cos(x)-x*sin(x)); } int ClenCurtStGew ( int n, REAL* StStelle, REAL* 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: * * int n n+1 is the number of nodes and weights: * * n > 1, n odd * * REAL StStelle [] computed nodes * * REAL Gewicht [] computed weights * * * * Return value : * * 0: all ok * * 1: n too small or odd * * * * Author: Uli Eggermann, 9.24.1990 * ***********************************************************************/ { int k, j, m; REAL p, f, g, h, i = ONE, d; if (n < 2 || n % 2 != 0) return 1; /* n - Test */ m = n / 2; /* m = n/2 */ d = (REAL)(n * n - 1); /* d = n*n+1 */ h = 2. / (n * d); /* h = 2/(n(n*n+1)) */ f = 4. / n; /* f = 4/n */ p = f * atan (1.); /* p = pi/n */ Gewicht [0] = Gewicht [n] = 1 / d; /* left end point */ StStelle [0] = -ONE; StStelle [n] = ONE; /* right end point */ for (k = 1; k <= m; k++, i = -i) /* Interval (-1,0] */ { for (g = ZERO, j = 1; j < m; j++) g += cos(2 * j * k * p) / (4 * j * j - 1); Gewicht [k] = h * (d + i) - f * g; StStelle [k] = - cos(k * p); } StStelle [m] = ZERO; /* center of interval */ for (k = m+1; k < n; k++) /* Interval (0, 1) */ { Gewicht [k] = Gewicht [n - k]; StStelle [k] = - StStelle [n - k]; } return 0; } // ClenCurtStGew int ClenCurt ( int m, REAL* t, int n, REAL* Tk, REAL* Ak, REAL* Resultat ) /*********************************************************************** * 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: * * int m number of subintervals * * REAL t [] partition * * int n n + 1 = number of nodes, n > 1, n even * * (n + 2 = global error order) * * REAL Tk [] Chebyshev nodes * * REAL Ak [] Weights * * Tk and Ak must be made available before * * calling this function by the procedure * * ClenCurtStGew for example * * REAL *Resultat 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 * ***********************************************************************/ { int j, k; REAL v, h, sum; if (n < 2 || n % 2 != 0) return 1; /* n positive ? n even ? */ if (m < 1) return 2; /* partition */ for (*Resultat = ZERO, j = 0; j < m; j++) /* loop over intervals */ { v = 0.5 * (t [j+1] - t [j]); /* half the interval size */ h = 0.5 * (t [j+1] + t [j]); /* Interval center */ for (sum = ZERO, k = 0; k <= n; k++) /* Chebyshev loop */ sum += Ak [k] * func (v * Tk [k] + h); *Resultat += v * sum; } return 0; } void main() { int error, m=4, n=6; REAL t[SIZE], Tk[SIZE], Ak[SIZE]; REAL Result; t[0]=ZERO; t[1]=PI/4.0; t[2]=PI/2.0; t[3]=3.0*PI/4.0; t[4]=PI; ClenCurtStGew(n,Tk,Ak); error = ClenCurt(m,t,n,Tk,Ak,&Result); printf("\n Integral = %f\n", Result); printf(" Error code = %d\n\n", error); } /* ------------------------ END clencurt.cpp ------------------------ */