$\newcommand{\dee}{\mathop{\mathrm{d}\!}}$

Integral Calculus in Python

This post shows how to perform integral calculus of continuous and limited real functions of real variables in Python through the use of common Python libraries frequently used in scientific applications. The integral calculation techniques here are both primarily numerical since this site deals with computation, however some analytical techniques are also shown.
The post is organized by examples: each paragraph contains an example of an integral to compute and the related Python code snippet that calculates it using an appropriate library.

All of the various code snippets described in this post require Python version 3 and the NumPy library, while individually they require an additional library (and its dependencies, if any) between SciPy and SymPy.

We thank Prof. Fausta D'Acunzo of Preparazione 2.0 for theoretical support provided on multi-variable integral calculus.

To get the code see the Full code download paragraph at the bottom of this post.

Integration via SciPy

Integral of function of one variable (with finite extremes)

In integral calculus, the definite integral is an operator that, given a real-valued function of a real-valued variable and an interval $[a,b]$ (subset of the domain), associates to the function the area subtended by its graph in the interval $[a,b]$.
The SciPy library provides several numerical methods for computing the integral of such functions; the purpose of this chapter is to present a series of demos to show "by examples" the use of such methods; for complete documentation of these methods, the reader is invited to consult the official documentation of SciPy.

Let the following integral of a function of one variable be given: $$\int_{1}^{5} 2 x e^{-x} \,dx$$ whose analytical solution is $\approx 1.3907$ verifiable online via Wolfram Alpha.

Below is the example of Python code that calculates the integral using the quad function of the SciPy library:

import scipy.integrate as spi
import numpy as np

print('Single integral computed by SciPy quad')
print('Integral of 2xe^-x from x=1 to x=5')

integrand = lambda x : 2 * x * np.exp(-x)
a = 1.
b = 5.

result, error = spi.quad(integrand, a, b)
print('Result is ', result, 'with error ', error)
whose output is:

Single integral computed by SciPy quad
Integral of 2xe^-x from x=1 to x=5
Result is  1.3906624006967436 with error  1.54394541673402e-14
Here is the link to the code on GitHub.

Below is the example of Python code that calculates the integral using the fixed_quad function of the SciPy library:

import scipy.integrate as spi
import numpy as np

print('Single integral computed by SciPy fixed_quad')
print('Integral of 2xe^-x from x=1 to x=5')

integrand = lambda x : 2 * x * np.exp(-x)
a = 1.
b = 5.

result, none = spi.fixed_quad(integrand, a, b, n=5)
print('Result is ', result)
whose output is:

Single integral computed by SciPy fixed_quad
Integral of 2xe^-x from x=1 to x=5
Result is  1.39066368239563
Here is the link to the code on GitHub.

Below is the example of Python code that calculates the integral using the quadrature function of the SciPy library:

import scipy.integrate as spi
import numpy as np

print('Single integral computed by SciPy quadrature')
print('Integral of 2xe^-x from x=1 to x=5')

integrand = lambda x : 2 * x * np.exp(-x)
a = 1.
b = 5.

result, error = spi.quadrature(integrand, a, b)
print('Result is ', result, 'with error ', error)
whose output is:

Single integral computed by SciPy quadrature
Integral of 2xe^-x from x=1 to x=5
Result is  1.3906624007789083 with error  1.225092050027854e-08
Here is the link to the code on GitHub.

Calculating the integral with romberg

Below is the example of Python code that calculates the integral using the romberg function of the SciPy library:

import scipy.integrate as spi
import numpy as np

print('Single integral computed by SciPy romberg')
print('Example 1-01 romberg')
print('Integral of 2xe^-x from x=1 to x=5')

integrand = lambda x : 2 * x * np.exp(-x)
a = 1.
b = 5.

result = spi.romberg(integrand, a, b)
print('Result is ', result)
whose output is:

Single integral computed by SciPy romberg
Example 1-01 romberg
Integral of 2xe^-x from x=1 to x=5
Result is  1.3906624006967394
Here is the link to the code on GitHub.

Calculating the integral with trapezoid

The trapezoid function is a fixed-sample function integration method, so the code first discretizes the integration interval evenly spaced and for all discrete values of $x$ it computes the corresponding values of $y$ and then passes the two sets of discrete values $x$ and $y$ to the integration method.

