Symbolic Integrals¶

The integrals module in SymPy implements methods to calculate definite and indefinite integrals of expressions.

Principal method in this module is integrate()

Examples¶

SymPy can integrate a vast array of functions. It can integrate polynomial functions:

>>> from sympy import *
>>> init_printing(use_unicode=False, wrap_line=False)
>>> x = Symbol('x')
>>> integrate(x**2 + x + 1, x)
 3    2
x    x
-- + -- + x
3    2

Rational functions:

>>> integrate(x/(x**2+2*x+1), x)
1
log(x + 1) + -----
x + 1

Exponential-polynomial functions. These multiplicative combinations of polynomials and the functions exp, cos and sin can be integrated by hand using repeated integration by parts, which is an extremely tedious process. Happily, SymPy will deal with these integrals.

>>> integrate(x**2 * exp(x) * cos(x), x)
 2  x           2  x                         x           x
x *e *sin(x)   x *e *cos(x)      x          e *sin(x)   e *cos(x)
------------ + ------------ - x*e *sin(x) + --------- - ---------
     2              2                           2           2

even a few nonelementary integrals (in particular, some integrals involving the error function) can be evaluated:

>>> integrate(exp(-x**2)*erf(x), x)
  ____    2
\/ pi *erf (x)
--------------
      4

Integral Transforms¶

SymPy has special support for definite integrals, and integral transforms.

sympy.integrals.transforms.mellin_transform(f, x, s, **hints)[source]¶

Compute the Mellin transform \(F(s)\) of \(f(x)\),

\[F(s) = \int_0^\infty x^{s-1} f(x) \mathrm{d}x.\]

For all “sensible” functions, this converges absolutely in a strip: \(a < \operatorname{Re}(s) < b\).

The Mellin transform is related via change of variables to the Fourier transform, and also to the (bilateral) Laplace transform.

This function returns (F, (a, b), cond) where F is the Mellin transform of f, (a, b) is the fundamental strip (as above), and cond are auxiliary convergence conditions.

If the integral cannot be computed in closed form, this function returns an unevaluated MellinTransform object.

For a description of possible hints, refer to the docstring of sympy.integrals.transforms.IntegralTransform.doit(). If noconds=False, then only \(F\) will be returned (i.e. not cond, and also not the strip (a, b)).

>>> from sympy.integrals.transforms import mellin_transform
>>> from sympy import exp
>>> from sympy.abc import x, s
>>> mellin_transform(exp(-x), x, s)
(gamma(s), (0, oo), True)

sympy.integrals.transforms.inverse_mellin_transform(F, s, x, strip, **hints)[source]¶

Compute the inverse Mellin transform of \(F(s)\) over the fundamental strip given by strip=(a, b).

This can be defined as

\[f(x) = \frac{1}{2\pi i} \int_{c - i\infty}^{c + i\infty} x^{-s} F(s) \mathrm{d}s,\]

for any \(c\) in the fundamental strip. Under certain regularity conditions on \(F\) and/or \(f\), this recovers \(f\) from its Mellin transform \(F\) (and vice versa), for positive real \(x\).

One of \(a\) or \(b\) may be passed as None; a suitable \(c\) will be inferred.

If the integral cannot be computed in closed form, this function returns an unevaluated InverseMellinTransform object.

Note that this function will assume x to be positive and real, regardless of the sympy assumptions!

For a description of possible hints, refer to the docstring of sympy.integrals.transforms.IntegralTransform.doit().

>>> from sympy.integrals.transforms import inverse_mellin_transform
>>> from sympy import oo, gamma
>>> from sympy.abc import x, s
>>> inverse_mellin_transform(gamma(s), s, x, (0, oo))
exp(-x)

The fundamental strip matters:

>>> f = 1/(s**2 - 1)
>>> inverse_mellin_transform(f, s, x, (-oo, -1))
(x/2 - 1/(2*x))*Heaviside(x - 1)
>>> inverse_mellin_transform(f, s, x, (-1, 1))
-x*Heaviside(1 - x)/2 - Heaviside(x - 1)/(2*x)
>>> inverse_mellin_transform(f, s, x, (1, oo))
(-x/2 + 1/(2*x))*Heaviside(1 - x)

sympy.integrals.transforms.laplace_transform(f, t, s, **hints)[source]¶

Compute the Laplace Transform \(F(s)\) of \(f(t)\),

\[F(s) = \int_0^\infty e^{-st} f(t) \mathrm{d}t.\]

For all “sensible” functions, this converges absolutely in a half plane \(a < \operatorname{Re}(s)\).

This function returns (F, a, cond) where F is the Laplace transform of f, \(\operatorname{Re}(s) > a\) is the half-plane of convergence, and cond are auxiliary convergence conditions.

If the integral cannot be computed in closed form, this function returns an unevaluated LaplaceTransform object.

For a description of possible hints, refer to the docstring of sympy.integrals.transforms.IntegralTransform.doit(). If noconds=True, only \(F\) will be returned (i.e. not cond, and also not the plane a).

>>> from sympy.integrals import laplace_transform
>>> from sympy.abc import t, s, a
>>> laplace_transform(t**a, t, s)
(s**(-a)*gamma(a + 1)/s, 0, re(a) > -1)

sympy.integrals.transforms.inverse_laplace_transform(F, s, t, plane=None, **hints)[source]¶

Compute the inverse Laplace transform of \(F(s)\), defined as

\[f(t) = \frac{1}{2\pi i} \int_{c-i\infty}^{c+i\infty} e^{st} F(s) \mathrm{d}s,\]

for \(c\) so large that \(F(s)\) has no singularites in the half-plane \(\operatorname{Re}(s) > c-\epsilon\).

The plane can be specified by argument plane, but will be inferred if passed as None.

Under certain regularity conditions, this recovers \(f(t)\) from its Laplace Transform \(F(s)\), for non-negative \(t\), and vice versa.

If the integral cannot be computed in closed form, this function returns an unevaluated InverseLaplaceTransform object.

Note that this function will always assume \(t\) to be real, regardless of the sympy assumption on \(t\).

For a description of possible hints, refer to the docstring of sympy.integrals.transforms.IntegralTransform.doit().

