Home > wafo > misc > ccquad.m

ccquad

PURPOSE ^

Numerical integration using a Clenshaw-Curtis quadrature.

SYNOPSIS ^

[int,tol] = ccquad(fun,a,b,n,trace)

DESCRIPTION ^

  CCQUAD Numerical integration using a Clenshaw-Curtis quadrature.
 
  CALL:     [int, tol] = ccquad(Fun,a,b,n) 
 
     int = evaluated integral
        tol = estimate of the absolute error (usually conservative).
     Fun = string containing the name of the function or a directly 
              given expression enclosed in parenthesis. 
        a,b = integration limits
          n = number of base points (abscissas). Default n=10
 
   The integral is exact for polynomials of degree n or less.
   Usually this routine gives accurate answers for smooth functions.  
 
 Example:% Integration of exp(x) from 0 to 1:
 
    a=0; b=1;
    ccquad('exp(x)',a,b) % Should be exp(1)-1
 
  See also  gaussq, qrule2d

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001  function [int,tol] = ccquad(fun,a,b,n,trace)
002 % CCQUAD Numerical integration using a Clenshaw-Curtis quadrature.
003 %
004 % CALL:     [int, tol] = ccquad(Fun,a,b,n) 
005 %
006 %    int = evaluated integral
007 %       tol = estimate of the absolute error (usually conservative).
008 %    Fun = string containing the name of the function or a directly 
009 %             given expression enclosed in parenthesis. 
010 %       a,b = integration limits
011 %         n = number of base points (abscissas). Default n=10
012 %
013 %  The integral is exact for polynomials of degree n or less.
014 %  Usually this routine gives accurate answers for smooth functions.  
015 %
016 %Example:% Integration of exp(x) from 0 to 1:
017 %
018 %   a=0; b=1;
019 %   ccquad('exp(x)',a,b) % Should be exp(1)-1
020 %
021 % See also  gaussq, qrule2d
022 
023 % References:
024 % [1] Goodwin, E.T. (1961),
025 % "Modern Computing Methods",
026 % 2nd edition, New yourk: Philosophical Library, pp. 78--79
027 %
028 % [2] Clenshaw, C.W. and Curtis, A.R. (1960),
029 % Numerische Matematik, Vol. 2, pp. 197--205
030 
031 % tested on: matlab 5.3
032 % history:
033 % revised pab 22Nov2004
034 % Added the possibility of using a function handle.
035 % by    Per A. Brodtkorb 26.07.1999
036 
037 if nargin<5|isempty(trace),
038    trace=0;
039 end
040 if nargin<4|isempty(n),
041    n = 10;
042 else
043    % make sure n is even
044    n = 2*ceil(n/2);
045 end;
046 
047 if (isa(fun,'char') &  any(fun=='(')), %  & any(fun=='x'),
048   exec_string=['f=',fun ';']; %the call function is already setup
049 else
050   if isa(fun,'function_handle')
051     fun = func2str(fun);
052   end
053   %setup string to call the function
054   exec_string=['f=feval(fun,x);'];
055 end
056 
057 
058 
059 
060 s = (0:n)';
061 s2 =(0:2:n)';
062 
063 x = cos(pi*s/n)*(b-a)/2+(b+a)/2;
064 eval(exec_string);
065 if trace==1,
066   plot(x,f,'+')
067 end
068 
069 % using a Gauss-Lobatto variant, i.e., first and last
070 % term f(a) and f(b) is multiplied with 0.5
071 f(1) = f(1)/2;
072 f(n+1) = f(n+1)/2;
073 
074 if 1,%fft for faster calculations
075  c=real(fft(f(1:n)));
076  c=2/n*(c(1:n/2+1)+f(n+1)*cos(pi*s2));
077 else %old call: slow for large n
078   c = 2/n * cos(s2*s'*pi/n) * f;
079 end
080 c(1) = c(1)/2;
081 c(n/2+1) = c(n/2+1)/2;
082 
083 int = (a-b)*sum(c./(s2-1)./(s2+1) );
084 tol = abs(c(n/2+1));
085 
086

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003