Below is the example of Python code that calculates the integral using the trapezoid function of the SciPy library:

import scipy.integrate as spi
import numpy as np

print('Single integral computed by SciPy trapezoid')
print('Example 1-01 trapezoid')
print('Integral of 2xe^-x from x=1 to x=5')

integrand = lambda x : 2 * x * np.exp(-x)
a = 1.
b = 5.
step = 1e-4

xs = np.arange(a, b, step)
ys = integrand(xs)

result = spi.trapezoid(ys, xs)
print('Result is ', result)
whose output is:

Single integral computed by SciPy trapezoid
Example 1-01 trapezoid
Integral of 2xe^-x from x=1 to x=5
Result is  1.3906556624352673
Here is the link to the code on GitHub.

Calculating the integral with cumulative_trapezoid

The function cumulative_trapezoid is also a fixed-sample function integration method, and so what was said about trapezoid applies.
Below is the example of Python code that calculates the integral using of the cumulative_trapezoid function of the SciPy library:

import scipy.integrate as spi
import numpy as np

print('Single integral computed by SciPy cumulative_trapezoid')
print('Example 1-01 cumulative_trapezoid')
print('Integral of 2xe^-x from x=1 to x=5')

integrand = lambda x : 2 * x * np.exp(-x)
a = 1.
b = 5.
step = 1e-4

xs = np.arange(a, b, step)
ys = integrand(xs)

result = spi.cumulative_trapezoid(ys, xs)
result = result[-1]
print('Result is ', result)
whose output is:

Single integral computed by SciPy cumulative_trapezoid
Example 1-01 cumulative_trapezoid
Integral of 2xe^-x from x=1 to x=5
Result is  1.3906556624352677
Here is the link to the code on GitHub.

Calculating the integral with simpson

The function simpson is also a fixed-sample function integration method, and so what was said about trapezoid applies.
Below is the example of Python code that calculates the integral using the simpson function of the SciPy library:

import scipy.integrate as spi
import numpy as np

print('Single integral computed by SciPy simpson')
print('Example 1-01 simpson')
print('Integral of 2xe^-x from x=1 to x=5')

integrand = lambda x : 2 * x * np.exp(-x)
a = 1.
b = 5.
step = 1e-4

xs = np.arange(a, b, step)
ys = integrand(xs)

result = spi.simpson(ys, xs)
print('Result is ', result)
whose output is:

Single integral computed by SciPy simpson
Example 1-01 simpson
Integral of 2xe^-x from x=1 to x=5
Result is  1.3906556624801614
Here is the link to the code on GitHub.

Integral of function of one variable with quad (with extreme infinity)

Set the following integral of a function of one variable be given: $$\int_{1}^{+\infty} 2 x e^{-x} \,dx$$ whose analytical solution is $\approx 1.4715$ verifiable online via Wolfram Alpha.

Below is the example of Python code that calculates the integral using the quad function of the SciPy library:

import scipy.integrate as spi
import numpy as np

print('Single integral computed by SciPy quad')
print('Integral of 2xe^-x from x=1 to x-->+inf')

integrand = lambda x : 2 * x * np.exp(-x)
a = 1.
b = np.inf

result, error = spi.quad(integrand, a, b)
print('Result is ', result, 'with error ', error)
whose output is:

Single integral computed by SciPy quad
Integral of 2xe^-x from x=1 to x-->+inf
Result is  1.4715177646857691 with error  3.7568301883294814e-10
Here is the link to the code on GitHub.

Calculation of the length of a planar curve arc

As is well known, the integrals of a function of one variable can also be used to calculate the length of an arc of a planar curve. If the curve is expressed in the form $y=f(x)$ and is continuous and derivable the formula for calculating the length of the;arc of the curve between $x=a$ and $x=b$ is as follows: $$\int_{a}^{b} \sqrt{1 + \left(\frac{\,dy}{\,dx}\right)^2} \,dx$$
If instead curve is expressed in the parametric form $x=f_x(t)$ and $y=f_y(t)$ with both $f_x$ and $f_y$ continuous and derivable, the formula for calculating the length of the arc of the curve between $t=a$ and $t=b$ is as follows: $$\int_{a}^{b} \sqrt{\left(\frac{\,dx}{\,dt}\right)^2 + \left(\frac{\,dy}{\,dt}\right)^2} \,dt$$
Note: in the two examples that follow, the prime derivatives are calculated using the AutoGrad library.

