Calcolo Integrale in Python

Questo post mostra come eseguire il calcolo integrale di funzioni reali continue e limitate a variabili reali in Python tramite l'utilizzo di comuni librerie per Python frequentemente adoperate nelle applicazioni scientifiche. Le tecniche di calcolo integrale qui sono sia principalmente numeriche in quanto questo sito si occupa di computazione, tuttavia è mostrata anche qualche tecnica analitica.
Il post è organizzato per esempi: ogni paragrafo contiene un esempio di un integrale da calcolare e il relativo frammento di codice Python che lo calcola utilizzando una opportuna libreria.

Tutti i vari frammenti di codice descritti in questo post richiedono la versione 3 di Python e la libreria NumPy, mentre singolarmente richiedono una ulteriore libreria addizionale (ed eventuali sue dipendenze) tra SciPy e SymPy.

Si ringrazia la prof.ssa Fausta D'Acunzo di Preparazione 2.0 per il supporto teorico fornito sul calcolo integrale a più variabili.

Per ottenere il codice si veda il paragrafo Download del codice completo in fondo a questo post.

Integrazione tramite SciPy

Integrale di funzione di una variabile (con estremi finiti)

In analisi matematica, l'integrale definito è un operatore che, dati una funzione reale di una variabile reale e un intervallo $ [a,b] $ (sottoinsieme del dominio), associa alla funzione l'area sottesa dal suo grafico nell'intervallo $[a,b]$.
La libreria SciPy fornisce una serie di metodi numerici per calcolare l'integrale di tali funzioni; scopo di questo capitolo è presentare una serie di demo per mostrare "per esempi" l'uso di tali metodi; per la documentazione esaustiva di essi si invita il lettore a consultare la documentazione ufficiale di SciPy.

Sia dato il seguente integrale di una funzione di una variabile: $$ \int_{1}^{5} 2 x e^{-x} \,dx $$ la cui soluzione analitica è $ \approx 1.3907 $ verificabile online tramite Wolfram Alpha.

Calcolo dell'integrale con quad

Qui di seguito l'esempio di codice Python che calcola l'integrale usando la funzione quad della libreria SciPy:

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

Calcolo dell'integrale con fixed_quad

Qui di seguito l'esempio di codice Python che calcola l'integrale usando la funzione fixed_quad della libreria SciPy:

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

Calcolo dell'integrale con quadrature

Qui di seguito l'esempio di codice Python che calcola l'integrale usando la funzione quadrature della libreria SciPy:

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

Calcolo dell'integrale con romberg

Qui di seguito l'esempio di codice Python che calcola l'integrale usando la funzione romberg della libreria SciPy:

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

Calcolo dell'integrale con trapezoid

La funzione trapezoid è un metodo di integrazione di funzione a campioni fissi, quindi il codice prima discretizza l'intervallo di integrazione in modo uniformemente spaziato e per tutti i valori discreti di $x$ calcola i corrispondenti valore di $y$ e successivamente passa i due insiemi di valori discreti $x$ e $y$ al metodo di integrazione.

Qui di seguito l'esempio di codice Python che calcola l'integrale usando la funzione trapezoid della libreria SciPy:

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

Calcolo dell'integrale con cumulative_trapezoid

Anche la funzione cumulative_trapezoid è un metodo di integrazione di funzione a campioni fissi, e quindi vale quanto detto per trapezoid.
Qui di seguito l'esempio di codice Python che calcola l'integrale usando la funzione cumulative_trapezoid della libreria SciPy:

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

Calcolo dell'integrale con simpson

Anche la funzione simpson è un metodo di integrazione di funzione a campioni fissi, e quindi vale quanto detto per trapezoid.
Qui di seguito l'esempio di codice Python che calcola l'integrale usando la funzione simpson della libreria SciPy:

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

Integrale di funzione di una variabile con quad (con estremo infinito)

