didattica, ricerca

I miti dell’integrazione numerica

Molti programmatori credono che l’uso di algoritmi d’integrazione di ordine superiore, combinato con un gran numero di divisioni dell’intervallo d’integrazione, sia utile (e talvolta necessario) per ottenere una buona precisione. In questo articolo dimostriamo che non è sempre vero.


Questo post è la traduzione in italiano di un articolo che ho scritto per medium.com, che ha raggiunto un gran numero di lettori (quasi 2000) in soli due giorni e ricevuto molti apprezzamenti, che interpreto come un segno che indica che l’articolo è stato utile e che, quindi, i miti di cui parlo sono molto diffusi. L’articolo in questione è stato ispirato da una storia scritta da Kazi Abu Rousan, che illustrava i vantaggi dell’utilizzo della formula dei trapezi per l’integrazione numerica, debitamente e magistralmente illustrata.

In effetti, l’articolo è corretto in ogni dettaglio, ma la conclusione cui può portare potrebbe non essere tale. Permettetemi d’iniziare con l’esempio fornito nell’articolo di Kazi in cui si calcola l’integrale di 1/(1+x²) tra 0 e 1, alla ricerca del numero n di passi in cui si deve dividere l’intervallo per ottenere una precisione dell’ordine di 10⁻⁶. Per l’occasione ho leggermente modificato il programma, allo scopo di iterare la procedura in modo da individuare il valore di n per il quale il resto (cioè la differenza tra il valore effettivo dell’integrale e la sua stima numerica) è al di sotto della tolleranza data. Utilizzando la formula del trapezio, ottengo la precisione data dividendo l’intervallo [0,1] in n=204 parti. Utilizzando la formula dei rettangoli, il numero di divisioni necessarie per ottenere la stessa precisione è n=145. Un momento…come può la divisione dell’intervallo in meno sotto-intervalli, e, per di più, utilizzando una regola d’integrazione di ordine inferiore, portare a un risultato coerente con un metodo di ordine superiore con più iterazioni?

Per comodità, riporto l’algoritmo d’integrazione, così come da me modificato.

original_value = np.pi/4def rectangles(f, a, b, tol):
    sum = 0
    n = 1
    dx = (b - a)/n
    while abs(sum - original_value) > tol:
        sum = 0
        dx = (b - a)/n
        x = a + 0.5*dx
        for i in range(0, n):
            sum += f(x)
            x += dx
        sum *= dx
        n += 1
    return sum, n-1

In effetti, un trucco c’è, e consiste nel fatto che ho usato una versione speciale della formula dei rettangoli, nota come regola del punto di mezzo. In questo caso, l’altezza del rettangolo, usata per approssimare la funzione, non è f(a), né f(b), ma f((a+b)/2). È facile dimostrare che l’area di un rettangolo la cui altezza è la media delle lunghezze delle basi di un trapezio è equivalente all’area di quest’ultimo. I due metodi sono, dunque, del tutto equivalenti. La differenza in n a vantaggio della formula del punto di mezzo è dovuta principalmente alle proprietà specifiche di f(x) (in particolare, al fatto che f(x) è monotòna decrescente).

Il metodo del punto di mezzo consiste nell’interpolazione di una curva in un intervallo [a,a+h], utilizzando una costante corrispondente a f(ξ), dove ξ è il punto medio dell’intervallo. È equivalente al metodo dei trapezi, ma il suo errore è inferiore se la funzione è monotòna nell’intervallo [figura fatta dall’autore e pubblicata su “Programmazione Scientifica”, ed. Pearson]

In effetti, la formula del punto di mezzo è un caso particolare dei cosiddetti metodi d’integrazione di Gauss, che sfruttano le simmetrie e altre proprietà della funzione da integrare, per fare una scelta furba dei punti in cui f(x)=P(x,m), dove P(x,m) è un polinomio di grado m.

Si può dimostrare che il resto di un algoritmo d’integrazione numerica, definito come la differenza tra i valori approssimati e quelli reali di un integrale, si può scrivere come

dove C è una costante che dipende dalla derivata m-esima di f(x), m è il grado del polinomio usato per interpolare f(x) e n il numero d’intervalli in cui dividere [a, b]. Per la formula dei rettangoli, m=1, mentre m=2 nel caso della formula dei trapezi e del punto di mezzo (che, come mostrato sopra, è equivalente a quella dei trapezi). Non è qui il caso di dare una dimostrazione formale di questo risultato, ma è facile capire che più grande è l’intervallo [a, b], più grande è l’errore, quindi 𝛿 aumenta con b-a e diminuisce con n; inoltre, la velocità con cui l’errore aumenta o diminuisce dipende dal grado m del polinomio utilizzato per approssimare f(x).

Al fine di ottenere un resto appena sotto 10⁻⁶, abbiamo avuto bisogno di n=145 nel caso più favorevole. Nella figura sotto, mostro come evolve l’integrale numerico I(n) con l’aumentare di n, se stimato col metodo del punto di mezzo per il quale la “m efficace” è m=2:

Il resto tende a zero come 1/n² per n che tende all’infinito. Quindi, se tracciamo I(n) in funzione di 1/n², secondo la formula del resto, ci aspettiamo che I(n) segua un andamento lineare. In effetti, ecco come I(n) cambia con 1/n²:

Quindi, per ottenere la stessa precisione di prima, è sufficiente calcolare numericamente lo stesso integrale per n=3, 5, 7 e 9, e poi interpolare i dati con una retta, la cui intercetta coinciderà con la stima dell’integrale per n=∞. Il resto che si ottiene, in questo caso, è 1.5×10⁻⁷, addirittura un ordine di grandezza migliore di quello ottenuto con n=145!

Eseguire l’interpolazione è facilissimo, in Python. Bastano due righe di codice:

from scipy import optimize
p, cov = optimize.curve_fit(line, xx, Integration_values)

Nel codice xx è una lista contenente i valori 1/n² e Integration_values è la corrispondente lista delle stime dell’integrale. line è una funzione definita come

def line(x, a, b):
return a * x + b

L’integrale sarà così dato da p[1]. Il suo resto si stima come la radice quadrata dell’elemento corrispondente della matrice di covarianza cov[1][1]:

I = params[1]
rem = np.sqrt(cov[1][1])
print("I: {:.2e}+-{:.2e}".format(I, rem))

Nel nostro caso otteniamo

I: 7.85e-01+-6.93e-08

Alla fine, valutando f(x) soltanto N=3+5+7+9=24 volte, abbiamo ottenuto una migliore precisione rispetto al caso n=145, che richiede di valutare f(x) 145 volte. Non male, no?

Naturalmente, l’apparente guadagno di un fattore 145/24≃6 nel tempo di calcolo è, almeno in parte, vanificato dalla necessità di interpolare i dati. Tuttavia, un fit lineare è, in realtà, estremamente semplice, e si può realizzare con un algoritmo molto rapido. In effetti, sono stati necessari 0,0066 s per calcolare l’integrale usando la formula del punto di mezzo con n=145, da confrontare con i 0,00012 s necessari con il metodo dell’interpolazione.

Rispondi

Inserisci i tuoi dati qui sotto o clicca su un'icona per effettuare l'accesso:

Logo di WordPress.com

Stai commentando usando il tuo account WordPress.com. Chiudi sessione /  Modifica )

Google photo

Stai commentando usando il tuo account Google. Chiudi sessione /  Modifica )

Foto Twitter

Stai commentando usando il tuo account Twitter. Chiudi sessione /  Modifica )

Foto di Facebook

Stai commentando usando il tuo account Facebook. Chiudi sessione /  Modifica )

Connessione a %s...