!************************************************************************ !* 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.14159265357197 * !* Error code = 0 * !* * !* (Exact value is -PI). * !* -------------------------------------------------------------------- * !* REF.: "Numerical Algorithms with C, By Gisela Engeln-Muellges * !* and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11]. * !* * !* F90 Release By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !************************************************************************ Program Test_Clencurt real*8, parameter :: ONE = 1.d0 real*8, parameter :: ZERO = 0.d0 integer, parameter :: SIZE = 10 integer error, m, n, ClenCurt real*8 t(0:SIZE), Tk(0:SIZE), Ak(0:SIZE) real*8 PI, Res PI = 4.d0*datan(ONE) m=4; n=6 t(0)=ZERO; t(1)=PI/4.d0; t(2)=PI/2.d0; t(3)=3.d0*PI/4.d0; t(4)=PI call ClenCurtStGew(n,Tk,Ak) error = ClenCurt(m,t,n,Tk,Ak,Res) print *,' ' print *,' Integral = ', Res print *,' Error code = ', error print *,' ' stop ' End of program.' END !Function to integrate real*8 Function func(x) real*8 x func = dcos(x)-x*dsin(x) return End Integer Function ClenCurtStGew(n, StStelle, Gewicht) integer n; real*8 StStelle(0:n), Gewicht(0:n) !************************************************************************ !* 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 * !* real*8 StStelle[] computed nodes * !* real*8 Gewicht[] computed weights * !* * !* Return value : * !* 0: all ok * !* 1: n too small or odd * !* * !* Author: Uli Eggermann, 9.24.1990 * !************************************************************************ integer k, j, m real*8 p, f, g, h, i, d, PI PI = 4.d0*datan(1.d0) i=1.d0 if (n < 2.or.Mod(n,2).ne.0) then ! n - Test ClenCurtStGew = 1 return end if m = n / 2 d = 1.d0*(n * n - 1) h = 2.d0 / (n * d) f = 4.d0 / n p = PI/n Gewicht(0) = 1.d0/d; Gewicht(n) = 1.d0/d ! left end point StStelle(0) = -1.d0; StStelle(n) = 1.d0 ! right end point do k = 1, m ! Interval [-1,0] g = 0.d0 do j = 1, m-1 g = g + dcos(2.d0 * j * k * p) / (4.d0 * j * j - 1.d0) end do Gewicht(k) = h * (d + i) - f * g StStelle(k) = - dcos(k * p) i = -i end do StStelle(m) = 0.d0 ! center of interval do k = m+1, n-1 ! Interval [0, 1] Gewicht(k) = Gewicht(n - k) StStelle(k) = -StStelle(n - k) end do ClenCurtStGew = 0 return End !ClenCurtStGew Integer Function ClenCurt(m, t, n, Tk, Ak, Res) integer m,n; real*8 t(0:m),Tk(0:n),Ak(0:n),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 * !* real*8 t(0:n) partition * !* integer n n + 1 = number of nodes, n > 1, n even * !* (n + 2 = global error order) * !* real*8 Tk(0:n) Chebyshev nodes * !* real*8 Ak(0:n) Weights * !* Tk and Ak must be made available before * !* calling this function by the procedure * !* ClenCurtStGew for example * !* real*8 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 * !************************************************************************ integer j, k real*8 func, v, h, sum if (n < 2.or.Mod(n,2).ne.0) then ! n positive ? n even ? ClenCurt = 1 return end if if (m < 1) then ! partition ClenCurt = 2 return end if Res = 0.d0 do j = 0, m-1 ! loop over intervals v = 0.5d0 * (t(j+1) - t(j)) ! half the interval size h = 0.5d0 * (t(j+1) + t(j)) ! Interval center sum = 0.d0 do k = 0, n ! Chebyshev loop sum = sum + Ak(k)*func(v*Tk(k) + h) end do Res = Res + v * sum end do ClenCurt = 0 return End ! ------------------------- END clencurt.f90 --------------------------