gusucode.com > symbolic工具箱matlab源码程序 > symbolic/@sym/vpaintegral.m
function y = vpaintegral(varargin)
%VPAINTEGRAL Numerical integration of symbolic expressions.
% Q = VPAINTEGRAL(FUN,XMIN,XMAX) approximates the integral of function
% FUN from XMIN to XMAX by variable precision arithmetic using global
% adaptive quadrature and default error tolerances.
%
% FUN must be a symbolic expression or an array of symbolic expressions.
% If FUN contains more than one symbolic variable, the input syntax
% Q = VPAINTEGRAL(FUN,X,XMIN,XMAX) should be used to indicate the
% variable X that serves as the integration variable.
% XMIN and XMAX can be symbolic numbers or double numbers or -Inf or
% Inf. If both are finite, they can be complex. If at least one is
% complex, VPAINTEGRAL approximates the path integral from XMIN to
% XMAX over a straight line path.
%
% The available options are
%
% 'AbsTol', absolute error tolerance
% 'RelTol', relative error tolerance
%
% VPAINTEGRAL attempts to satisfy |Q-I|<=max(AbsTol,RelTol*|Q|),
% where I denotes the exact value of the integral. Usually RelTol
% determines the accuracy of the integration. However, if |Q| is
% sufficiently small, AbsTol determines the accuracy of the
% integration, instead. The default value of AbsTol is 1e-10, and
% the default value of RelTol is 1e-6. The value of AbsTol must
% be nonnegative, the value of RelTol must be positive.
%
% 'Waypoints', vector of integration waypoints
%
% If FUN has discontinuities or singularities in the interval of
% integration, the locations should be supplied as a 'Waypoints'
% vector. If XMIN, XMAX, or any entry of the waypoints vector is
% complex, the integration is performed over a sequence of straight
% line paths in the complex plane, from XMIN to the first waypoint,
% from the first waypoint to the second, and so forth, and finally
% from the last waypoint to XMAX.
%
% 'MaxFunctionCalls', maximal number of evaluations of FUN
%
% The numerical integration gives up and throws an error if the
% number of internal FUN evaluations exceeds this value. The value
% of this option must be a positive integer or Inf. The default
% value is 10^5.
%
% If any of FUN, XMIN, XMAX, or the 'Waypoints' vector contain a
% symbolic variable other than the integration variable, then a
% symbolic call of VPAINTEGRAL is returned. When several symbolic
% variables are present, it is recommended to use the input syntax
% VPAINTEGRAL(FUN,X,XMIN,XMAX) in which the integration variable X
% is explicitly specified. Higher dimensional integrals can be
% computed by nested calls such as
% vpaintegral(vpaintegral(FUN(X,Y),Y,YMIN(X),YMAX(X)),X,XMIN,XMAX).
%
% Examples:
% % Integrate f(x) = exp(-x^2)*log(x)^2 from 0 to infinity.
% syms x;
% f = exp(-x^2)*log(x)^2;
% Q = vpaintegral(f,0,Inf)
%
% Q =
%
% 1.94752
%
% % Specify tolerances.
%
% Q = vpaintegral(sin(x)/x^3,1,Inf,'AbsTol',1e-16,'RelTol',1e-10)
%
% Q =
%
% 0.3785300171
%
% % Integrate f(x) = 1/(2*x-1) in the complex plane over the
% % closed triangular path from 0 to 1+1i to 1-1i to 0. The
% % result changes its sign when the orientation of the path
% % is reverted.
% Q = vpaintegral(1/(2*x-1),0,0,'Waypoints',[1+1i,1-1i])
%
% Q =
%
% 1.11022e-16 - 3.14159i
%
% Q = vpaintegral(1/(2*x-1),0,0,'Waypoints',[1-1i,1+1i])
%
% Q =
%
% 1.11022e-16 + 3.14159i
%
% % Integrate the vector-valued function sin((1:5)*x) from 0 to pi.
% Q = vpaintegral(sin((1:5)*x),0,pi)
%
% Q =
%
% [ 2.0, 0, 0.666667, -1.08996e-16, 0.4]
%
% % Use a parameter c in the integrand. The integration variable x
% % is specified explicitly. The result cannot be computed
% % numerically because of the symbolic parameter.
% syms f(x,c)
% f(x,c) = 1/(x^3-2*x-c);
% Q = vpaintegral(f,x,0,2)
%
% Q(c) =
%
% vpaintegral(-1/(- x^3 + 2*x + c), x, 0, 2)
%
% % Multiple integration over the triangular region 0<=x<=1,|y|<x.
%
% syms x y;
% Q = vpaintegral(vpaintegral(sin(x-y)/(x-y),y,[-x,x]),x,[0,1])
%
% Q =
%
% 0.89734
%
% See also INT, INTEGRAL, DSOLVE, VPA.
% Copyright 2015 The MathWorks, Inc.
narginchk(2,Inf);
k = find(cellfun(@ischar, varargin), 1);
if isempty(k)
[fun, x, a, b] = parseFunctionVariableBounds(varargin{:});
elseif k >= 3
[fun, x, a, b] = parseFunctionVariableBounds(varargin{1:min(k-1,4)});
else % k = 2
error(message('symbolic:sym:vpaintegral:UnexpectedCharInput'));
end
if ~sym.isVariable(x)
error(message('symbolic:sym:int:InvalidIntegrationVariable', char(x)));
end
p = inputParser;
p.addParameter('RelTol', 1.0e-06, @(x) ~ischar(x) && isscalar(x) && isnumeric(double(x)) && logical(x>0));
p.addParameter('AbsTol', 1.0e-10, @(x) ~ischar(x) && isscalar(x) && isnumeric(double(x)) && logical(x>=0));
p.addParameter('Waypoints', sym([]), @(x) ~ischar(x) && (isempty(x) || isvector(sym(x))));
p.addParameter('MaxFunctionCalls', sym(100000), ...
@(x) ~ischar(x) && isscalar(x) && (sym(x) == sym(Inf) || isPositiveInteger(sym(x))));
% Allow option 'ArrayValued' of the Core Matlab 'integral' routine, but ignore it.
p.addParameter('ArrayValued', false);
p.parse(varargin{k:end});
abstol = double(p.Results.AbsTol);
reltol = double(p.Results.RelTol);
waypoints = sym(p.Results.Waypoints);
maxfuncalls = sym(p.Results.MaxFunctionCalls);
fun = privResolveArgs(sym(fun));
fun = fun{1};
issymfun = isa(fun, 'symfun');
fargs = sym([]);
if issymfun
fargs = argnames(fun);
% Eliminate x from the arguments of the symfun
fargs = privsubsref(fargs, logical(fargs ~= x));
fun = formula(fun);
end
if isa(a, 'symfun')
fargs = union(fargs, argnames(a));
a = formula(a);
issymfun = true;
end
if isa(b, 'symfun')
fargs = union(fargs, argnames(b));
b = formula(b);
issymfun = true;
end
if isa(waypoints, 'symfun')
fargs = union(fargs, argnames(waypoints));
waypoints = formula(waypoints);
issymfun = true;
end
if isempty(a) || isempty(b)
% consistent with integral(@sin, [])
error(message('MATLAB:narginchk:notEnoughInputs'));
end
if (isinf(a) && isempty(symvar(a)) && ~isreal(a)) || ...
(isinf(b) && isempty(symvar(b)) && ~isreal(b))
error(message('symbolic:sym:vpaintegral:ComplexContourMustBeFinite'));
end
newDigits = ceil(max(2,log10(1/reltol)));
oldDigits = digits(newDigits + 2);
cleanupObj = onCleanup(@() digits(oldDigits));
if (any(isnan(reshape(fun,[],1))) || isnan(a) || isnan(b))
ysym = sym(NaN);
elseif ~isfinite(fun)
ysym = sym(fun);
else
ysym = feval(symengine, 'vpaintegral', fun,...
[x.s '=' a.s '..' b.s],...
['RelativeError = ' num2str(reltol/100)],...
['AbsoluteError = ' num2str(abstol)],...
['"Waypoints" = ' waypoints.s],...
['MaxCalls = ' maxfuncalls.s]);
if ysym == evalin(symengine, 'FAIL')
error(message('symbolic:sym:vpaintegral:IncreaseMaxFunctionCalls'))
end
end
y = privResolveOutput(ysym, fun);
% Reduce the precision of the result to reltol. Thus,
% only the reliable leading digits become visible.
digits(newDigits);
y = vpa(y);
if issymfun && ~isempty(fargs)
y = symfun(y, fargs);
end
% end of function vpaintegral
% ====================================
function [f, x, a, b] = parseFunctionVariableBounds(f, x, a, b, varargin)
%parseFunctionVariableBounds Helper method for argument parsing
% It assists in parsing inputs to a function S with calling syntax
% S(f, 'Option',...)
% S(f, x, 'Option',...)
% S(f, [a b], 'Option',...)
% S(f, a, b, 'Option',...)
% S(f, x, [a b], 'Option',...)
% S(f, x, a, b, 'Option',...)
% where x defaults to symvar(f, 1).
f = sym(f);
x = sym(x);
switch nargin
case 2
[a, b] = rangeVector(x);
if ~isempty(a)
x = symv(f);
end
case 3
a = sym(a);
[tmp, b] = rangeVector(a);
if isa(tmp, 'symfun')
tmp = formula(tmp);
end
if isempty(tmp)
b = a;
a = x;
x = symv(f);
else
a = tmp;
end
case 4
a = sym(a);
b = sym(b);
end
% ====================================
function X = symv(f)
X = symvar(f, 1);
if isempty(X)
X = sym('x');
end
% ====================================
function [a, b] = rangeVector(x)
if isa(x, 'symfun')
fargs = argnames(x);
issymfun = true;
else
issymfun = false;
end
x = formula(x);
if isvector(x) && numel(x) == 2
a = privsubsref(x,1);
b = privsubsref(x,2);
else
a = [];
b = [];
end
if issymfun && ~isempty(fargs)
a = symfun(a, fargs);
b = symfun(b, fargs);
end