日志

quadgk

已有 1363 次阅读2008-5-26 12:32 |个人分类:【MATLAB Function Reference】

Numerically evaluate integral, adaptive Gauss-Kronrod quadrature

Syntax

q = quadgk(fun,a,b) [q,errbnd] = quadgk(fun,a,b,tol) [q,errbnd] = quadgk(fun,a,b,param1,val1,param2,val2,...)

Description

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. The integrand evaluated at each element of x. fun must be a function handle. See Function Handles in the MATLAB Programming documentation for more information. 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.

, in the MATLAB Mathematics documentation, explains how to provide additional parameters to the function fun, if necessary.

[q,errbnd] = quadgk(fun,a,b,tol) returns an approximate bound on the absolute error, |Q - I|, where I denotes the exact value of the integral.

[q,errbnd] = quadgk(fun,a,b,param1,val1,param2,val2,...) performs the integration with specified values of optional parameters. The available parameters are

Parameter	Description
`'AbsTol'`	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'RelTol'= 0`. For pure relative error control use `'AbsTol' = 0`. Except when using pure absolute error control, the minimum relative tolerance is `'RelTol' >= 100eps(class(Q))`. and
`'RelTol'`	Relative error tolerance. The default value of `'RelTol'` is `1.e-6` (double), `1.e-4` (single).
`'Waypoints'`	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, the waypoints must be supplied in strictly increasing or strictly decreasing order, and only the waypoints between `a` and `b` are used. 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`.
`'MaxIntervalCount'`	Maximum number of intervals allowed. The default value is `650`.	The `'MaxIntervalCount'quadgk` uses at any one time after the first iteration. A warning is issued if `quadgkerrbnd` is small enough that the desired accuracy has nearly been achieved. parameter limits the number of intervals that returns early because of this limit. Routinely increasing this value is not recommended, but it may be appropriate when

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 xquadgk 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. approaches infinity, and
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.

Examples

Integrand with a singularity at an integration end point

Write an M-file 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:

[q,errbnd] = quadgk(@(x)x.^5.*exp(-x).*sin(x),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:

Q = quadgk(@(z)1./(2*z - 1),-1-i,-1-i,'Waypoints',[1-i,1+i,-1+i])

Q =

   0.0000 + 3.1416i

Algorithm

quadgk implements adaptive quadrature based on a Gauss-Kronrod pair (15^th and 7^th order formulas).

Diagnostics

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.