>>> from sympy.integrals.transforms import inverse_laplace_transform
>>> from sympy import exp, Symbol
>>> from sympy.abc import s, t
>>> a = Symbol('a', positive=True)
>>> inverse_laplace_transform(exp(-a*s)/s, s, t)
Heaviside(-a + t)

sympy.integrals.transforms.fourier_transform(f, x, k, **hints)[source]¶

Compute the unitary, ordinary-frequency Fourier transform of \(f\), defined as

\[F(k) = \int_{-\infty}^\infty f(x) e^{-2\pi i x k} \mathrm{d} x.\]

If the transform cannot be computed in closed form, this function returns an unevaluated FourierTransform object.

For other Fourier transform conventions, see the function sympy.integrals.transforms._fourier_transform().

For a description of possible hints, refer to the docstring of sympy.integrals.transforms.IntegralTransform.doit(). Note that for this transform, by default noconds=True.

>>> from sympy import fourier_transform, exp
>>> from sympy.abc import x, k
>>> fourier_transform(exp(-x**2), x, k)
sqrt(pi)*exp(-pi**2*k**2)
>>> fourier_transform(exp(-x**2), x, k, noconds=False)
(sqrt(pi)*exp(-pi**2*k**2), True)

sympy.integrals.transforms.inverse_fourier_transform(F, k, x, **hints)[source]¶

Compute the unitary, ordinary-frequency inverse Fourier transform of \(F\), defined as

\[f(x) = \int_{-\infty}^\infty F(k) e^{2\pi i x k} \mathrm{d} k.\]

If the transform cannot be computed in closed form, this function returns an unevaluated InverseFourierTransform object.

For other Fourier transform conventions, see the function sympy.integrals.transforms._fourier_transform().

For a description of possible hints, refer to the docstring of sympy.integrals.transforms.IntegralTransform.doit(). Note that for this transform, by default noconds=True.

>>> from sympy import inverse_fourier_transform, exp, sqrt, pi
>>> from sympy.abc import x, k
>>> inverse_fourier_transform(sqrt(pi)*exp(-(pi*k)**2), k, x)
exp(-x**2)
>>> inverse_fourier_transform(sqrt(pi)*exp(-(pi*k)**2), k, x, noconds=False)
(exp(-x**2), True)

sympy.integrals.transforms.sine_transform(f, x, k, **hints)[source]¶

Compute the unitary, ordinary-frequency sine transform of \(f\), defined as

\[F(k) = \sqrt{\frac{2}{\pi}} \int_{0}^\infty f(x) \sin(2\pi x k) \mathrm{d} x.\]

If the transform cannot be computed in closed form, this function returns an unevaluated SineTransform object.

For a description of possible hints, refer to the docstring of sympy.integrals.transforms.IntegralTransform.doit(). Note that for this transform, by default noconds=True.

>>> from sympy import sine_transform, exp
>>> from sympy.abc import x, k, a
>>> sine_transform(x*exp(-a*x**2), x, k)
sqrt(2)*k*exp(-k**2/(4*a))/(4*a**(3/2))
>>> sine_transform(x**(-a), x, k)
2**(1/2 - a)*k**(a - 1)*gamma(1 - a/2)/gamma(a/2 + 1/2)

sympy.integrals.transforms.inverse_sine_transform(F, k, x, **hints)[source]¶

Compute the unitary, ordinary-frequency inverse sine transform of \(F\), defined as

\[f(x) = \sqrt{\frac{2}{\pi}} \int_{0}^\infty F(k) \sin(2\pi x k) \mathrm{d} k.\]

If the transform cannot be computed in closed form, this function returns an unevaluated InverseSineTransform object.

For a description of possible hints, refer to the docstring of sympy.integrals.transforms.IntegralTransform.doit(). Note that for this transform, by default noconds=True.

>>> from sympy import inverse_sine_transform, exp, sqrt, gamma, pi
>>> from sympy.abc import x, k, a
>>> inverse_sine_transform(2**((1-2*a)/2)*k**(a - 1)*
...     gamma(-a/2 + 1)/gamma((a+1)/2), k, x)
x**(-a)
>>> inverse_sine_transform(sqrt(2)*k*exp(-k**2/(4*a))/(4*sqrt(a)**3), k, x)
x*exp(-a*x**2)

sympy.integrals.transforms.cosine_transform(f, x, k, **hints)[source]¶

Compute the unitary, ordinary-frequency cosine transform of \(f\), defined as

\[F(k) = \sqrt{\frac{2}{\pi}} \int_{0}^\infty f(x) \cos(2\pi x k) \mathrm{d} x.\]

If the transform cannot be computed in closed form, this function returns an unevaluated CosineTransform object.

For a description of possible hints, refer to the docstring of sympy.integrals.transforms.IntegralTransform.doit(). Note that for this transform, by default noconds=True.

>>> from sympy import cosine_transform, exp, sqrt, cos
>>> from sympy.abc import x, k, a
>>> cosine_transform(exp(-a*x), x, k)
sqrt(2)*a/(sqrt(pi)*(a**2 + k**2))
>>> cosine_transform(exp(-a*sqrt(x))*cos(a*sqrt(x)), x, k)
a*exp(-a**2/(2*k))/(2*k**(3/2))

sympy.integrals.transforms.inverse_cosine_transform(F, k, x, **hints)[source]¶

Compute the unitary, ordinary-frequency inverse cosine transform of \(F\), defined as

\[f(x) = \sqrt{\frac{2}{\pi}} \int_{0}^\infty F(k) \cos(2\pi x k) \mathrm{d} k.\]

If the transform cannot be computed in closed form, this function returns an unevaluated InverseCosineTransform object.

For a description of possible hints, refer to the docstring of sympy.integrals.transforms.IntegralTransform.doit(). Note that for this transform, by default noconds=True.

>>> from sympy import inverse_cosine_transform, exp, sqrt, pi
>>> from sympy.abc import x, k, a
>>> inverse_cosine_transform(sqrt(2)*a/(sqrt(pi)*(a**2 + k**2)), k, x)
exp(-a*x)
>>> inverse_cosine_transform(1/sqrt(k), k, x)
1/sqrt(x)

