using QUAD
Sriram Srinivasan
I needed to use the quadrature function QUAD('F',a,b,tol).
which requires that F be a function (an m-file) which accepts x
as input and outputs F(x).
However, in my application, the function to be integrated changes
while the application is running so basically, I cannot write a single
function (m-file). Is there a way to get around this?
thanks
Sriram (ss...@kelvin.seas.virginia.edu)
Loren Shure
ss...@kelvin.seas.Virginia.EDU writes:
> I needed to use the quadrature function QUAD('F',a,b,tol).
> which requires that F be a function (an m-file) which accepts x
> as input and outputs F(x).
>
> However, in my application, the function to be integrated
changes
> while the application is running so basically, I cannot write a
single
> function (m-file). Is there a way to get around this?
I'm enclosing new versions of quad.m and quadstp.m (quad's helper) so
you can include optional arguments for your function calls. To
integrate the Jacobi elliptic function sn (first output argument of
ellipj), between 0 and pi, for modulus m=.3, use: y =
quad('ellipj',0,pi,[],[],.3)
Suppose you need your function F to change during the integration.
Then you need to write your function F to detect when to switch from
one functional form to another.
Loren Shure | Every creative act is a sudden
lo...@mathworks.com | cessation of stupidity. - Edwin Land
-----------
quad.m
-----------
function [Q,cnt] =
quad(funfcn,a,b,tol,trace,p1,p2,p3,p4,p5,p6,p7,p8,p9)
%QUAD Numerical evaluation of an integral, low order method.
% Q = QUAD('F',A,B) approximates the integral of F(X) from A
to B
% to within a relative error of 1e-3. 'F' is a string
% containing the name of the function. Function F must return
a
% vector of output values if given a vector of input values.
% Q = QUAD(F,A,B,TOL) integrates to a relative error of TOL.
% Q = Inf is returned if an excessive recursion level is
reached,
% indicating a possibly singular integral.
% Q = QUAD(F,A,B,TOL,TRACE) integrates to a relative error of
TOL and
% traces the function evaluations with a point plot of the
integrand.
%
% QUAD uses an adaptive recursive Simpson's rule.
%
% Q = QUAD('F',A,B,TOL,TRACE,P1,P2,...) allows coefficients
P1, P2, ...
% to be passed directly to function F: G = F(X,P1,P2,...).
% To use default values for TOL or TRACE, you may pass in the
empty matrix ([]).
% See also QUAD8.
% C.B. Moler, 3-22-87.
% Copyright (c) 1984-93 by The MathWorks, Inc.
% [Q,cnt] = quad(F,a,b,tol) also returns a function evaluation count.
if nargin < 4, tol = 1.e-3; trace = 0; end
if nargin < 5, trace = 0; end
c = (a + b)/2;
if isempty(tol), tol = 1.e-3; end
if isempty(trace), trace = 0; end
% Top level initialization
x = [a b c a:(b-a)/10:b];
% set up function call
args = '(x';
args1 = [];
for n = 1:nargin-5
args1 = [args1,',p',int2str(n)];
end
args1 = [args1,')'];
args = [args args1];
y = eval([funfcn,args]);
fa = y(1);
fb = y(2);
fc = y(3);
if trace
lims = [min(x) max(x) min(y) max(y)];
if any(imag(y))
lims(3) = min(min(real(y)),min(imag(y)));
lims(4) = max(max(real(y)),max(imag(y)));
end
ind = find(~finite(lims));
if ~isempty(ind)
[mind,nind] = size(ind);
lims(ind) = 1.e30*(-ones(mind,nind) .^ rem(ind,2));
end
axis(lims);
plot([c b],[fc fb],'.')
hold on
if any(imag(y))
plot([c b],imag([fc fb]),'+')
end
end
lev = 1;
% Adaptive, recursive Simpson's quadrature
if any(imag(y))
Q0 = 1e30;
else
Q0 = inf;
end
[Q,cnt] =
eval(['quadstp(funfcn,a,b,tol,lev,fa,fc,fb,Q0,trace',args1]);
cnt = cnt + 3;
if trace
hold off
axis('auto');
end
----------
quadstp.m
----------
function [Q,cnt] =
quadstp(FunFcn,a,b,tol,lev,fa,fc,fb,Q0,trace,p1,p2,p3,p4,p5,p6,p7,
8,p9)
%QUADSTP Recursive function used by QUAD.
% [Q,cnt] = quadstp(F,a,b,tol,lev,fa,fc,fb,Q0) tries to
% approximate the integral of f(x) from a to b to within a
% relative error of tol. F is a string containing the name
% of f. The remaining arguments are generated by quad or by
% the recursion. lev is the recursion level.
% fa = f(a). fc = f((a+b)/2). fb = f(b).
% Q0 is an approximate value of the integral.
% See also QUAD and QUAD8.
% C.B. Moler, 3-22-87.
% Copyright (c) 1984-93 by The MathWorks, Inc.
LEVMAX = 10;
args1 = [];
for n = 1:nargin-10
args1 = [args1,',p',int2str(n)];
end
args1 = [args1,')'];
if lev > LEVMAX
disp('Recursion level limit reached in quad. Singularity
likely.')
Q = Q0;
cnt = 0;
c = (a + b)/2;
if trace
args = '(c';
args = [args args1];
yc = eval([FunFcn,args]);
% yc = feval(FunFcn,c);
plot(c,yc,'*');
if any(imag(yc))
plot(c,imag(yc),'*');
end
end
else
% Evaluate function at midpoints of left and right half
intervals.
h = b - a;
c = (a + b)/2;
x = [a+h/4, b-h/4];
args = '(x';
args = [args args1];
f = eval([FunFcn,args]);
% f = feval(FunFcn,x);
if trace
plot(x,f,'.');
if any(imag(f))
plot(x,imag(f),'+')
end
end
cnt = 2;
% Simpson's rule for half intervals.
Q1 = h*(fa + 4*f(1) + fc)/12;
Q2 = h*(fc + 4*f(2) + fb)/12;
Q = Q1 + Q2;
% Recursively refine approximations.
if abs(Q - Q0) > tol*abs(Q)
[Q1,cnt1] =
eval(['quadstp(FunFcn,a,c,tol/2,lev+1,fa,f(1),fc,Q1,trace',args1]);
[Q2,cnt2] =
eval(['quadstp(FunFcn,c,b,tol/2,lev+1,fc,f(2),fb,Q2,trace',args1]);
% [Q1,cnt1] =
quadstp(FunFcn,a,c,tol/2,lev+1,fa,f(1),fc,Q1,trace);
% [Q2,cnt2] =
quadstp(FunFcn,c,b,tol/2,lev+1,fc,f(2),fb,Q2,trace);
Q = Q1 + Q2;
cnt = cnt + cnt1 + cnt2;
end
end