Calculation of the length of a planar curve arc expressed in explicit form with quad

Let the following planar curve be given in explicit form: $$y=e^{-x^2}$$ by applying the corresponding formula above to calculate the length of the arc between $x=-1$ and $x = 1$, the integral is obtained: $$\int_{-1}^{1} \sqrt{1 + \left(\frac{\,d(e^{-x^2})}{\,dx}\right)^2} \,dx$$ whose analytical solution is $\approx 2.4089$ verifiable online via Wolfram Alpha.

Here is the example of Python code that calculates length of a planar curve arc expressed in explicit form using the quad function of SciPy library:

import scipy.integrate as spi
import autograd.numpy as anp   # Thinly-wrapped version of Numpy

print('Compute length of an arc of planar curve by SciPy quad')
print('Length of arc of curve y=e^(-x^2) from x=-1 to x=1')

y = lambda x : anp.exp(-x**2)
integrand = lambda x : anp.sqrt(1 + dy_dx(x) ** 2)
a = -1.
b = 1.

result, error = spi.quad(integrand, a, b)
print('Result is ', result, 'with error ', error)
whose output is:

Compute length of an arc of planar curve by SciPy quad
Length of arc of curve y=e^(-x^2) from x=-1 to x=1
Result is  2.408882141747104 with error  1.3501578327579367e-12
Here is the link to the code on GitHub.

Calculation of the length of a planar curve arc expressed in parametric form with quad

Let the following planar curve be given in parametric form: $$x(t)=cos^3 t$$ $$y(t)=sin^3 t$$ by applying the corresponding formula above to calculate the length of the arc between $t=0$ e $t = 2\pi$, the integral is obtained: $$\int_{0}^{2 \pi} \sqrt{\left(\frac{\,d (cos^3 t)}{\,dt}\right)^2 + \left(\frac{\,d (sin^3 t)}{\,dt}\right)^2} \,dt$$ whose analytical solution is $6$ verifiable online via Wolfram Alpha.

Here is the example of Python code that calculates length of a planar curve arc expressed in parametric form using the quad function of SciPy library:

import scipy.integrate as spi
import autograd.numpy as anp   # Thinly-wrapped version of Numpy

print('Compute length of an arc of planar parametric curve by SciPy quad')
print('Length of arc of parametric curve x(t)=cos^3(t) and y(t)=sin^3(t) from t=0 to t=2pi')

x = lambda t : anp.cos(t) ** 3
y = lambda t : anp.sin(t) ** 3

integrand = lambda t : anp.sqrt(dx_dt(t) ** 2 + dy_dt(t) ** 2)
a = 0.
b = 2 * anp.pi

result, error = spi.quad(integrand, a, b)
print('Result is ', result, 'with error ', error)
whose output is:

Compute length of an arc of planar parametric curve by SciPy quad
Length of arc of parametric curve x(t)=cos^3(t) and y(t)=sin^3(t) from t=0 to t=2pi
Result is  6.0 with error  6.616929226765933e-14
Here is the link to the code on GitHub.

Double integral of a function of two variables

In integral calculus, the definite double integral is an operator that, given a real-valued function of two real-valued variables and a set included in the domain, associates to the function the volume of the solid (called cylindroid) between the surface described by the function and the plane containing the given set.

Let the following double integral of a function of two variables be given: $$\int_{1}^{5} \int_{y-1}^{y+1} 2 x y e^{-x y} \,dx dy$$ whose analytical solution is $\approx 1.0273$ verifiable online via Wolfram Alpha.

Below is the example of Python code that calculates the integral using the dblquad function of the SciPy library:

import scipy.integrate as spi
import numpy as np

print('Double integral computed by SciPy dblquad')
print('Integral of 2xye^-xy from y=1 to y=5 and from x=y-1 to x=y+1')

integrand = lambda x, y : 2 * x * y * np.exp(-x * y)
ya = 1.
yb = 5.

result, error = spi.dblquad(integrand, ya, yb, lambda y: y-1, lambda y: y+1)
print('Result is ', result, 'with error ', error)
whose output is:

