Home > wafo > misc > ccquad.m

## PURPOSE

Numerical integration using a Clenshaw-Curtis quadrature.

## DESCRIPTION

```  CCQUAD Numerical integration using a Clenshaw-Curtis quadrature.

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;

## CROSS-REFERENCE INFORMATION

This function calls:
 fft Discrete Fourier transform. func2str Construct a string from a function handle. isa True if object is a given class. plot Linear plot.
This function is called by:

## SOURCE CODE

```001  function [int,tol] = ccquad(fun,a,b,n,trace)
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 %
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