mpmath.quad()

Algorithms

    Mpmath presently implements two integration algorithms: tanh-sinh
    quadrature and Gauss-Legendre quadrature. These can be selected
    using *method='tanh-sinh'* or *method='gauss-legendre'* or by
    passing the classes *method=TanhSinh*, *method=GaussLegendre*.
    The functions :func:`~mpmath.quadts` and :func:`~mpmath.quadgl` are also available
    as shortcuts.

    Both algorithms have the property that doubling the number of
    evaluation points roughly doubles the accuracy, so both are ideal
    for high precision quadrature (hundreds or thousands of digits).

    At high precision, computing the nodes and weights for the
    integration can be expensive (more expensive than computing the
    function values). To make repeated integrations fast, nodes
    are automatically cached.

    The advantages of the tanh-sinh algorithm are that it tends to
    handle endpoint singularities well, and that the nodes are cheap
    to compute on the first run. For these reasons, it is used by
    :func:`~mpmath.quad` as the default algorithm.

    Gauss-Legendre quadrature often requires fewer function
    evaluations, and is therefore often faster for repeated use, but
    the algorithm does not handle endpoint singularities as well and
    the nodes are more expensive to compute. Gauss-Legendre quadrature
    can be a better choice if the integrand is smooth and repeated
    integrations are required (e.g. for multiple integrals).

Examples of 1D integrals

    Intervals may be infinite or half-infinite. The following two
    examples evaluate the limits of the inverse tangent function
    (`\int 1/(1+x^2) = \tan^{-1} x`), and the Gaussian integral
    `\int_{\infty}^{\infty} \exp(-x^2)\,dx = \sqrt{\pi}`::

        >>> mp.dps = 15
        >>> quad(lambda x: 2/(x**2+1), [0, inf])
        3.14159265358979
        >>> quad(lambda x: exp(-x**2), [-inf, inf])**2
        3.14159265358979

    Integrals can typically be resolved to high precision.
    The following computes 50 digits of `\pi` by integrating the
    area of the half-circle defined by `x^2 + y^2 \le 1`,
    `-1 \le x \le 1`, `y \ge 0`::

        >>> mp.dps = 50
        >>> 2*quad(lambda x: sqrt(1-x**2), [-1, 1])
        3.1415926535897932384626433832795028841971693993751

    One can just as well compute 1000 digits (output truncated)::

        >>> mp.dps = 1000
        >>> 2*quad(lambda x: sqrt(1-x**2), [-1, 1])  #doctest:+ELLIPSIS
        3.141592653589793238462643383279502884...216420199

    Complex integrals are supported. The following computes
    a residue at `z = 0` by integrating counterclockwise along the
    diamond-shaped path from `1` to `+i` to `-1` to `-i` to `1`::

        >>> mp.dps = 15
        >>> chop(quad(lambda z: 1/z, [1,j,-1,-j,1]))
        (0.0 + 6.28318530717959j)

    **Examples of 2D and 3D integrals**

    Here are several nice examples of analytically solvable
    2D integrals (taken from MathWorld [1]) that can be evaluated
    to high precision fairly rapidly by :func:`~mpmath.quad`::

        >>> mp.dps = 30
        >>> f = lambda x, y: (x-1)/((1-x*y)*log(x*y))
        >>> quad(f, [0, 1], [0, 1])
        0.577215664901532860606512090082
        >>> +euler
        0.577215664901532860606512090082

        >>> f = lambda x, y: 1/sqrt(1+x**2+y**2)
        >>> quad(f, [-1, 1], [-1, 1])
        3.17343648530607134219175646705
        >>> 4*log(2+sqrt(3))-2*pi/3
        3.17343648530607134219175646705

        >>> f = lambda x, y: 1/(1-x**2 * y**2)
        >>> quad(f, [0, 1], [0, 1])
        1.23370055013616982735431137498
        >>> pi**2 / 8
        1.23370055013616982735431137498

        >>> quad(lambda x, y: 1/(1-x*y), [0, 1], [0, 1])
        1.64493406684822643647241516665
        >>> pi**2 / 6
        1.64493406684822643647241516665

    Multiple integrals may be done over infinite ranges::

        >>> mp.dps = 15
        >>> print(quad(lambda x,y: exp(-x-y), [0, inf], [1, inf]))
        0.367879441171442
        >>> print(1/e)
        0.367879441171442

    For nonrectangular areas, one can call :func:`~mpmath.quad` recursively.
    For example, we can replicate the earlier example of calculating
    `\pi` by integrating over the unit-circle, and actually use double
    quadrature to actually measure the area circle::

        >>> f = lambda x: quad(lambda y: 1, [-sqrt(1-x**2), sqrt(1-x**2)])
        >>> quad(f, [-1, 1])
        3.14159265358979

    Here is a simple triple integral::

        >>> mp.dps = 15
        >>> f = lambda x,y,z: x*y/(1+z)
        >>> quad(f, [0,1], [0,1], [1,2], method='gauss-legendre')
        0.101366277027041
        >>> (log(3)-log(2))/4
        0.101366277027041