Sia dato il seguente integrale di una funzione di una variabile: $$ \int_{1}^{+\infty} 2 x e^{-x} \,dx $$ la cui soluzione analitica è $ \approx 1.4715 $ verificabile online tramite Wolfram Alpha.

Qui di seguito l'esempio di codice Python che calcola l'integrale usando la funzione quad della libreria SciPy:

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

Calcolo della lunghezza di un arco di curva planare

Come è noto, gli integrali di una funzione di una variabile possono essere adoperati anche per calcolare la lunghezza di un arco di curva planare. Se la curva è espressa nella forma $y=f(x)$ ed è continua e derivabile la formula del calcolo della lunghezza dell'arco di curva tra $ x=a $ e $ x=b $ è la seguente: $$ \int_{a}^{b} \sqrt{1 + \left(\frac{\,dy}{\,dx}\right)^2} \,dx $$
Se invece curva è espressa nella forma parametrica $x=f_x(t)$ e $y=f_y(t)$ con entrambe $f_x$ e $f_y$ continue e derivabili, la formula del calcolo della lunghezza dell'arco di curva tra $ t=a $ e $ t=b $ è la seguente: $$ \int_{a}^{b} \sqrt{\left(\frac{\,dx}{\,dt}\right)^2 + \left(\frac{\,dy}{\,dt}\right)^2} \,dt $$
Nota: Nei due esempi che seguiranno, le derivate prime sono calcolate tramite la libreria AutoGrad.

Calcolo della lunghezza di un arco di curva planare espressa in forma esplicita con quad

Sia data la seguente curva planare in forma esplicita: $$ y=e^{-x^2} $$ applicando la corrispondente formula di cui sopra per calcolare la lunghezza dell'arco tra $ x=-1 $ e $ x = 1 $, si ottiene l'integrale: $$ \int_{-1}^{1} \sqrt{1 + \left(\frac{\,d(e^{-x^2})}{\,dx}\right)^2} \,dx $$ la cui soluzione analitica è $ \approx 2.4089 $ verificabile online tramite Wolfram Alpha.

Qui di seguito l'esempio di codice Python che calcola lunghezza di un arco di curva planare espressa in forma esplicita usando la funzione quad della libreria SciPy:

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

Calcolo della lunghezza di un arco di curva planare espressa in forma parametrica con quad

Sia data la seguente curva planare in forma parametrica: $$ x(t)=cos^3 t $$ $$ y(t)=sin^3 t $$ applicando la corrispondente formula di cui sopra per calcolare la lunghezza dell'arco tra $ t=0 $ e $ t = 2\pi $, si ottiene l'integrale: $$ \int_{0}^{2 \pi} \sqrt{\left(\frac{\,d (cos^3 t)}{\,dt}\right)^2 + \left(\frac{\,d (sin^3 t)}{\,dt}\right)^2} \,dt $$ la cui soluzione analitica è $ 6 $ verificabile online tramite Wolfram Alpha.

Qui di seguito l'esempio di codice Python che calcola lunghezza di un arco di curva planare espressa in forma parametrica usando la funzione quad della libreria SciPy:

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

Integrale doppio di funzione di due variabili

In analisi matematica, l'integrale definito doppio è un operatore che, dati una funzione reale di due variabili reali e un insieme compreso nel dominio, associa alla funzione il volume del solido (detto cilindroide) compreso tra la superficie descritta dalla funzione e il piano contenente l'insieme dato.

Sia dato il seguente integrale doppio di una funzione di due variabili: $$ \int_{1}^{5} \int_{y-1}^{y+1} 2 x y e^{-x y} \,dx dy $$ la cui soluzione analitica è $ \approx 1.0273 $ verificabile online tramite Wolfram Alpha.

Calcolo dell'integrale con dblquad

Qui di seguito l'esempio di codice Python che calcola l'integrale usando la funzione dblquad della libreria SciPy:

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

Calcolo dell'integrale con nquad