sympy.integrals.transforms.hankel_transform(f, r, k, nu, **hints)[source]¶

Compute the Hankel transform of \(f\), defined as

\[F_\nu(k) = \int_{0}^\infty f(r) J_\nu(k r) r \mathrm{d} r.\]

If the transform cannot be computed in closed form, this function returns an unevaluated HankelTransform object.

For a description of possible hints, refer to the docstring of sympy.integrals.transforms.IntegralTransform.doit(). Note that for this transform, by default noconds=True.

>>> from sympy import hankel_transform, inverse_hankel_transform
>>> from sympy import gamma, exp, sinh, cosh
>>> from sympy.abc import r, k, m, nu, a

>>> ht = hankel_transform(1/r**m, r, k, nu)
>>> ht
2*2**(-m)*k**(m - 2)*gamma(-m/2 + nu/2 + 1)/gamma(m/2 + nu/2)

>>> inverse_hankel_transform(ht, k, r, nu)
r**(-m)

>>> ht = hankel_transform(exp(-a*r), r, k, 0)
>>> ht
a/(k**3*(a**2/k**2 + 1)**(3/2))

>>> inverse_hankel_transform(ht, k, r, 0)
exp(-a*r)

sympy.integrals.transforms.inverse_hankel_transform(F, k, r, nu, **hints)[source]¶

Compute the inverse Hankel transform of \(F\) defined as

\[f(r) = \int_{0}^\infty F_\nu(k) J_\nu(k r) k \mathrm{d} k.\]

If the transform cannot be computed in closed form, this function returns an unevaluated InverseHankelTransform object.

For a description of possible hints, refer to the docstring of sympy.integrals.transforms.IntegralTransform.doit(). Note that for this transform, by default noconds=True.

>>> from sympy import hankel_transform, inverse_hankel_transform, gamma
>>> from sympy import gamma, exp, sinh, cosh
>>> from sympy.abc import r, k, m, nu, a

>>> ht = hankel_transform(1/r**m, r, k, nu)
>>> ht
2*2**(-m)*k**(m - 2)*gamma(-m/2 + nu/2 + 1)/gamma(m/2 + nu/2)

>>> inverse_hankel_transform(ht, k, r, nu)
r**(-m)

>>> ht = hankel_transform(exp(-a*r), r, k, 0)
>>> ht
a/(k**3*(a**2/k**2 + 1)**(3/2))

>>> inverse_hankel_transform(ht, k, r, 0)
exp(-a*r)

Internals¶

There is a general method for calculating antiderivatives of elementary functions, called the Risch algorithm. The Risch algorithm is a decision procedure that can determine whether an elementary solution exists, and in that case calculate it. It can be extended to handle many nonelementary functions in addition to the elementary ones.

SymPy currently uses a simplified version of the Risch algorithm, called the Risch-Norman algorithm. This algorithm is much faster, but may fail to find an antiderivative, although it is still very powerful. SymPy also uses pattern matching and heuristics to speed up evaluation of some types of integrals, e.g. polynomials.

For non-elementary definite integrals, SymPy uses so-called Meijer G-functions. Details are described here:

API reference¶

sympy.integrals.integrate(f, var, ...)[source]¶

Compute definite or indefinite integral of one or more variables using Risch-Norman algorithm and table lookup. This procedure is able to handle elementary algebraic and transcendental functions and also a huge class of special functions, including Airy, Bessel, Whittaker and Lambert.

var can be:

a symbol – indefinite integration
a tuple (symbol, a) – indefinite integration with result
given with \(a\) replacing \(symbol\)
a tuple (symbol, a, b) – definite integration

Several variables can be specified, in which case the result is multiple integration. (If var is omitted and the integrand is univariate, the indefinite integral in that variable will be performed.)

Indefinite integrals are returned without terms that are independent of the integration variables. (see examples)

Definite improper integrals often entail delicate convergence conditions. Pass conds=’piecewise’, ‘separate’ or ‘none’ to have these returned, respectively, as a Piecewise function, as a separate result (i.e. result will be a tuple), or not at all (default is ‘piecewise’).

Strategy

SymPy uses various approaches to definite integration. One method is to find an antiderivative for the integrand, and then use the fundamental theorem of calculus. Various functions are implemented to integrate polynomial, rational and trigonometric functions, and integrands containing DiracDelta terms.

SymPy also implements the part of the Risch algorithm, which is a decision procedure for integrating elementary functions, i.e., the algorithm can either find an elementary antiderivative, or prove that one does not exist. There is also a (very successful, albeit somewhat slow) general implementation of the heuristic Risch algorithm. This algorithm will eventually be phased out as more of the full Risch algorithm is implemented. See the docstring of Integral._eval_integral() for more details on computing the antiderivative using algebraic methods.

The option risch=True can be used to use only the (full) Risch algorithm. This is useful if you want to know if an elementary function has an elementary antiderivative. If the indefinite Integral returned by this function is an instance of NonElementaryIntegral, that means that the Risch algorithm has proven that integral to be non-elementary. Note that by default, additional methods (such as the Meijer G method outlined below) are tried on these integrals, as they may be expressible in terms of special functions, so if you only care about elementary answers, use risch=True. Also note that an unevaluated Integral returned by this function is not necessarily a NonElementaryIntegral, even with risch=True, as it may just be an indication that the particular part of the Risch algorithm needed to integrate that function is not yet implemented.

Another family of strategies comes from re-writing the integrand in terms of so-called Meijer G-functions. Indefinite integrals of a single G-function can always be computed, and the definite integral of a product of two G-functions can be computed from zero to infinity. Various strategies are implemented to rewrite integrands as G-functions, and use this information to compute integrals (see the meijerint module).

The option manual=True can be used to use only an algorithm that tries to mimic integration by hand. This algorithm does not handle as many integrands as the other algorithms implemented but may return results in a more familiar form. The manualintegrate module has functions that return the steps used (see the module docstring for more information).

In general, the algebraic methods work best for computing antiderivatives of (possibly complicated) combinations of elementary functions. The G-function methods work best for computing definite integrals from zero to infinity of moderately complicated combinations of special functions, or indefinite integrals of very simple combinations of special functions.

