Integrals¶
The integrals
module in Diofant implements methods to calculate definite and indefinite integrals of expressions.
Principal method in this module is integrate()
integrate(f, x)
returns the indefinite integral \(\int f\,dx\)integrate(f, (x, a, b))
returns the definite integral \(\int_{a}^{b} f\,dx\)
Examples¶
Diofant can integrate a vast array of functions. It can integrate polynomial functions:
>>> init_printing(pretty_print=True, use_unicode=False, wrap_line=False, no_global=True)
>>> 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
Exponentialpolynomial 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, Diofant will deal with these integrals.
>>> integrate(x**2 * exp(x) * cos(x), x)
x 2 x 2 x x
E *x *sin(x) E *x *cos(x) x E *sin(x) E *cos(x)
 +   E *x*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¶
Diofant has special support for definite integrals, and integral transforms.

diofant.integrals.transforms.
mellin_transform
(f, x, s, **hints)[source]¶ Compute the Mellin transform \(F(s)\) of \(f(x)\),
\[F(s) = \int_0^\infty x^{s1} 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)
whereF
is the Mellin transform off
,(a, b)
is the fundamental strip (as above), andcond
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
diofant.integrals.transforms.IntegralTransform.doit()
. Ifnoconds=False
, then only \(F\) will be returned (i.e. notcond
, and also not the strip(a, b)
).>>> from diofant.abc import s >>> mellin_transform(exp(x), x, s) (gamma(s), (0, oo), True)

diofant.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) = \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 diofant assumptions!
For a description of possible hints, refer to the docstring of
diofant.integrals.transforms.IntegralTransform.doit()
.>>> from diofant.abc import s >>> inverse_mellin_transform(gamma(s), s, x, (0, oo)) E**(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(x + 1)/2  Heaviside(x  1)/(2*x) >>> inverse_mellin_transform(f, s, x, (1, oo)) (x/2 + 1/(2*x))*Heaviside(x + 1)

diofant.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)
whereF
is the Laplace transform off
, \(\operatorname{Re}(s) > a\) is the halfplane of convergence, andcond
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
diofant.integrals.transforms.IntegralTransform.doit()
. Ifnoconds=True
, only \(F\) will be returned (i.e. notcond
, and also not the planea
).>>> from diofant.abc import t, s, a >>> laplace_transform(t**a, t, s) (s**(a)*gamma(a + 1)/s, 0, re(a) < 1)

