5 views

### Sriram Srinivasan

Mar 23, 1993, 10:07:50 PM3/23/93
to
Hi,

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

Mar 31, 1993, 11:41:47 AM3/31/93
to
In article <C4DI1...@murdoch.acc.Virginia.EDU> Sriram Srinivasan,

ss...@kelvin.seas.Virginia.EDU writes:
> 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?

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 =

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

-----------
-----------
function [Q,cnt] =
%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.
%
%
% 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 ([]).

% 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;

if any(imag(y))
Q0 = 1e30;
else
Q0 = inf;
end
[Q,cnt] =
cnt = cnt + 3;
if trace
hold off
axis('auto');
end

----------
----------
function [Q,cnt] =
8,p9)
% [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.

% 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] =
[Q2,cnt2] =