标题: 积分与微分方程 [打印本页] 作者: 李秋锐 时间: 2015-2-4 20:17 标题: 积分与微分方程 function [Q,cnt] = 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 = [rel_tol abs_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 $
% [Q,cnt] = quad8(F,a,b,tol) also returns a function evaluation count.
if nargin < 4, tol = [1.e-3 0]; trace = 0; end
if nargin < 5, trace = 0; end
if isempty(tol), tol = [1.e-3 0]; end
if length(tol)==1, tol = [tol 0]; 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 = [3956 23552 -3712 41984 -18160 41984 -3712 23552 3956]/14175;
x = a + (0:8)*(b-a)/8;
y = feval(funfcn,x,varargin{:});
yllow = min([min(real(y)) min(imag(y))]);
ylhi = max([max(real(y)) max(imag(y))]);
lims = [min(x) max(x) yllow ylhi];
ind = find(~isfinite(lims));
if ~isempty(ind)
[mind,nind] = 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([a b],[real(y(1)) real(y(9))],'.'), hold on
if ~isreal(y)
plot([a b],[imag(y(1)) imag(y(9))],'+')
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;
[Q,cnt,recur_lev_excess] = 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 [Q,cnt,recur_lev_excess] = quad8stp(FunFcn,a,b,tol,lev,w,x0,f0,Q0,trace,recur_lev_excess,varargin)
%QUAD8STP Recursive function used by QUAD8.
% [Q,cnt] = 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;
[Q1,cnt1,recur_lev_excess] = quad8stp(FunFcn,a,c,tol/2,lev+1,w,x(1:9),f(1:9),Q1,trace,recur_lev_excess,varargin{:});
[Q2,cnt2,recur_lev_excess] = 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