diofant.integrals.transforms.
inverse_laplace_transform
(F, s, t, plane=None, **hints)[source]¶ Compute the inverse Laplace transform of \(F(s)\), defined as
\[f(t) = \int_{ci\infty}^{c+i\infty} e^{st} F(s) \mathrm{d}s,\]for \(c\) so large that \(F(s)\) has no singularites in the halfplane \(\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 nonnegative \(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 diofant assumption on \(t\).
For a description of possible hints, refer to the docstring of
diofant.integrals.transforms.IntegralTransform.doit()
.>>> from diofant.abc import s, t >>> a = Symbol('a', positive=True) >>> inverse_laplace_transform(exp(a*s)/s, s, t) Heaviside(a + t)

diofant.integrals.transforms.
fourier_transform
(f, x, k, **hints)[source]¶ Compute the unitary, ordinaryfrequency 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 a description of possible hints, refer to the docstring of
diofant.integrals.transforms.IntegralTransform.doit()
. Note that for this transform, by defaultnoconds=True
.>>> fourier_transform(exp(x**2), x, k) E**(pi**2*k**2)*sqrt(pi) >>> fourier_transform(exp(x**2), x, k, noconds=False) (E**(pi**2*k**2)*sqrt(pi), True)

diofant.integrals.transforms.
inverse_fourier_transform
(F, k, x, **hints)[source]¶ Compute the unitary, ordinaryfrequency 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 a description of possible hints, refer to the docstring of
diofant.integrals.transforms.IntegralTransform.doit()
. Note that for this transform, by defaultnoconds=True
.>>> inverse_fourier_transform(sqrt(pi)*exp((pi*k)**2), k, x) E**(x**2) >>> inverse_fourier_transform(sqrt(pi)*exp((pi*k)**2), k, x, noconds=False) (E**(x**2), True)

diofant.integrals.transforms.
sine_transform
(f, x, k, **hints)[source]¶ Compute the unitary, ordinaryfrequency 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
diofant.integrals.transforms.IntegralTransform.doit()
. Note that for this transform, by defaultnoconds=True
.>>> from diofant.abc import a >>> sine_transform(x*exp(a*x**2), x, k) sqrt(2)*E**(k**2/(4*a))*k/(4*a**(3/2)) >>> sine_transform(x**(a), x, k) 2**(a + 1/2)*k**(a  1)*gamma(a/2 + 1)/gamma(a/2 + 1/2)

diofant.integrals.transforms.
inverse_sine_transform
(F, k, x, **hints)[source]¶ Compute the unitary, ordinaryfrequency 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
diofant.integrals.transforms.IntegralTransform.doit()
. Note that for this transform, by defaultnoconds=True
.>>> from diofant.abc import a >>> inverse_sine_transform(2**((12*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) E**(a*x**2)*x

diofant.integrals.transforms.
cosine_transform
(f, x, k, **hints)[source]¶ Compute the unitary, ordinaryfrequency 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
diofant.integrals.transforms.IntegralTransform.doit()
. Note that for this transform, by defaultnoconds=True
.>>> from diofant.abc import 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) E**(a**2/(2*k))*a/(2*k**(3/2))

diofant.integrals.transforms.
inverse_cosine_transform
(F, k, x, **hints)[source]¶ Compute the unitary, ordinaryfrequency 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
diofant.integrals.transforms.IntegralTransform.doit()
. Note that for this transform, by defaultnoconds=True
.>>> from diofant.abc import a >>> inverse_cosine_transform(sqrt(2)*a/(sqrt(pi)*(a**2 + k**2)), k, x) E**(a*x) >>> inverse_cosine_transform(1/sqrt(k), k, x) 1/sqrt(x)

diofant.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
diofant.integrals.transforms.IntegralTransform.doit()
. Note that for this transform, by defaultnoconds=True
.>>> from diofant.abc import r, nu, a, k
>>> 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) E**(a*r)

diofant.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
diofant.integrals.transforms.IntegralTransform.doit()
. Note that for this transform, by defaultnoconds=True
.>>> from diofant.abc import r, nu, a, k
>>> 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) E**(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.
Diofant currently uses a simplified version of the Risch algorithm, called the RischNorman algorithm. This algorithm is much faster, but may fail to find an antiderivative, although it is still very powerful. Diofant also uses pattern matching and heuristics to speed up evaluation of some types of integrals, e.g. polynomials.
For nonelementary definite integrals, Diofant uses socalled Meijer Gfunctions. Details are described here.
API reference¶

diofant.integrals.integrals.
integrate
(f, var, ...)[source]¶ Compute definite or indefinite integral of one or more variables using RischNorman 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
Diofant 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.
Diofant 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 nonelementary. 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 rewriting the integrand in terms of socalled Meijer Gfunctions. Indefinite integrals of a single Gfunction can always be computed, and the definite integral of a product of two Gfunctions can be computed from zero to infinity. Various strategies are implemented to rewrite integrands as Gfunctions, and use this information to compute integrals (see the
meijerint
module).In general, the algebraic methods work best for computing antiderivatives of (possibly complicated) combinations of elementary functions. The Gfunction 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 Gfunction 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 Gfunctions second to last, and heuristic Risch last).
 If still not successful, try Gfunctions irrespective of the limits.
The option meijerg=True, False, None can be used to, respectively: always use Gfunction methods and no others, never use Gfunction methods, or use all available methods (in order as described above). It defaults to None.
Examples
>>> from diofant.abc import a
>>> 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:
>>> 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(E**(x)*x**a, (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)

diofant.integrals.integrals.
line_integrate
(field, Curve, variables)[source]¶ Compute the line integral.
Examples
>>> from diofant.abc import t >>> C = Curve([E**t + 1, E**t  1], (t, 0, ln(2))) >>> line_integrate(x + y, C, [x, y]) 3*sqrt(2)

diofant.integrals.deltafunctions.
deltaintegrate
(f, x)[source]¶ The idea for integration is the following:
If we are dealing with a DiracDelta expression, i.e. DiracDelta(g(x)), we try to simplify it.
If we could simplify it, then we integrate the resulting expression. We already know we can integrate a simplified expression, because only simple DiracDelta expressions are involved.
If we couldn’t simplify it, there are two cases:
 The expression is a simple expression: we return the integral, taking care if we are dealing with a Derivative or with a proper DiracDelta.
 The expression is not simple (i.e. DiracDelta(cos(x))): we can do nothing at all.
If the node is a multiplication node having a DiracDelta term:
First we expand it.
If the expansion did work, then we try to integrate the expansion.
If not, we try to extract a simple DiracDelta term, then we have two cases:
 We have a simple DiracDelta term, so we return the integral.
 We didn’t have a simple term, but we do have an expression with simplified DiracDelta terms, so we integrate this expression.
Examples
>>> deltaintegrate(x*sin(x)*cos(x)*DiracDelta(x  1), x) sin(1)*cos(1)*Heaviside(x  1) >>> deltaintegrate(y**2*DiracDelta(x  z)*DiracDelta(y  z), y) z**2*DiracDelta(x  z)*Heaviside(y  z)

diofant.integrals.rationaltools.
ratint
(f, x, **flags)[source]¶ Performs indefinite integration of rational functions.
Given a field \(K\) and a rational function \(f = p/q\), where \(p\) and \(q\) are polynomials in \(K[x]\), returns a function \(g\) such that \(f = g'\).
>>> ratint(36/(x**5  2*x**4  2*x**3 + 4*x**2 + x  2), x) (12*x + 6)/(x**2  1) + 4*log(x  2)  4*log(x + 1)
See also
diofant.integrals.integrals.Integral.doit
,diofant.integrals.rationaltools.ratint_logpart
,diofant.integrals.rationaltools.ratint_ratpart
References
[Bro05388] M. Bronstein, Symbolic Integration I: Transcendental Functions, Second Edition, SpringerVerlag, 2005, pp. 3570

diofant.integrals.rationaltools.
ratint_logpart
(f, g, x, t=None)[source]¶ LazardRiobooTrager algorithm.
Given a field K and polynomials f and g in K[x], such that f and g are coprime, deg(f) < deg(g) and g is squarefree, returns a list of tuples (s_i, q_i) of polynomials, for i = 1..n, such that s_i in K[t, x] and q_i in K[t], and:
___ ___ d f d \ ` \ `   =  ) ) a log(s_i(a, x)) dx g dx /__, /__, i=1..n a  q_i(a) = 0
Examples
>>> ratint_logpart(Poly(1, x), Poly(x**2 + x + 1, x), x) [(Poly(x + 3*_t/2 + 1/2, x, domain='QQ[_t]'), Poly(3*_t**2 + 1, _t, domain='ZZ'))] >>> ratint_logpart(Poly(12, x), Poly(x**2  x  2, x), x) [(Poly(x  3*_t/8  1/2, x, domain='QQ[_t]'), Poly(_t**2 + 16, _t, domain='ZZ'))]

diofant.integrals.rationaltools.
ratint_ratpart
(f, g, x)[source]¶ HorowitzOstrogradsky algorithm.
Given a field K and polynomials f and g in K[x], such that f and g are coprime and deg(f) < deg(g), returns fractions A and B in K(x), such that f/g = A’ + B and B has squarefree denominator.
Examples
>>> ratint_ratpart(Poly(1, x), Poly(x + 1, x), x) (0, 1/(x + 1)) >>> ratint_ratpart(Poly(1, x, domain='EX'), ... Poly(x**2 + y**2, x, domain='EX'), x) (0, 1/(x**2 + y**2)) >>> ratint_ratpart(Poly(36, x), ... Poly(x**5  2*x**4  2*x**3 + 4*x**2 + x  2, x), x) ((12*x + 6)/(x**2  1), 12/(x**2  x  2))

diofant.integrals.heurisch.
components
(f, x)[source]¶ Returns a set of all functional components of the given expression which includes symbols, function applications and compositions and noninteger powers. Fractional powers are collected with with minimal, positive exponents.
>>> components(sin(x)*cos(x)**2, x) {x, sin(x), cos(x)}
See also

diofant.integrals.heurisch.
heurisch
(f, x, rewrite=False, hints=None, mappings=None, retries=3, degree_offset=0, unnecessary_permutations=None)[source]¶ Compute indefinite integral using heuristic Risch algorithm.
This is a heuristic approach to indefinite integration in finite terms using the extended heuristic (parallel) Risch algorithm, based on Manuel Bronstein’s “Poor Man’s Integrator” [R390].
The algorithm supports various classes of functions including transcendental elementary or special functions like Airy, Bessel, Whittaker and Lambert.
Note that this algorithm is not a decision procedure. If it isn’t able to compute the antiderivative for a given function, then this is not a proof that such a functions does not exist. One should use recursive Risch algorithm in such case. It’s an open question if this algorithm can be made a full decision procedure.
This is an internal integrator procedure. You should use toplevel ‘integrate’ function in most cases, as this procedure needs some preprocessing steps and otherwise may fail.
Parameters: f : Expr
expression
x : Symbol
variable
rewrite : Boolean, optional
force rewrite ‘f’ in terms of ‘tan’ and ‘tanh’, default False.
hints : None or list
a list of functions that may appear in antiderivate. If None (default)  no suggestions at all, if empty list  try to figure out.
See also
diofant.integrals.integrals.Integral.doit
,diofant.integrals.integrals.Integral
,diofant.integrals.heurisch.components
References
[R390] (1, 2) Manuel Bronstein’s “Poor Man’s Integrator”, http://wwwsop.inria.fr/cafe/Manuel.Bronstein/pmint/index.html [R391] K. Geddes, L. Stefanus, On the RischNorman Integration Method and its Implementation in Maple, Proceedings of ISSAC‘89, ACM Press, 212217. [R392] J. H. Davenport, On the Parallel Risch Algorithm (I), Proceedings of EUROCAM‘82, LNCS 144, Springer, 144157. [R393] J. H. Davenport, On the Parallel Risch Algorithm (III): Use of Tangents, SIGSAM Bulletin 16 (1982), 36. [R394] J. H. Davenport, B. M. Trager, On the Parallel Risch Algorithm (II), ACM Transactions on Mathematical Software 11 (1985), 356362. Examples
>>> heurisch(y*tan(x), x) y*log(tan(x)**2 + 1)/2

diofant.integrals.heurisch.
heurisch_wrapper
(f, x, rewrite=False, hints=None, mappings=None, retries=3, degree_offset=0, unnecessary_permutations=None)[source]¶ A wrapper around the heurisch integration algorithm.
This method takes the result from heurisch and checks for poles in the denominator. For each of these poles, the integral is reevaluated, and the final integration result is given in terms of a Piecewise.
See also
Examples
>>> heurisch(cos(n*x), x) sin(n*x)/n >>> heurisch_wrapper(cos(n*x), x) Piecewise((x, Eq(n, 0)), (sin(n*x)/n, true))

diofant.integrals.trigonometry.
trigintegrate
(f, x, conds='piecewise')[source]¶ Integrate f = Mul(trig) over x
>>> trigintegrate(sin(x)*cos(x), x) sin(x)**2/2
>>> trigintegrate(sin(x)**2, x) x/2  sin(x)*cos(x)/2
>>> trigintegrate(tan(x)*sec(x), x) 1/cos(x)
>>> trigintegrate(sin(x)*tan(x), x) log(sin(x)  1)/2 + log(sin(x) + 1)/2  sin(x)
References
[R395] https//en.wikibooks.org/wiki/Calculus/Integration_techniques
The class \(Integral\) represents an unevaluated integral and has some methods that help in the integration of an expression.

class
diofant.integrals.integrals.
Integral
[source]¶ Represents unevaluated integral.

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

as_sum
(n, method='midpoint')[source]¶ Approximates the definite integral by a sum.
method … one of: left, right, midpoint, trapezoid
These are all basically the rectangle method [1], the only difference is where the function value is taken in each interval to define the rectangle.
See also
diofant.integrals.integrals.Integral.doit
 Perform the integration using any hints
References
[R396] https//en.wikipedia.org/wiki/Rectangle_method Examples
>>> 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 lefthand 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 righthand 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
All but the trapexoid method may be used when dealing with a function with a discontinuity. Here, the discontinuity at x = 0 can be avoided by using the midpoint or righthand 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 oo:
>>> e.as_sum(5, 'left') oo >>> e.as_sum(5, 'trapezoid') oo

doit
(**hints)[source]¶ Perform the integration using any hints given.
See also
diofant.integrals.trigonometry.trigintegrate
,diofant.integrals.heurisch.heurisch
,diofant.integrals.rationaltools.ratint
diofant.integrals.integrals.Integral.as_sum
 Approximate the integral using a sum
Examples
>>> from diofant.abc import i >>> Integral(x**i, (i, 1, 3)).doit() Piecewise((2, Eq(log(x), 0)), (x**3/log(x)  x/log(x), true))

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.
See also
diofant.concrete.expr_with_limits.ExprWithLimits.function
,diofant.concrete.expr_with_limits.ExprWithLimits.limits
,diofant.concrete.expr_with_limits.ExprWithLimits.variables
Examples
>>> Integral(x, (x, y, 1)).free_symbols {y}

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 usubstitution.
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.
See also
diofant.concrete.expr_with_limits.ExprWithLimits.variables
 Lists the integration variables
diofant.concrete.expr_with_limits.ExprWithLimits.as_dummy
 Replace integration variables with dummy ones
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 diofant.abc import a, b, c, d, u
>>> 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 usubstitution 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 “usubstitution”:
>>> 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, a + 1)) >>> Integral(x, (x, 0, 1)).transform(x, (u + a, a)) Integral(a + u, (a, u, u + 1))


