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
Discrete Fourier transform. | |
Construct a string from a function handle. | |
True if object is a given class. | |
Linear plot. |
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
Comments or corrections to the WAFO group