The strategy employed by the integration code is as follows:

If computing a definite integral, and both limits are real, and at least one limit is +- oo, try the G-function method of definite integration first.
Try to find an antiderivative, using all available methods, ordered by performance (that is try fastest method first, slowest last; in particular polynomial integration is tried first, Meijer G-functions second to last, and heuristic Risch last).
If still not successful, try G-functions irrespective of the limits.

The option meijerg=True, False, None can be used to, respectively: always use G-function methods and no others, never use G-function methods, or use all available methods (in order as described above). It defaults to None.

Examples

>>> from sympy import integrate, log, exp, oo
>>> from sympy.abc import a, x, y

>>> integrate(x*y, x)
x**2*y/2

>>> integrate(log(x), x)
x*log(x) - x

>>> integrate(log(x), (x, 1, a))
a*log(a) - a + 1

>>> integrate(x)
x**2/2

Terms that are independent of x are dropped by indefinite integration:

>>> from sympy import sqrt
>>> integrate(sqrt(1 + x), (x, 0, x))
2*(x + 1)**(3/2)/3 - 2/3
>>> integrate(sqrt(1 + x), x)
2*(x + 1)**(3/2)/3

>>> integrate(x*y)
Traceback (most recent call last):
...
ValueError: specify integration variables to integrate x*y

Note that integrate(x) syntax is meant only for convenience in interactive sessions and should be avoided in library code.

>>> integrate(x**a*exp(-x), (x, 0, oo)) # same as conds='piecewise'
Piecewise((gamma(a + 1), re(a) > -1),
    (Integral(x**a*exp(-x), (x, 0, oo)), True))

>>> integrate(x**a*exp(-x), (x, 0, oo), conds='none')
gamma(a + 1)

>>> integrate(x**a*exp(-x), (x, 0, oo), conds='separate')
(gamma(a + 1), -re(a) < 1)

See also

Integral, Integral.doit

sympy.integrals.line_integrate(field, Curve, variables)[source]¶

Compute the line integral.

Examples

>>> from sympy import Curve, line_integrate, E, ln
>>> from sympy.abc import x, y, t
>>> C = Curve([E**t + 1, E**t - 1], (t, 0, ln(2)))
>>> line_integrate(x + y, C, [x, y])
3*sqrt(2)

See also

integrate, Integral

The class \(Integral\) represents an unevaluated integral and has some methods that help in the integration of an expression.

class sympy.integrals.Integral[source]¶

Represents unevaluated integral.

is_commutative¶: Returns whether all the free symbols in the integral are commutative.

as_sum(n=None, method='midpoint', evaluate=True)[source]¶

Approximates a definite integral by a sum.

Arguments

n: The number of subintervals to use, optional.
method: One of: ‘left’, ‘right’, ‘midpoint’, ‘trapezoid’.
evaluate: If False, returns an unevaluated Sum expression. The default is True, evaluate the sum.

These methods of approximate integration are described in [1].

[1] https://en.wikipedia.org/wiki/Riemann_sum#Methods

Examples

>>> from sympy import sin, sqrt
>>> from sympy.abc import x, n
>>> from sympy.integrals import Integral
>>> e = Integral(sin(x), (x, 3, 7))
>>> e
Integral(sin(x), (x, 3, 7))

For demonstration purposes, this interval will only be split into 2 regions, bounded by [3, 5] and [5, 7].

The left-hand rule uses function evaluations at the left of each interval:

>>> e.as_sum(2, 'left')
2*sin(5) + 2*sin(3)

The midpoint rule uses evaluations at the center of each interval:

>>> e.as_sum(2, 'midpoint')
2*sin(4) + 2*sin(6)

The right-hand rule uses function evaluations at the right of each interval:

>>> e.as_sum(2, 'right')
2*sin(5) + 2*sin(7)

The trapezoid rule uses function evaluations on both sides of the intervals. This is equivalent to taking the average of the left and right hand rule results:

>>> e.as_sum(2, 'trapezoid')
2*sin(5) + sin(3) + sin(7)
>>> (e.as_sum(2, 'left') + e.as_sum(2, 'right'))/2 == _
True

Here, the discontinuity at x = 0 can be avoided by using the midpoint or right-hand method:

>>> e = Integral(1/sqrt(x), (x, 0, 1))
>>> e.as_sum(5).n(4)
1.730
>>> e.as_sum(10).n(4)
1.809
>>> e.doit().n(4)  # the actual value is 2
2.000

The left- or trapezoid method will encounter the discontinuity and return infinity:

>>> e.as_sum(5, 'left')
zoo

The number of intervals can be symbolic. If omitted, a dummy symbol will be used for it. >>> e = Integral(x**2, (x, 0, 2)) >>> e.as_sum(n, ‘right’).expand() 8/3 + 4/n + 4/(3*n**2)

This shows that the midpoint rule is more accurate, as its error term decays as the square of n: >>> e.as_sum(method=’midpoint’).expand() 8/3 - 2/(3*_n**2)

A symbolic sum is returned with evaluate=False: >>> e.as_sum(n, ‘midpoint’, evaluate=False) 2*Sum((2*_k/n - 1/n)**2, (_k, 1, n))/n

See also

Integral.doit: Perform the integration using any hints

doit(**hints)[source]¶

Perform the integration using any hints given.

Examples

>>> from sympy import Integral
>>> from sympy.abc import x, i
>>> Integral(x**i, (i, 1, 3)).doit()
Piecewise((x**3/log(x) - x/log(x),
    (x > 1) | ((x >= 0) & (x < 1))), (2, True))

See also

sympy.integrals.trigonometry.trigintegrate, sympy.integrals.risch.heurisch, sympy.integrals.rationaltools.ratint

as_sum: Approximate the integral using a sum

free_symbols¶

This method returns the symbols that will exist when the integral is evaluated. This is useful if one is trying to determine whether an integral depends on a certain symbol or not.

Examples

>>> from sympy import Integral
>>> from sympy.abc import x, y
>>> Integral(x, (x, y, 1)).free_symbols
{y}