class
diofant.integrals.transforms.
IntegralTransform
[source]¶ Base class for integral transforms.
This class represents unevaluated transforms.
To implement a concrete transform, derive from this class and implement the
_compute_transform(f, x, s, **hints)
and_as_integral(f, x, s)
functions. If the transform cannot be computed, raise IntegralTransformError.Also set
cls._name
.Implement self._collapse_extra if your function returns more than just a number and possibly a convergence condition.

doit
(**hints)[source]¶ Try to evaluate the transform in closed form.
This general function handles linearity, but apart from that leaves pretty much everything to _compute_transform.
Standard hints are the following:
simplify
: whether or not to simplify the resultnoconds
: if True, don’t return convergence conditionsneedeval
: if True, raise IntegralTransformError instead of returning IntegralTransform objects
The default values of these hints depend on the concrete transform, usually the default is
(simplify, noconds, needeval) = (True, False, False)
.

free_symbols
¶ This method returns the symbols that will exist when the transform is evaluated.

function
¶ The function to be transformed.

function_variable
¶ The dependent variable of the function to be transformed.

transform_variable
¶ The independent transform variable.


class
diofant.integrals.transforms.
MellinTransform
[source]¶ Class representing unevaluated Mellin transforms.
See also

class
diofant.integrals.transforms.
InverseMellinTransform
[source]¶ Class representing unevaluated inverse Mellin transforms.
See also