Double integral computed by SciPy dblquad
Integral of xye^-xy from y=1 to y=5 and from x=y-1 to x=y+1
Result is  1.0273038469260316 with error  1.3420097685964054e-14
Here is the link to the code on GitHub.

Below is the example of Python code that calculates the integral using the nquad function of the SciPy library:

import scipy.integrate as spi
import numpy as np

print('Double integral computed by SciPy nquad')
print('Integral of 2xye^-xy from y=1 to y=5 and from x=y-1 to x=y+1')

integrand = lambda x, y : 2 * x * y * np.exp(-x * y)

bounds_y = lambda : [1., 5.]
bounds_x = lambda y : [y-1, y+1]

result, error = spi.nquad(integrand, [bounds_x, bounds_y])
print('Result is ', result, 'with error ', error)
whose output is:

Double integral computed by SciPy nquad
Integral of 2xye^-xy from y=1 to y=5 and from x=y-1 to x=y+1
Result is  1.0273038469260316 with error  1.3420097685964054e-14
Here is the link to the code on GitHub.

Double integral of a function of two variables with nquad (other example)

Let the following double integral of a function of two variables be given: $$\int_{1}^{+\infty} \int_{1}^{+\infty} 2 x y e^{-x y} \,dx dy$$ whose analytical solution is $\approx 1.17453$ verifiable online via Wolfram Alpha.

Below is the example of Python code that calculates the integral using the nquad function of the SciPy library:

import scipy.integrate as spi
import numpy as np

print('Double integral computed by SciPy nquad')
print('Integral of 2xye^-xy from y=1 to x-->+inf and from x=1 to x-->+inf')

integrand = lambda x, y : 2 * x * y * np.exp(-x * y)
ya = 1.
yb = 5.

result, error = spi.nquad(integrand, [[1, np.inf],[1, np.inf]])
print('Result is ', result, 'with error ', error)
whose output is:

Double integral computed by SciPy nquad
Integral of 2xye^-xy from y=1 to y-->+inf and from x=1 to x-->+inf
Result is  1.1745267511339024 with error  1.6321550842479074e-08
Here is the link to the code on GitHub.

Triple Integral of a function of three variables

Let the following triple integral of a function of three variables be given: $$\int_{1}^{2} \int_{z+1}^{z+2} \int_{y+z}^{2(y+z)} x + yz^2 \,dx dy dz$$ whose analytical solution is $\approx 65.7194$ verifiable online via Wolfram Alpha.

Below is the example of Python code that calculates the integral using the tplquad function of the SciPy library:

import scipy.integrate as spi
import numpy as np

print('Triple integral computed by SciPy tplquad')
print('Integral of x + yz^2 from z=1 to z=2, y=z+1 to y=z+2 and from x=y+x to x=2(y+z)')

integrand = lambda x, y, z : x + y * z ** 2
za = 1.
zb = 2.
ya=lambda z: z + 1
yb=lambda z: z + 2
xa=lambda z, y : y + z
xb=lambda z, y : 2 * (y + z)

result, error = spi.tplquad(integrand, za, zb, ya, yb, xa, xb)
print('Result is ', result, 'with error ', error)
whose output is:

Triple integral computed by SciPy tplquad
Integral of x + yz^2 from z=1 to z=2, y=z+1 to y=z+2 and from x=y+x to x=2(y+z)
Result is  65.71944444444445 with error  1.659412309590769e-12
Here is the link to the code on GitHub.

Below is the example of Python code that calculates the integral using the nquad function of the SciPy library:

import scipy.integrate as spi
import numpy as np

print('Triple integral computed by SciPy nquad')
print('Integral of x + yz^2 from z=1 to z=2, y=z+1 to y=z+2 and from x=y+x to x=2(y+z)')

integrand = lambda x, y, z : x + y * z ** 2

bounds_z = lambda : [1., 2.]
bounds_y = lambda z : [z+1, z+2]
bounds_x = lambda z, y : [y+z, 2 * (y+z)]
ya=lambda z: z + 1
yb=lambda z: z + 2
xa=lambda z, y : y + z
xb=lambda z, y : 2 * (y + z)

result, error = spi.nquad(integrand, [bounds_x, bounds_y, bounds_z])
print('Result is ', result, 'with error ', error)
whose output is:

