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.
Calculating the integral with quad
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('Example 1-01 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)
Single integral computed by SciPy quad
Example 1-01 quad
Integral of 2xe^-x from x=1 to x=5
Result is 1.3906624006967436 with error 1.54394541673402e-14
Calculating the integral with fixed_quad
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('Example 1-01 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)
Single integral computed by SciPy fixed_quad
Example 1-01 fixed_quad
Integral of 2xe^-x from x=1 to x=5
Result is 1.39066368239563
Calculating the integral with quadrature
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('Example 1-01 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)
Single integral computed by SciPy quadrature
Example 1-01 quadrature
Integral of 2xe^-x from x=1 to x=5
Result is 1.3906624007789083 with error 1.225092050027854e-08
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)
Single integral computed by SciPy romberg
Example 1-01 romberg
Integral of 2xe^-x from x=1 to x=5
Result is 1.3906624006967394
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)
Single integral computed by SciPy trapezoid
Example 1-01 trapezoid
Integral of 2xe^-x from x=1 to x=5
Result is 1.3906556624352673
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)
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
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)
Single integral computed by SciPy simpson
Example 1-01 simpson
Integral of 2xe^-x from x=1 to x=5
Result is 1.3906556624801614
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('Example 1-02 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)
Single integral computed by SciPy quad
Example 1-02 quad
Integral of 2xe^-x from x=1 to x-->+inf
Result is 1.4715177646857691 with error 3.7568301883294814e-10
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
import autograd as ag
print('Compute length of an arc of planar curve by SciPy quad')
print('Example 1-03 quad arclength')
print('Length of arc of curve y=e^(-x^2) from x=-1 to x=1')
y = lambda x : anp.exp(-x**2)
dy_dx = ag.grad(y)
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)
Compute length of an arc of planar curve by SciPy quad
Example 1-03 quad arclength
Length of arc of curve y=e^(-x^2) from x=-1 to x=1
Result is 2.408882141747104 with error 1.3501578327579367e-12
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
import autograd as ag
print('Compute length of an arc of planar parametric curve by SciPy quad')
print('Example 1-04 quad arclength')
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
dx_dt = ag.grad(x)
dy_dt = ag.grad(y)
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)
Compute length of an arc of planar parametric curve by SciPy quad
Example 1-04 quad arclength
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
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.
Calculating the integral with dblquad
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('Example 2-01 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)
Double integral computed by SciPy dblquad
Example 2-01 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
Calculating the integral with nquad
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('Example 2-01 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)
Double integral computed by SciPy nquad
Example 2-01 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
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('Example 2-02 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)
Double integral computed by SciPy nquad
Example 2-02 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
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.
Calculating the integral with tplquad
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('Example 3-01 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)
Triple integral computed by SciPy tplquad
Example 3-01 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
Calculating the integral with nquad
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('Example 3-01 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)
Triple integral computed by SciPy nquad
Example 3-01 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
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)
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
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)
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
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)
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
Download of the complete code
The complete code is available at GitHub.
These materials are distributed under MIT license; feel free to use, share, fork and adapt these materials as you see fit.
Also please feel free to submit pull-requests and bug-reports to this GitHub repository or contact me on my social media channels available on the top right corner of this page.