Qui di seguito l'esempio di codice Python che calcola l'integrale usando la funzione nquad della libreria SciPy:

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

Integrale doppio di funzione di due variabili con nquad (altro esempio)

Sia dato il seguente integrale di una funzione due variabili: $$ \int_{1}^{+\infty} \int_{1}^{+\infty} 2 x y e^{-x y} \,dx dy $$ la cui soluzione analitica è $ \approx 1.17453 $ verificabile online tramite Wolfram Alpha.

Qui di seguito l'esempio di codice Python che calcola l'integrale usando la funzione nquad della libreria SciPy:

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

Integrale triplo di funzione di tre variabili

Sia dato il seguente integrale triplo di una funzione di due variabili: $$ \int_{1}^{2} \int_{z+1}^{z+2} \int_{y+z}^{2(y+z)} x + yz^2 \,dx dy dz $$ la cui soluzione analitica è $ \approx 65.7194 $ verificabile online tramite Wolfram Alpha.

Calcolo dell'integrale con tplquad

Qui di seguito l'esempio di codice Python che calcola l'integrale usando la funzione tplquad della libreria SciPy:

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

Calcolo dell'integrale con nquad

Qui di seguito l'esempio di codice Python che calcola l'integrale usando la funzione nquad della libreria SciPy:

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

Integrazione tramite SymPy

Integrale di funzione di una variabile (con estremi finiti)

In analisi matematica, l'integrale definito è un operatore che, dati una funzione reale di una variabile reale e un intervallo $ [a,b] $ (sottoinsieme del dominio), associa alla funzione l'area sottesa dal suo grafico nell'intervallo $[a,b]$.
La libreria SymPy fornisce una serie di metodi simbolici/analitici per calcolare l'integrale di tali funzioni; scopo di questo capitolo è presentare una serie di demo per mostrare "per esempi" l'uso di tali metodi; per la documentazione esaustiva di essi si invita il lettore a consultare la documentazione ufficiale di SymPy.

Sia dato il seguente integrale di una funzione di una variabile: $$ \int_{1}^{5} 2 x e^{-x} \,dx $$ la cui soluzione analitica è $ \approx 1.3907 $ verificabile online tramite Wolfram Alpha.

Integrale di funzione di una variabile con integrate(f, (x, a, b))

Qui di seguito l'esempio di codice Python che calcola l'integrale usando la funzione integrate(f, (x, a, b)) della libreria SymPy:

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

Integrale di funzione di una variabile con integrate(f, (x, a, b)) tramite integrale indefinito

Qui di seguito l'esempio di codice Python che calcola l'integrale usando la funzione integrate(f, x) della libreria SymPy:

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

Integrale doppio di funzione di due variabili

In analisi matematica, l'integrale definito doppio è un operatore che, dati una funzione reale di due variabili reali e un insieme compreso nel dominio, associa alla funzione il volume del solido (detto cilindroide) compreso tra la superficie descritta dalla funzione e il piano contenente l'insieme dato.

Sia dato il seguente integrale doppio di una funzione di due variabili: $$ \int_{1}^{4} \int_{y-1}^{y+2} x y e^{-x} e^{-y} \,dx dy $$ la cui soluzione analitica è $ \approx 0.396134 $ verificabile online tramite Wolfram Alpha.

Integrale doppio di funzione di due variabili con integrate(f, (x, xa, xb), (y, ya, yb))

Qui di seguito l'esempio di codice Python che calcola l'integrale usando la funzione integrate(f, (x, xa, xb), (y, ya, yb)) della libreria SymPy:

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 del codice completo

Il codice completo è disponibile su GitHub.
Questo materiale è distribuito su licenza MIT; sentiti libero di usare, condividere, "forkare" e adattare tale materiale come credi.
Sentiti anche libero di pubblicare pull-request e bug-report su questo repository di GitHub oppure di contattarmi sui miei canali social disponibili nell'angolo in alto a destra di questa pagina.