Triple integral computed by SciPy nquad
Integral of x + yz^2 from z=1 to z=2, y=z+1 to y=z+2 and from x=y+x to x=2(y+z)
Result is  65.71944444444445 with error  1.659412309590769e-12
Here is the link to the code on GitHub.

Integration via SymPy

Integral of function of one variable (with finite extremes)

In integral calculus, the definite integral is an operator that, given a real-valued function of a real-valued variable and an interval $[a,b]$ (subset of the domain), associates to the function the area subtended by its graph in the interval $[a,b]$.
The SciPy library provides several numerical methods for computing the integral of such functions; the purpose of this chapter is to present a series of demos to show "by examples" the use of such methods; for complete documentation of these methods, the reader is invited to consult the official documentation of SciPy.

Let the following integral of a function of one variable be given: $$\int_{1}^{5} 2 x e^{-x} \,dx$$ whose analytical solution is $\approx 1.3907$ verifiable online via Wolfram Alpha.

Integral of a function of one variable with integrate(f, (x, a, b))

Below is the example of Python code that calculates the integral using the integrate(f, (x, a, b)) function of the SymPy library:

import sympy as sp

print('Single integral computed by SymPy definite integrate')
print('Example 1-01 definite')
print('Integral of 2xe^-x from x=1 to x=5')

x = sp.Symbol('x')
f = 2 * x * sp.exp(-x)
a = 1.
b = 5.
integral = sp.integrate(f, (x, a, b))
integral = integral.evalf()
print('Result is ', integral)
whose output is:

Single integral computed by SymPy definite integrate
Example 1-01 definite
Integral of 2xe^-x from x=1 to x=5
Result is  1.39066240069674
Here is the link to the code on GitHub.

Integral of a function of one variable with integrate(f, (x, a, b)) via indefinite integral

Below is the example of Python code that calculates the integral using the integrate(f, x) function of the SymPy library:

import sympy as sp

print('Single integral computed by SymPy indefinite integrate')
print('Example 1-01 indefinite integrate')
print('Integral of 2xe^-x from x=1 to x=5')

x = sp.Symbol('x')
f = 2 * x * sp.exp(-x)
primitive = sp.integrate(f, x)
print('Primitive is ', primitive)

primitive_lambda = sp.lambdify(x, primitive)
a = 1.
b = 5.
integral = primitive_lambda(b) - primitive_lambda(a)
print('Result is ', integral)
whose output is:

Single integral computed by SymPy indefinite integrate
Example 1-01 indefinite integrate
Integral of 2xe^-x from x=1 to x=5
Primitive is  (-2*x - 2)*exp(-x)
Result is  1.3906624006967436
The program first calculates the indefinite integral and then applying the fundamental theorem of integral calculus, calculates the value of the integral.

Here is the link to the code on GitHub.

Double integral of a function of two variables

In integral calculus, the definite double integral is an operator that, given a real-valued function of two real-valued variables and a set included in the domain, associates to the function the volume of the solid (called cylindroid) between the surface described by the function and the plane containing the given set.

Let the following double integral of a function of two variables be given: $$\int_{1}^{4} \int_{y-1}^{y+2} x y e^{-x} e^{-y} \,dx dy$$ whose analytical solution is $\approx 0.396134$ verifiable online via Wolfram Alpha.

Double Integral of a function of two variables with integrate(f, (x, xa, xb), (y, ya, yb))

Below is the example of Python code that calculates the integral using the integrate(f, (x, xa, xb), (y, ya, yb)) of the SymPy library:

import sympy as sp

print('Double integral computed by SymPy definite integrate')
print('Example 2-03 definite')
print('Integral of xye^-xy from y=1 to y=4 and from x=y-1 to x=y+1')

x, y = sp.symbols('x y')
f = x * y * sp.exp(-x) * sp.exp(-y)
ya = 0.
yb = 4.
xa = y-1.
xb = y+1.

integral = sp.integrate(f, (x, xa, xb), (y, ya, yb))
integral = integral.evalf()
print('Result is ', integral)
whose output is:

Double integral computed by SymPy definite integrate
Example 2-03 definite
Integral of xye^-xy from y=1 to y=4 and from x=y-1 to x=y+1
Result is  0.396134380699524
Here is the link to the code on GitHub.