**Computational Mindset**by Ettore Messina

# 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)
```

whose output is:```
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
```

Here is the link to the code on GitHub.
#### 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)
```

whose output is:```
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
```

Here is the link to the code on GitHub.
#### 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)
```

whose output is:```
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
```

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('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)
```

whose output is:```
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
```

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
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)
```

whose output is:```
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
```

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
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)
```

whose output is:```
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
```

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.

#### 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)
```

whose output is:```
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
```

Here is the link to the code on GitHub.
#### 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)
```

whose output is:```
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
```

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('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)
```

whose output is:```
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
```

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.

#### 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)
```

whose output is:```
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
```

Here is the link to the code on GitHub.
#### 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)
```

whose output is:```
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
```

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.
## 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.