See also

function, limits, variables

principal_value(**kwargs)[source]¶

Compute the Cauchy Principal Value of the definite integral of a real function in the given interval on the real axis. In mathematics, the Cauchy principal value, is a method for assigning values to certain improper integrals which would otherwise be undefined.

Examples

>>> from sympy import Dummy, symbols, integrate, limit, oo
>>> from sympy.integrals.integrals import Integral
>>> from sympy.calculus.singularities import singularities
>>> x = symbols('x')
>>> Integral(x+1, (x, -oo, oo)).principal_value()
oo
>>> f = 1 / (x**3)
>>> Integral(f, (x, -oo, oo)).principal_value()
0
>>> Integral(f, (x, -10, 10)).principal_value()
0
>>> Integral(f, (x, -10, oo)).principal_value() + Integral(f, (x, -oo, 10)).principal_value()
0

References

R435: https://en.wikipedia.org/wiki/Cauchy_principal_value
R436: http://mathworld.wolfram.com/CauchyPrincipalValue.html

transform(x, u)[source]¶

Performs a change of variables from \(x\) to \(u\) using the relationship given by \(x\) and \(u\) which will define the transformations \(f\) and \(F\) (which are inverses of each other) as follows:

If \(x\) is a Symbol (which is a variable of integration) then \(u\) will be interpreted as some function, f(u), with inverse F(u). This, in effect, just makes the substitution of x with f(x).
If \(u\) is a Symbol then \(x\) will be interpreted as some function, F(x), with inverse f(u). This is commonly referred to as u-substitution.

Once f and F have been identified, the transformation is made as follows:

\[\int_a^b x \mathrm{d}x \rightarrow \int_{F(a)}^{F(b)} f(x) \frac{\mathrm{d}}{\mathrm{d}x}\]

where \(F(x)\) is the inverse of \(f(x)\) and the limits and integrand have been corrected so as to retain the same value after integration.

Notes

The mappings, F(x) or f(u), must lead to a unique integral. Linear or rational linear expression, \(2*x\), \(1/x\) and \(sqrt(x)\), will always work; quadratic expressions like \(x**2 - 1\) are acceptable as long as the resulting integrand does not depend on the sign of the solutions (see examples).

The integral will be returned unchanged if \(x\) is not a variable of integration.

\(x\) must be (or contain) only one of of the integration variables. If \(u\) has more than one free symbol then it should be sent as a tuple (\(u\), \(uvar\)) where \(uvar\) identifies which variable is replacing the integration variable. XXX can it contain another integration variable?

Examples

>>> from sympy.abc import a, b, c, d, x, u, y
>>> from sympy import Integral, S, cos, sqrt

>>> i = Integral(x*cos(x**2 - 1), (x, 0, 1))

transform can change the variable of integration

>>> i.transform(x, u)
Integral(u*cos(u**2 - 1), (u, 0, 1))

transform can perform u-substitution as long as a unique integrand is obtained:

>>> i.transform(x**2 - 1, u)
Integral(cos(u)/2, (u, -1, 0))

This attempt fails because x = +/-sqrt(u + 1) and the sign does not cancel out of the integrand:

>>> Integral(cos(x**2 - 1), (x, 0, 1)).transform(x**2 - 1, u)
Traceback (most recent call last):
...
ValueError:
The mapping between F(x) and f(u) did not give a unique integrand.

transform can do a substitution. Here, the previous result is transformed back into the original expression using “u-substitution”:

>>> ui = _
>>> _.transform(sqrt(u + 1), x) == i
True

We can accomplish the same with a regular substitution:

>>> ui.transform(u, x**2 - 1) == i
True

If the \(x\) does not contain a symbol of integration then the integral will be returned unchanged. Integral \(i\) does not have an integration variable \(a\) so no change is made:

>>> i.transform(a, x) == i
True

When \(u\) has more than one free symbol the symbol that is replacing \(x\) must be identified by passing \(u\) as a tuple:

>>> Integral(x, (x, 0, 1)).transform(x, (u + a, u))
Integral(a + u, (u, -a, 1 - a))
>>> Integral(x, (x, 0, 1)).transform(x, (u + a, a))
Integral(a + u, (a, -u, 1 - u))

See also

variables: Lists the integration variables
as_dummy: Replace integration variables with dummy ones

TODO and Bugs¶

There are still lots of functions that SymPy does not know how to integrate. For bugs related to this module, see https://github.com/sympy/sympy/issues?q=label%3AIntegration

Numeric Integrals¶

SymPy has functions to calculate points and weights for Gaussian quadrature of any order and any precision:

sympy.integrals.quadrature.gauss_legendre(n, n_digits)[source]¶

Computes the Gauss-Legendre quadrature [R437] points and weights.

The Gauss-Legendre quadrature approximates the integral:

\[\int_{-1}^1 f(x)\,dx \approx \sum_{i=1}^n w_i f(x_i)\]

The nodes \(x_i\) of an order \(n\) quadrature rule are the roots of \(P_n\) and the weights \(w_i\) are given by:

