Numerically evaluate integral, adaptive Gauss-Kronrod quadrature
q = quadgk(fun,a,b)
[q,errbnd] = quadgk(fun,a,b)
[q,errbnd] = quadgk(fun,a,b,param1,val1,param2,val2,...)
q = quadgk(fun,a,b) attempts to approximate the integral of a scalar-valued function fun from a to b using high-order global adaptive quadrature and default error tolerances. The function y = fun(x) should accept a vector argument x and return a vector result y, where y is the integrand evaluated at each element of x. fun must be a function handle. Limits a and b can be -Inf or Inf. If both are finite, they can be complex. If at least one is complex, the integral is approximated over a straight line path from a to b in the complex plane.
Parameterizing Functions explains how to provide additional parameters to the function fun, if necessary.
Absolute error tolerance.
The default value of 'AbsTol' is 1.e-10 (double), 1.e-5 (single).
quadgk attempts to satisfy errbnd <= max(AbsTol,RelTol*|Q|). This is absolute error control when |Q| is sufficiently small and relative error control when |Q| is larger. For pure absolute error control use 'AbsTol' > 0 and'RelTol'= 0. For pure relative error control use 'AbsTol' = 0. Except when using pure absolute error control, the minimum relative tolerance is 'RelTol' >= 100*eps(class(Q)).
Relative error tolerance.
The default value of 'RelTol' is 1.e-6 (double), 1.e-4 (single).
Vector of integration waypoints.
If fun(x) has discontinuities in the interval of integration, the locations should be supplied as a Waypoints vector. When a, b, and the waypoints are all real, only the waypoints between a and b are used, and they are used in sorted order. Note that waypoints are not intended for singularities in fun(x). Singular points should be handled by making them endpoints of separate integrations and adding the results.
If a, b, 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 a to the first waypoint, from the first waypoint to the second, and so forth, and finally from the last waypoint to b.
Maximum number of intervals allowed.
The default value is 650.
The 'MaxIntervalCount' parameter limits the number of intervals that quadgk uses at any one time after the first iteration. A warning is issued if quadgk returns early because of this limit. Routinely increasing this value is not recommended, but it may be appropriate when errbnd is small enough that the desired accuracy has nearly been achieved.
The list below contains information to help you determine which quadrature function in MATLAB to use:
The quad function may be most efficient for low accuracies with nonsmooth integrands.
The quadl function may be more efficient than quad at higher accuracies with smooth integrands.
The quadgk function may be most efficient for high accuracies and oscillatory integrands. It supports infinite intervals and can handle moderate singularities at the endpoints. It also supports contour integration along piecewise linear paths.
The quadv function vectorizes quad for an array-valued fun.
If the interval is infinite, [a,Inf), then for the integral of fun(x) to exist, fun(x) must decay as x approaches infinity, and quadgk requires it to decay rapidly. Special methods should be used for oscillatory functions on infinite intervals, but quadgk can be used if fun(x) decays fast enough.
The quadgk function will integrate functions that are singular at finite endpoints if the singularities are not too strong. For example, it will integrate functions that behave at an endpoint c like log|x-c| or |x-c|p for p >= -1/2. If the function is singular at points inside (a,b), write the integral as a sum of integrals over subintervals with the singular points as endpoints, compute them with quadgk, and add the results.
Integrand with a singularity at an integration end point
Write a function myfun that computes the integrand:
function y = myfun(x) y = exp(x).*log(x);
Then pass @myfun, a function handle to myfun, to quadgk, along with the limits of integration, 0 to 1:
q = quadgk(@myfun,0,1) q = -1.3179
Alternatively, you can pass the integrand to quadgk as an anonymous function handle F:
f = (@(x)exp(x).*log(x)); q = quadgk(f,0,1);
Oscillatory integrand on a semi-infinite interval
Integrate over a semi-infinite interval with specified tolerances, and return the approximate error bound:
f = @(x)x.^5.*exp(-x).*sin(x); [q,errbnd] = quadgk(f,0,inf,'RelTol',1e-8,'AbsTol',1e-12) q = -15.0000 errbnd = 9.4386e-009
Contour integration around a pole
Use Waypoints to integrate around a pole using a piecewise linear contour:
f = @(z)1./(2*z - 1); q = quadgk(f,-1-i,-1-i,'Waypoints',[1-i,1+i,-1+i]) q = 0.0000 + 3.1416i
quadgk may issue one of the following warnings:
'Minimum step size reached' indicates that interval subdivision has produced a subinterval whose length is on the order of roundoff error in the length of the original interval. A nonintegrable singularity is possible.
'Reached the limit on the maximum number of intervals in use' indicates that the integration was terminated before meeting the tolerance requirements and that continuing the integration would require more than MaxIntervalCount subintervals. The integral may not exist, or it may be difficult to approximate numerically. Increasing MaxIntervalCount usually does not help unless the tolerance requirements were nearly met when the integration was previously terminated.
'Infinite or Not-a-Number function value encountered' indicates a floating point overflow or division by zero during the evaluation of the integrand in the interior of the interval.
 L.F. Shampine "Vectorized Adaptive Quadrature in MATLAB," Journal of Computational and Applied Mathematics, 211, 2008, pp.131–140.