Singularities

    Both tanh-sinh and Gauss-Legendre quadrature are designed to
    integrate smooth (infinitely differentiable) functions. Neither
    algorithm copes well with mid-interval singularities (such as
    mid-interval discontinuities in `f(x)` or `f'(x)`).
    The best solution is to split the integral into parts::

        >>> mp.dps = 15
        >>> quad(lambda x: abs(sin(x)), [0, 2*pi])   # Bad
        3.99900894176779
        >>> quad(lambda x: abs(sin(x)), [0, pi, 2*pi])  # Good
        4.0

    The tanh-sinh rule often works well for integrands having a
    singularity at one or both endpoints::

        >>> mp.dps = 15
        >>> quad(log, [0, 1], method='tanh-sinh')  # Good
        -1.0
        >>> quad(log, [0, 1], method='gauss-legendre')  # Bad
        -0.999932197413801

    However, the result may still be inaccurate for some functions::

        >>> quad(lambda x: 1/sqrt(x), [0, 1], method='tanh-sinh')
        1.99999999946942

    This problem is not due to the quadrature rule per se, but to
    numerical amplification of errors in the nodes. The problem can be
    circumvented by temporarily increasing the precision::

        >>> mp.dps = 30
        >>> a = quad(lambda x: 1/sqrt(x), [0, 1], method='tanh-sinh')
        >>> mp.dps = 15
        >>> +a
        2.0

Highly variable functions

    For functions that are smooth (in the sense of being infinitely
    differentiable) but contain sharp mid-interval peaks or many
    "bumps", :func:`~mpmath.quad` may fail to provide full accuracy. For
    example, with default settings, :func:`~mpmath.quad` is able to integrate
    `\sin(x)` accurately over an interval of length 100 but not over
    length 1000::

        >>> quad(sin, [0, 100]); 1-cos(100)   # Good
        0.137681127712316
        0.137681127712316
        >>> quad(sin, [0, 1000]); 1-cos(1000)   # Bad
        -37.8587612408485
        0.437620923709297

    One solution is to break the integration into 10 intervals of
    length 100::

        >>> quad(sin, linspace(0, 1000, 10))   # Good
        0.437620923709297

    Another is to increase the degree of the quadrature::

        >>> quad(sin, [0, 1000], maxdegree=10)   # Also good
        0.437620923709297

    Whether splitting the interval or increasing the degree is
    more efficient differs from case to case. Another example is the
    function `1/(1+x^2)`, which has a sharp peak centered around
    `x = 0`::

        >>> f = lambda x: 1/(1+x**2)
        >>> quad(f, [-100, 100])   # Bad
        3.64804647105268
        >>> quad(f, [-100, 100], maxdegree=10)   # Good
        3.12159332021646
        >>> quad(f, [-100, 0, 100])   # Also good
        3.12159332021646

References

    1. http://mathworld.wolfram.com/DoubleIntegral.html

存在的问题

    无法处理函数表达式中存在系数的情况。
    如 f(x,a,b) : return a * x + b
    使用lambda函数的方式解决：
    def integral_f():
    a = 3
    b = 5
    f = lamda x : a * x + b
    ret = quad(f,[-100,100])