class
diofant.integrals.transforms.
LaplaceTransform
[source]¶ Class representing unevaluated Laplace transforms.
See also

class
diofant.integrals.transforms.
InverseLaplaceTransform
[source]¶ Class representing unevaluated inverse Laplace transforms.
See also

class
diofant.integrals.transforms.
FourierTransform
[source]¶ Class representing unevaluated Fourier transforms.
See also

class
diofant.integrals.transforms.
InverseFourierTransform
[source]¶ Class representing unevaluated inverse Fourier transforms.
See also

class
diofant.integrals.transforms.
SineTransform
[source]¶ Class representing unevaluated sine transforms.
See also

class
diofant.integrals.transforms.
InverseSineTransform
[source]¶ Class representing unevaluated inverse sine transforms.
See also

class
diofant.integrals.transforms.
CosineTransform
[source]¶ Class representing unevaluated cosine transforms.
See also

class
diofant.integrals.transforms.
InverseCosineTransform
[source]¶ Class representing unevaluated inverse cosine transforms.
See also

class
diofant.integrals.transforms.
HankelTransform
[source]¶ Class representing unevaluated Hankel transforms.
See also
Numeric Integrals¶
Diofant has functions to calculate points and weights for Gaussian quadrature of any order and any precision:

diofant.integrals.quadrature.
gauss_legendre
(n, n_digits)[source]¶ Computes the GaussLegendre quadrature [R397] points and weights.
The GaussLegendre 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(1x_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
andw
are lists of points and weights as Floats.The points \(x_i\) and weights \(w_i\) are returned as
(x, w)
tuple of lists.See also
gauss_laguerre
,gauss_gen_laguerre
,gauss_hermite
,gauss_chebyshev_t
,gauss_chebyshev_u
,gauss_jacobi
References
[R397] (1, 2) https//en.wikipedia.org/wiki/Gaussian_quadrature [R398] http://people.sc.fsu.edu/~jburkardt/cpp_src/legendre_rule/legendre_rule.html Examples
>>> 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]

