www.gusucode.com > funfun工具箱matlab源码程序 > funfun/quadl.m
function [Q,fcnt] = quadl(funfcn,a,b,tol,trace,varargin) %QUADL Numerically evaluate integral, adaptive Lobatto quadrature. % Q = QUADL(FUN,A,B) tries to approximate the integral of scalar-valued % function FUN from A to B to within an error of 1.e-6 using high order % recursive adaptive quadrature. FUN is a function handle. The function % Y=FUN(X) should accept a vector argument X and return a vector result % Y, the integrand evaluated at each element of X. % % Q = QUADL(FUN,A,B,TOL) uses an absolute error tolerance of TOL % instead of the default, which is 1.e-6. Larger values of TOL % result in fewer function evaluations and faster computation, % but less accurate results. % % Q = QUADL(FUN,A,B,TOL,TRACE) with non-zero TRACE shows the values % of [fcnt a b-a Q] during the recursion. Use [] as a placeholder to % obtain the default value of TOL. % % [Q,FCNT] = QUADL(...) returns the number of function evaluations. % % Use array operators .*, ./ and .^ in the definition of FUN % so that it can be evaluated with a vector argument. % % QUADL will be removed in a future release. Use INTEGRAL instead. % % Example: % Q = quadl(@myfun,0,2); % where the file myfun.m defines the function: % %-------------------% % function y = myfun(x) % y = 1./(x.^3-2*x-5); % %-------------------% % % or, use a parameter for the constant: % Q = quadl(@(x)myfun2(x,5),0,2); % where the file myfun2.m defines the function: % %----------------------% % function y = myfun2(x,c) % y = 1./(x.^3-2*x-c); % %----------------------% % % Class support for inputs A, B, and the output of FUN: % float: double, single % % See also INTEGRAL, INTEGRAL2, INTEGRAL3, QUADGK, QUAD2D, TRAPZ, ... % FUNCTION_HANDLE. % Based on "adaptlob" by Walter Gautschi. % Ref: W. Gander and W. Gautschi, "Adaptive Quadrature Revisited", 1998. % http://www.inf.ethz.ch/personal/gander % Copyright 1984-2013 The MathWorks, Inc. f = fcnchk(funfcn); if nargin < 4 || isempty(tol), tol = 1.e-6; end; if nargin < 5 || isempty(trace), trace = 0; end; % Initialize with 13 function evaluations. c = (a + b)/2; h = (b - a)/2; s = [.942882415695480 sqrt(2/3) .641853342345781 1/sqrt(5) .236383199662150]; x = [a c-h*s c c+h*fliplr(s) b]; y = feval(f,x,varargin{:}); y = y(:).'; fcnt = 13; if numel(y) ~= fcnt error(message('MATLAB:quadl:funNotVectorized')); end % Fudge endpoints to avoid infinities. if ~isfinite(y(1)) y(1) = feval(f,a+eps(superiorfloat(a,b))*(b-a),varargin{:}); fcnt = fcnt+1; end if ~isfinite(y(13)) y(13) = feval(f,b-eps(superiorfloat(a,b))*(b-a),varargin{:}); fcnt = fcnt+1; end % Increase tolerance if refinement appears to be effective. Q1 = (h/6)*[1 5 5 1]*y(1:4:13).'; Q2 = (h/1470)*[77 432 625 672 625 432 77]*y(1:2:13).'; s = [.0158271919734802 .094273840218850 .155071987336585 ... .188821573960182 .199773405226859 .224926465333340]; w = [s .242611071901408 fliplr(s)]; Q0 = h*w*y.'; r = abs(Q2-Q0)/abs(Q1-Q0+realmin(class(Q0))); if r > 0 && r < 1 tol = tol/r; end; % Call the recursive core integrator. hmin = eps(b-a)/1024; [Q,fcnt,warn] = quadlstep(f,a,b,y(1),y(13),tol,trace,fcnt,hmin,varargin{:}); switch warn case 1 warning(message('MATLAB:quadl:MinStepSize')) case 2 warning(message('MATLAB:quadl:MaxFcnCount')) case 3 warning(message('MATLAB:quadl:ImproperFcnValue')) otherwise % No warning. end % ------------------------------------------------------------------------ function [Q,fcnt,warn] = quadlstep(f,a,b,fa,fb,tol,trace,fcnt,hmin,varargin) %QUADLSTEP Recursive core routine for function QUADL. maxfcnt = 10000; % Evaluate integrand five times in interior of subinterval [a,b]. c = (a + b)/2; h = (b - a)/2; if abs(h) < hmin || c == a || c == b % Minimum step size reached; singularity possible. Q = h*(fa+fb); warn = 1; return end alpha = sqrt(2/3); beta = 1/sqrt(5); x = [c-alpha*h c-beta*h c c+beta*h c+alpha*h]; y = feval(f,x,varargin{:}); fcnt = fcnt + 5; if fcnt > maxfcnt % Maximum function count exceeded; singularity likely. Q = h*(fa+fb); warn = 2; return end x = [a x b]; y = [fa y(:).' fb]; % Four point Lobatto quadrature. Q1 = (h/6)*[1 5 5 1]*y(1:2:7).'; % Seven point Kronrod refinement. Q2 = (h/1470)*[77 432 625 672 625 432 77]*y.'; Q = Q2; if ~isfinite(Q) % Infinite or Not-a-Number function value encountered. warn = 3; return end if trace fprintf('%8.0f %16.10f %18.8e %16.10f\n',fcnt,a,h,Q); end % Check accuracy of integral over this subinterval. if abs(Q1 - Q2) <= tol warn = 0; return % Subdivide into six subintervals. else Q = 0; warn = 0; for k = 1:6 [Qk,fcnt,wk] = quadlstep(f,x(k),x(k+1),y(k),y(k+1), ... tol,trace,fcnt,hmin,varargin{:}); Q = Q + Qk; warn = max(warn,wk); end end