积分与微分方程
function = quad8(funfcn,a,b,tol,trace,varargin)%数值积分
%z=quad8('Fun',A,B,Tol,trace,p1,p2,L)
% 其中:"Fun"-表示被积函数的M函数名.
% A,B-上﹑下限.
% Tol-精度,缺省值为1e-3.
% trace-非零时显示计算过程,缺省值为0.
% p1,p2,L-参变量,一般缺省.
%例如求 exp(-x^2)在[-1，1]积分
% 先写M函数quadeg3fun.m
% function y=quadeg3fun(x)
% y=exp(-x.^2);
% 用
% z=quad8('quadeg3fun',-1,1)
%
%QUAD8 Numerically evaluate integral, higher order method.
% Q = QUAD8('F',A,B) approximates the integral of F(X) from A to B to
% within a relative error of 1e-3 using an adaptive recursive Newton
% Cotes 8 panel rule. 'F' is a string containing the name of the
% function. The function must return a vector of output values if
% given a vector of input values. Q = Inf is returned if an excessive
% recursion level is reached, indicating a possibly singular integral.
%
% Q = QUAD8('F',A,B,TOL) integrates to a relative error of TOL. Use
% a two element tolerance, TOL = , to specify a
% combination of relative and absolute error.
%
% Q = QUAD8('F',A,B,TOL,TRACE) integrates to a relative error of TOL and
% for non-zero TRACE traces the function evaluations with a point plot
% of the integrand.
%
% Q = QUAD8('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 QUAD, DBLQUAD.
% Cleve Moler, 5-08-88.
% Copyright (c) 1984-98 by The MathWorks, Inc.
% $Revision: 5.16 $ $Date: 1997/11/21 23:31:14 $
% = quad8(F,a,b,tol) also returns a function evaluation count.
if nargin < 4, tol = ; trace = 0; end
if nargin < 5, trace = 0; end
if isempty(tol), tol = ; end
if length(tol)==1, tol = ; end
if isempty(trace), trace = 0; end
% QUAD8 usually does better than the default 1e-3.
h = b - a;
% Top level initialization, Newton-Cotes weights
w = /14175;
x = a + (0:8)*(b-a)/8;
y = feval(funfcn,x,varargin{:});
yllow = min();
ylhi = max();
lims = ;
ind = find(~isfinite(lims));
if ~isempty(ind)
= size(ind);
lims(ind) = 1.e30*(-ones(mind,nind) .^ rem(ind,2));
end
if trace
axis(lims);
% doesn't take care of complex case
plot(,,'.'), hold on
if ~isreal(y)
plot(,,'+')
end
drawnow
end
lev = 1;
% Adaptive, recursive Newton-Cotes 8 panel quadrature
if ~isreal(y)
Q0 = 1e30;
else
Q0 = inf;
end
recur_lev_excess = 0;
= quad8stp(funfcn,a,b,tol,lev,w,x,y,Q0,trace,recur_lev_excess,varargin{:});
cnt = cnt + 9;
if trace
hold off
axis('auto');
end
if (recur_lev_excess > 1)
warning(sprintf('Recursion level limit reached %d times.', ...
recur_lev_excess ))
end
%------------------------------------------------------------------
function = quad8stp(FunFcn,a,b,tol,lev,w,x0,f0,Q0,trace,recur_lev_excess,varargin)
%QUAD8STP Recursive function used by QUAD8.
% = quad8stp(F,a,b,tol,lev,w,f,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 quad8 or by the recursion. lev is the recursion level.
% w is the weights in the 8 panel Newton Cotes formula.
% x0 is a vector of 9 equally spaced abscissa is the interval.
% f0 is a vector of the 9 function values at x.
% Q0 is an approximate value of the integral.
% See also QUAD8 and QUAD.
% Cleve Moler, 5-08-88.
LEVMAX = 10;
% Evaluate function at midpoints of left and right half intervals.
x = zeros(1,17);
f = zeros(1,17);
x(1:2:17) = x0;
f(1:2:17) = f0;
x(2:2:16) = (x0(1:8) + x0(2:9))/2;
f(2:2:16) = feval(FunFcn,x(2:2:16),varargin{:});
if trace
plot(x(2:2:16),f(2:2:16),'.');
if ~isreal(f)
plot(x(2:2:16),imag(f(2:2:16)),'+');
end
drawnow
end
cnt = 8;
% Integrate over half intervals.
h = (b-a)/16;
Q1 = h*w*f(1:9).';
Q2 = h*w*f(9:17).';
Q = Q1 + Q2;
% Recursively refine approximations.
if (abs(Q - Q0) > tol(1)*abs(Q)+tol(2)) & (lev <= LEVMAX)
c = (a+b)/2;
= quad8stp(FunFcn,a,c,tol/2,lev+1,w,x(1:9),f(1:9),Q1,trace,recur_lev_excess,varargin{:});
= quad8stp(FunFcn,c,b,tol/2,lev+1,w,x(9:17),f(9:17),Q2,trace,recur_lev_excess,varargin{:});
Q = Q1 + Q2;
cnt = cnt + cnt1 + cnt2;
elseif (lev > LEVMAX)
if ~recur_lev_excess
warning('Recursion level limit reached in quad8. Singularity likely.')
recur_lev_excess = 1;
else
recur_lev_excess = recur_lev_excess + 1;
end
end
不错不错。。。。。。。。。。。。。。。。。。。。。。。。。。。。。
页:
[1]