diofant.integrals.quadrature.
gauss_laguerre
(n, n_digits)[source]¶ Computes the GaussLaguerre quadrature [R399] points and weights.
The GaussLaguerre 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
andw
are lists of points and weights as Floats.The points \(x_i\) and weights \(w_i\) are returned as
(x, w)
tuple of lists.See also
gauss_legendre
,gauss_gen_laguerre
,gauss_hermite
,gauss_chebyshev_t
,gauss_chebyshev_u
,gauss_jacobi
References
[R399] (1, 2) https//en.wikipedia.org/wiki/Gauss%E2%80%93Laguerre_quadrature [R400] http://people.sc.fsu.edu/~jburkardt/cpp_src/laguerre_rule/laguerre_rule.html Examples
>>> 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.9855e7]

diofant.integrals.quadrature.
gauss_hermite
(n, n_digits)[source]¶ Computes the GaussHermite quadrature [R401] points and weights.
The GaussHermite 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^{n1} n! \sqrt{\pi}}{n^2 \left(H_{n1}(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
andw
are lists of points and weights as Floats.The points \(x_i\) and weights \(w_i\) are returned as
(x, w)
tuple of lists.See also
gauss_legendre
,gauss_laguerre
,gauss_gen_laguerre
,gauss_chebyshev_t
,gauss_chebyshev_u
,gauss_jacobi
References
[R401] (1, 2) https//en.wikipedia.org/wiki/GaussHermite_Quadrature [R402] http://people.sc.fsu.edu/~jburkardt/cpp_src/hermite_rule/hermite_rule.html [R403] http://people.sc.fsu.edu/~jburkardt/cpp_src/gen_hermite_rule/gen_hermite_rule.html Examples
>>> 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]

diofant.integrals.quadrature.
gauss_gen_laguerre
(n, alpha, n_digits)[source]¶ Computes the generalized GaussLaguerre quadrature [R404] points and weights.
The generalized GaussLaguerre 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}_{n1}(x_i) L^{\alpha+1}_{n1}(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
andw
are lists of points and weights as Floats.The points \(x_i\) and weights \(w_i\) are returned as
(x, w)
tuple of lists.See also
gauss_legendre
,gauss_laguerre
,gauss_hermite
,gauss_chebyshev_t
,gauss_chebyshev_u
,gauss_jacobi
References
[R404] (1, 2) https//en.wikipedia.org/wiki/Gauss%E2%80%93Laguerre_quadrature [R405] http://people.sc.fsu.edu/~jburkardt/cpp_src/gen_laguerre_rule/gen_laguerre_rule.html Examples
>>> x, w = gauss_gen_laguerre(3, 0.5, 5) >>> x [0.19016, 1.7845, 5.5253] >>> w [1.4493, 0.31413, 0.00906]
>>> x, w = gauss_gen_laguerre(4, 1.5, 5) >>> x [0.97851, 2.9904, 6.3193, 11.712] >>> w [0.53087, 0.67721, 0.11895, 0.0023152]

diofant.integrals.quadrature.
gauss_chebyshev_t
(n, n_digits)[source]¶ Computes the GaussChebyshev quadrature [R406] points and weights of the first kind.
The GaussChebyshev quadrature of the first kind approximates the integral:
\[\int_{1}^{1} \frac{1}{\sqrt{1x^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
andw
are lists of points and weights as Floats.The points \(x_i\) and weights \(w_i\) are returned as
(x, w)
tuple of lists.See also
gauss_legendre
,gauss_laguerre
,gauss_hermite
,gauss_gen_laguerre
,gauss_chebyshev_u
,gauss_jacobi
References
[R406] (1, 2) https//en.wikipedia.org/wiki/Chebyshev%E2%80%93Gauss_quadrature [R407] http://people.sc.fsu.edu/~jburkardt/cpp_src/chebyshev1_rule/chebyshev1_rule.html Examples
>>> 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]

diofant.integrals.quadrature.
gauss_chebyshev_u
(n, n_digits)[source]¶ Computes the GaussChebyshev quadrature [R408] points and weights of the second kind.
The GaussChebyshev quadrature of the second kind approximates the integral:
\[\int_{1}^{1} \sqrt{1x^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
andw
are lists of points and weights as Floats.The points \(x_i\) and weights \(w_i\) are returned as
(x, w)
tuple of lists.See also
gauss_legendre
,gauss_laguerre
,gauss_hermite
,gauss_gen_laguerre
,gauss_chebyshev_t
,gauss_jacobi
References
[R408] (1, 2) https//en.wikipedia.org/wiki/Chebyshev%E2%80%93Gauss_quadrature [R409] http://people.sc.fsu.edu/~jburkardt/cpp_src/chebyshev2_rule/chebyshev2_rule.html Examples
>>> 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]

diofant.integrals.quadrature.
gauss_jacobi
(n, alpha, beta, n_digits)[source]¶ Computes the GaussJacobi quadrature [R410] points and weights.
The GaussJacobi quadrature of the first kind approximates the integral:
\[\int_{1}^1 (1x)^\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
andw
are lists of points and weights as Floats.The points \(x_i\) and weights \(w_i\) are returned as
(x, w)
tuple of lists.See also
gauss_legendre
,gauss_laguerre
,gauss_hermite
,gauss_gen_laguerre
,gauss_chebyshev_t
,gauss_chebyshev_u
References
[R410] (1, 2) https//en.wikipedia.org/wiki/Gauss%E2%80%93Jacobi_quadrature [R411] http://people.sc.fsu.edu/~jburkardt/cpp_src/jacobi_rule/jacobi_rule.html [R412] http://people.sc.fsu.edu/~jburkardt/cpp_src/gegenbauer_rule/gegenbauer_rule.html Examples
>>> x, w = gauss_jacobi(3, 0.5, 0.5, 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]