\[w_i = \frac{2}{\left(1-x_i^2\right) \left(P'_n(x_i)\right)^2}\]

Parameters

n : the order of quadrature

n_digits : number of significant digits of the points and weights to return

Returns

(x, w) : the x and w are lists of points and weights as Floats.

Examples

>>> from sympy.integrals.quadrature import gauss_legendre
>>> x, w = gauss_legendre(3, 5)
>>> x
[-0.7746, 0, 0.7746]
>>> w
[0.55556, 0.88889, 0.55556]
>>> x, w = gauss_legendre(4, 5)
>>> x
[-0.86114, -0.33998, 0.33998, 0.86114]
>>> w
[0.34785, 0.65215, 0.65215, 0.34785]

References

R437(1,2): https://en.wikipedia.org/wiki/Gaussian_quadrature
R438: http://people.sc.fsu.edu/~jburkardt/cpp_src/legendre_rule/legendre_rule.html

sympy.integrals.quadrature.gauss_laguerre(n, n_digits)[source]¶

Computes the Gauss-Laguerre quadrature [R439] points and weights.

The Gauss-Laguerre quadrature approximates the integral:

\[\int_0^{\infty} e^{-x} f(x)\,dx \approx \sum_{i=1}^n w_i f(x_i)\]

The nodes \(x_i\) of an order \(n\) quadrature rule are the roots of \(L_n\) and the weights \(w_i\) are given by:

\[w_i = \frac{x_i}{(n+1)^2 \left(L_{n+1}(x_i)\right)^2}\]

Parameters

n : the order of quadrature

n_digits : number of significant digits of the points and weights to return

Returns

(x, w) : the x and w are lists of points and weights as Floats.

Examples

>>> from sympy.integrals.quadrature import gauss_laguerre
>>> x, w = gauss_laguerre(3, 5)
>>> x
[0.41577, 2.2943, 6.2899]
>>> w
[0.71109, 0.27852, 0.010389]
>>> x, w = gauss_laguerre(6, 5)
>>> x
[0.22285, 1.1889, 2.9927, 5.7751, 9.8375, 15.983]
>>> w
[0.45896, 0.417, 0.11337, 0.010399, 0.00026102, 8.9855e-7]

References

R439(1,2): https://en.wikipedia.org/wiki/Gauss%E2%80%93Laguerre_quadrature
R440: http://people.sc.fsu.edu/~jburkardt/cpp_src/laguerre_rule/laguerre_rule.html

sympy.integrals.quadrature.gauss_hermite(n, n_digits)[source]¶

Computes the Gauss-Hermite quadrature [R441] points and weights.

The Gauss-Hermite quadrature approximates the integral:

\[\int_{-\infty}^{\infty} e^{-x^2} f(x)\,dx \approx \sum_{i=1}^n w_i f(x_i)\]

The nodes \(x_i\) of an order \(n\) quadrature rule are the roots of \(H_n\) and the weights \(w_i\) are given by:

\[w_i = \frac{2^{n-1} n! \sqrt{\pi}}{n^2 \left(H_{n-1}(x_i)\right)^2}\]

Parameters

n : the order of quadrature

n_digits : number of significant digits of the points and weights to return

Returns

(x, w) : the x and w are lists of points and weights as Floats.

Examples

>>> from sympy.integrals.quadrature import gauss_hermite
>>> x, w = gauss_hermite(3, 5)
>>> x
[-1.2247, 0, 1.2247]
>>> w
[0.29541, 1.1816, 0.29541]

>>> x, w = gauss_hermite(6, 5)
>>> x
[-2.3506, -1.3358, -0.43608, 0.43608, 1.3358, 2.3506]
>>> w
[0.00453, 0.15707, 0.72463, 0.72463, 0.15707, 0.00453]

References

R441(1,2): https://en.wikipedia.org/wiki/Gauss-Hermite_Quadrature
R442: http://people.sc.fsu.edu/~jburkardt/cpp_src/hermite_rule/hermite_rule.html
R443: http://people.sc.fsu.edu/~jburkardt/cpp_src/gen_hermite_rule/gen_hermite_rule.html

sympy.integrals.quadrature.gauss_gen_laguerre(n, alpha, n_digits)[source]¶

Computes the generalized Gauss-Laguerre quadrature [R444] points and weights.

The generalized Gauss-Laguerre quadrature approximates the integral:

\[\int_{0}^\infty x^{\alpha} e^{-x} f(x)\,dx \approx \sum_{i=1}^n w_i f(x_i)\]

The nodes \(x_i\) of an order \(n\) quadrature rule are the roots of \(L^{\alpha}_n\) and the weights \(w_i\) are given by:

\[w_i = \frac{\Gamma(\alpha+n)} {n \Gamma(n) L^{\alpha}_{n-1}(x_i) L^{\alpha+1}_{n-1}(x_i)}\]

Parameters

n : the order of quadrature

alpha : the exponent of the singularity, \(\alpha > -1\)

n_digits : number of significant digits of the points and weights to return

Returns

(x, w) : the x and w are lists of points and weights as Floats.

Examples

>>> from sympy import S
>>> from sympy.integrals.quadrature import gauss_gen_laguerre
>>> x, w = gauss_gen_laguerre(3, -S.Half, 5)
>>> x
[0.19016, 1.7845, 5.5253]
>>> w
[1.4493, 0.31413, 0.00906]

>>> x, w = gauss_gen_laguerre(4, 3*S.Half, 5)
>>> x
[0.97851, 2.9904, 6.3193, 11.712]
>>> w
[0.53087, 0.67721, 0.11895, 0.0023152]

References

R444(1,2): https://en.wikipedia.org/wiki/Gauss%E2%80%93Laguerre_quadrature
R445: http://people.sc.fsu.edu/~jburkardt/cpp_src/gen_laguerre_rule/gen_laguerre_rule.html

sympy.integrals.quadrature.gauss_chebyshev_t(n, n_digits)[source]¶

Computes the Gauss-Chebyshev quadrature [R446] points and weights of the first kind.

The Gauss-Chebyshev quadrature of the first kind approximates the integral:

\[\int_{-1}^{1} \frac{1}{\sqrt{1-x^2}} f(x)\,dx \approx \sum_{i=1}^n w_i f(x_i)\]

The nodes \(x_i\) of an order \(n\) quadrature rule are the roots of \(T_n\) and the weights \(w_i\) are given by:

\[w_i = \frac{\pi}{n}\]

Parameters

n : the order of quadrature

n_digits : number of significant digits of the points and weights to return

Returns

(x, w) : the x and w are lists of points and weights as Floats.

Examples

>>> from sympy import S
>>> from sympy.integrals.quadrature import gauss_chebyshev_t
>>> x, w = gauss_chebyshev_t(3, 5)
>>> x
[0.86602, 0, -0.86602]
>>> w
[1.0472, 1.0472, 1.0472]

>>> x, w = gauss_chebyshev_t(6, 5)
>>> x
[0.96593, 0.70711, 0.25882, -0.25882, -0.70711, -0.96593]
>>> w
[0.5236, 0.5236, 0.5236, 0.5236, 0.5236, 0.5236]

References

R446(1,2): https://en.wikipedia.org/wiki/Chebyshev%E2%80%93Gauss_quadrature
R447: http://people.sc.fsu.edu/~jburkardt/cpp_src/chebyshev1_rule/chebyshev1_rule.html

sympy.integrals.quadrature.gauss_chebyshev_u(n, n_digits)[source]¶

Computes the Gauss-Chebyshev quadrature [R448] points and weights of the second kind.

The Gauss-Chebyshev quadrature of the second kind approximates the integral:

\[\int_{-1}^{1} \sqrt{1-x^2} f(x)\,dx \approx \sum_{i=1}^n w_i f(x_i)\]

The nodes \(x_i\) of an order \(n\) quadrature rule are the roots of \(U_n\) and the weights \(w_i\) are given by:

\[w_i = \frac{\pi}{n+1} \sin^2 \left(\frac{i}{n+1}\pi\right)\]

Parameters

n : the order of quadrature

n_digits : number of significant digits of the points and weights to return

Returns

(x, w) : the x and w are lists of points and weights as Floats.

Examples

>>> from sympy import S
>>> from sympy.integrals.quadrature import gauss_chebyshev_u
>>> x, w = gauss_chebyshev_u(3, 5)
>>> x
[0.70711, 0, -0.70711]
>>> w
[0.3927, 0.7854, 0.3927]

>>> x, w = gauss_chebyshev_u(6, 5)
>>> x
[0.90097, 0.62349, 0.22252, -0.22252, -0.62349, -0.90097]
>>> w
[0.084489, 0.27433, 0.42658, 0.42658, 0.27433, 0.084489]

References

R448(1,2): https://en.wikipedia.org/wiki/Chebyshev%E2%80%93Gauss_quadrature
R449: http://people.sc.fsu.edu/~jburkardt/cpp_src/chebyshev2_rule/chebyshev2_rule.html

sympy.integrals.quadrature.gauss_jacobi(n, alpha, beta, n_digits)[source]¶

Computes the Gauss-Jacobi quadrature [R450] points and weights.

The Gauss-Jacobi quadrature of the first kind approximates the integral:

\[\int_{-1}^1 (1-x)^\alpha (1+x)^\beta f(x)\,dx \approx \sum_{i=1}^n w_i f(x_i)\]

The nodes \(x_i\) of an order \(n\) quadrature rule are the roots of \(P^{(\alpha,\beta)}_n\) and the weights \(w_i\) are given by:

\[w_i = -\frac{2n+\alpha+\beta+2}{n+\alpha+\beta+1} \frac{\Gamma(n+\alpha+1)\Gamma(n+\beta+1)} {\Gamma(n+\alpha+\beta+1)(n+1)!} \frac{2^{\alpha+\beta}}{P'_n(x_i) P^{(\alpha,\beta)}_{n+1}(x_i)}\]

Parameters

n : the order of quadrature

alpha : the first parameter of the Jacobi Polynomial, \(\alpha > -1\)

beta : the second parameter of the Jacobi Polynomial, \(\beta > -1\)

n_digits : number of significant digits of the points and weights to return

Returns

(x, w) : the x and w are lists of points and weights as Floats.

Examples

>>> from sympy import S
>>> from sympy.integrals.quadrature import gauss_jacobi
>>> x, w = gauss_jacobi(3, S.Half, -S.Half, 5)
>>> x
[-0.90097, -0.22252, 0.62349]
>>> w
[1.7063, 1.0973, 0.33795]

>>> x, w = gauss_jacobi(6, 1, 1, 5)
>>> x
[-0.87174, -0.5917, -0.2093, 0.2093, 0.5917, 0.87174]
>>> w
[0.050584, 0.22169, 0.39439, 0.39439, 0.22169, 0.050584]

References

R450(1,2): https://en.wikipedia.org/wiki/Gauss%E2%80%93Jacobi_quadrature
R451: http://people.sc.fsu.edu/~jburkardt/cpp_src/jacobi_rule/jacobi_rule.html
R452: http://people.sc.fsu.edu/~jburkardt/cpp_src/gegenbauer_rule/gegenbauer_rule.html

sympy.integrals.quadrature.gauss_lobatto(n, n_digits)[source]¶

Computes the Gauss-Lobatto quadrature [R453] points and weights.

The Gauss-Lobatto quadrature approximates the integral:

\[\int_{-1}^1 f(x)\,dx \approx \sum_{i=1}^n w_i f(x_i)\]

The nodes \(x_i\) of an order \(n\) quadrature rule are the roots of \(P'_(n-1)\) and the weights \(w_i\) are given by:

\[\begin{split}&w_i = \frac{2}{n(n-1) \left[P_{n-1}(x_i)\right]^2},\quad x\neq\pm 1\\ &w_i = \frac{2}{n(n-1)},\quad x=\pm 1\end{split}\]

Parameters

n : the order of quadrature

n_digits : number of significant digits of the points and weights to return

Returns

(x, w) : the x and w are lists of points and weights as Floats.

Examples

>>> from sympy.integrals.quadrature import gauss_lobatto
>>> x, w = gauss_lobatto(3, 5)
>>> x
[-1, 0, 1]
>>> w
[0.33333, 1.3333, 0.33333]
>>> x, w = gauss_lobatto(4, 5)
>>> x
[-1, -0.44721, 0.44721, 1]
>>> w
[0.16667, 0.83333, 0.83333, 0.16667]

References

R453(1,2): https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Lobatto_rules
R454: http://people.math.sfu.ca/~cbm/aands/page_888.htm

Integration over Polytopes¶

The intpoly module in SymPy implements methods to calculate the integral of a polynomial over 2/3-Polytopes. Uses evaluation techniques as described in Chin et al. (2015) [1].

The input for 2-Polytope or Polygon uses the already existing Polygon data structure in SymPy. See sympy.geometry.polygon for how to create a polygon.

For the 3-Polytope or Polyhedron, the most economical representation is to specify a list of vertices and then to provide each constituting face(Polygon) as a list of vertex indices.

For example, consider the unit cube. Here is how it would be represented.

unit_cube = [[(0, 0, 0), (0, 0, 1), (0, 1, 0), (0, 1, 1), (1, 0, 0),(1, 0, 1), (1, 1, 0), (1, 1, 1)],: [3, 7, 6, 2], [1, 5, 7, 3], [5, 4, 6, 7], [0, 4, 5, 1], [2, 0, 1, 3], [2, 6, 4, 0]]

Here, the first sublist is the list of vertices. The other smaller lists such as [3, 7, 6, 2] represent a 2D face of the polyhedra with vertices having index 3, 7, 6 and 2 in the first sublist(in that order).

Principal method in this module is polytope_integrate()

References¶

[1] : Chin, Eric B., Jean B. Lasserre, and N. Sukumar. “Numerical integration of homogeneous functions on convex and nonconvex polygons and polyhedra.” Computational Mechanics 56.6 (2015): 967-981

PDF link : http://dilbert.engr.ucdavis.edu/~suku/quadrature/cls-integration.pdf

Examples¶

For 2D Polygons¶

Single Polynomial:

>>> from sympy.integrals.intpoly import *
>>> init_printing(use_unicode=False, wrap_line=False)
>>> polytope_integrate(Polygon((0, 0), (0, 1), (1, 0)), x)
1/6
>>> polytope_integrate(Polygon((0, 0), (0, 1), (1, 0)), x + x*y + y**2)
7/24

List of specified polynomials:

>>> polytope_integrate(Polygon((0, 0), (0, 1), (1, 0)), [3, x*y + y**2, x**4], max_degree=4)
          4               2
{3: 3/2, x : 1/30, x*y + y : 1/8}
>>> polytope_integrate(Polygon((0, 0), (0, 1), (1, 0)), [1.125, x, x**2, 6.89*x**3, x*y + y**2, x**4], max_degree=4)
                       2              3  689    4               2
{1.125: 9/16, x: 1/6, x : 1/12, 6.89*x : ----, x : 1/30, x*y + y : 1/8}
                                         2000

Computing all monomials up to a maximum degree:

>>> polytope_integrate(Polygon((0, 0), (0, 1), (1, 0)),max_degree=3)
                        2         3                 2         3                      2         2
{0: 0, 1: 1/2, x: 1/6, x : 1/12, x : 1/20, y: 1/6, y : 1/12, y : 1/20, x*y: 1/24, x*y : 1/60, x *y: 1/60}

For 3-Polytopes/Polyhedra¶

Single Polynomial:

>>> from sympy.integrals.intpoly import *
>>> cube = [[(0, 0, 0), (0, 0, 5), (0, 5, 0), (0, 5, 5), (5, 0, 0), (5, 0, 5), (5, 5, 0), (5, 5, 5)], [2, 6, 7, 3], [3, 7, 5, 1], [7, 6, 4, 5], [1, 5, 4, 0], [3, 1, 0, 2], [0, 4, 6, 2]]
>>> polytope_integrate(cube, x**2 + y**2 + z**2 + x*y + y*z + x*z)
-21875/4
>>> octahedron = [[(S(-1) / sqrt(2), 0, 0), (0, S(1) / sqrt(2), 0), (0, 0, S(-1) / sqrt(2)), (0, 0, S(1) / sqrt(2)), (0, S(-1) / sqrt(2), 0), (S(1) / sqrt(2), 0, 0)], [3, 4, 5], [3, 5, 1], [3, 1, 0], [3, 0, 4], [4, 0, 2], [4, 2, 5], [2, 0, 1], [5, 2, 1]]
>>> polytope_integrate(octahedron, x**2 + y**2 + z**2 + x*y + y*z + x*z)
  ___
\/ 2
-----
  20

List of specified polynomials:

>>> polytope_integrate(Polygon((0, 0), (0, 1), (1, 0)), [3, x*y + y**2, x**4], max_degree=4)
          4               2
{3: 3/2, x : 1/30, x*y + y : 1/8}
>>> polytope_integrate(Polygon((0, 0), (0, 1), (1, 0)), [1.125, x, x**2, 6.89*x**3, x*y + y**2, x**4], max_degree=4)
                       2              3  689    4               2
{1.125: 9/16, x: 1/6, x : 1/12, 6.89*x : ----, x : 1/30, x*y + y : 1/8}
                                         2000

Computing all monomials up to a maximum degree:

>>> polytope_integrate(Polygon((0, 0), (0, 1), (1, 0)),max_degree=3)
                        2         3                 2         3                      2         2
{0: 0, 1: 1/2, x: 1/6, x : 1/12, x : 1/20, y: 1/6, y : 1/12, y : 1/20, x*y: 1/24, x*y : 1/60, x *y: 1/60}

API reference¶

sympy.integrals.intpoly.polytope_integrate(poly, expr=None, **kwargs)[source]¶

Integrates polynomials over 2/3-Polytopes.

This function accepts the polytope in \(poly\) and the function in \(expr\) (uni/bi/trivariate polynomials are implemented) and returns the exact integral of \(expr\) over \(poly\).

Parameters

poly : The input Polygon.

expr : The input polynomial.

clockwise : Binary value to sort input points of 2-Polytope clockwise.(Optional)

max_degree : The maximum degree of any monomial of the input polynomial.(Optional)

Examples

>>> from sympy.abc import x, y
>>> from sympy.geometry.polygon import Polygon
>>> from sympy.geometry.point import Point
>>> from sympy.integrals.intpoly import polytope_integrate
>>> polygon = Polygon(Point(0, 0), Point(0, 1), Point(1, 1), Point(1, 0))
>>> polys = [1, x, y, x*y, x**2*y, x*y**2]
>>> expr = x*y
>>> polytope_integrate(polygon, expr)
1/4
>>> polytope_integrate(polygon, polys, max_degree=3)
{1: 1, x: 1/2, y: 1/2, x*y: 1/4, x*y**2: 1/6, x**2*y: 1/6}

Symbolic Integrals¶

Examples¶

Integral Transforms¶

Internals¶

API reference¶

TODO and Bugs¶

Numeric Integrals¶

Integration over Polytopes¶

References¶

Examples¶

For 2D Polygons¶

For 3-Polytopes/Polyhedra¶

API reference¶

Table of Contents

Previous topic

Next topic

This Page