Follow
In parecchie applicazioni informatiche, specie in quelle simulative, saper integrare una funzione ricopre un ruolo fondamentale. Purtroppo però le tecniche di integrazione che ci hanno insegnato a scuola o nei corsi di Analisi I sono del tutto inutili. Raramente, infatti, ci troveremo a che fare con funzioni integrabili, anzi, spesso dovremo integrare funzioni non esprimibili tramite funzioni analitiche.
Supponete ad esempio di voler implementare un piccolo gioco di guida. La prima cosa da fare è progettare il nostro sistema di guida. Ovviamente, dato che siamo dei tipi tosti, vogliamo che sia il più realistico possibile e non c’è nulla di più realistico che usare le leggi della fisica.
Supponiamo che il giocatore controlli direttamente l’accelerazione della macchina (una supposizione molto ragionevole) e che tale accelerazione sia descritta da una funzione a(t). A meno che il vostro giocatore non sia un androide o una specie di divinità è plausibile pensare che l’accelerazione imposta alla macchina non sia una precisa funzione analitica bensì un guazzabuglio di accelerate e frenate e, anche se il pilota è talmente bravo da accelerare seguendo una funzione analitica, voi non potreste sapere in anticipo quale funzione sia. Nonostante ciò dovete derivare lo stesso questa funzione sconosciuta per ricavare la velocità (e integrare di nuovo per ottenere la posizione della macchinina): non avete altra scelta che usare dei metodi di integrazione numerica.
I metodi di integrazione numerica permettono di stimare l’integrale di una funzione conoscendo solamente il suo valore in un insieme finito di punti. Ho selezionato tre metodi elencati in ordine di complessità. Esistono anche alcuni metodi più precisi (come, ad esempio, i metodi di integrazione gaussiana) ma necessitano di molti più calcoli e non sono quindi adatti all’uso “informatico” se non in alcuni campi specifici.
Metodo dei Rettangoli
Questo metodo è sicuramente il più semplice e, di contro, il meno preciso. Supponiamo di voler integrare una funzione f(x) nell’intervallo [a,b]. Per fare ciò dividiamo l’intervallo in n intervalli di dimensione delta.
Ogni rettangolo avrà un area di f(x1)*delta. Per stimare il nostro integrale ci basta sommare i contributi di ogni rettangolo. Ovviamente, più il numero di intervalli è alto più il nostro integrale sarà ben approssimato.
Una versione delmetodo dei rettangoli in python è data dal seguente codice:
total = 0
i = a
while i <= b :
total += function2(i)
i += delta
return total*delta
Total è l’accumulatore del risultato. Il puntatore i va da a a b con passo delta scandendo quindi intervallo per intervallo la funzione da integrare (rappresentata da function2).
Metodo dei Trapezi
Il metodo dei trapezi è la logica estensione del metodo dei rettangoli. Invece di calcolare l’area del rettangolo di base delta e altezza f(x_n) calcoliamo l’area del trapezio il cui lato ubiquo unisce f(a+n*delta) con f(a+(n+1)*delta).
Il metodo è molto più preciso e può essere scritto in python con una sola riga
return (b-a) * ( f(a)/2 + f(b)/2 + sum([f(a + (b-a)*k/N) for k in xrange(1,N)]) ) / N
Metodo di Cavalieri-Simpson
Il metodo migliore è però il metodo di Cavalieri-Simpson. Il metodo approssima due intervalli con una parabola passante per il punto in comune fra i due.
Ogni coppia di intervalli ha come integrale
Dove h è la dimensione di un intervallo, y0 è il valore della funzione nel primo estremo del primo intervallo, y2 è il valore della funzione nel secondo estremo e y1 è il valore della funzione nell’estremo in comune.
L’algoritmo in python è il seguente
total = 0
i = 0
while i+2*delta <= 10 :
value1 = function2(i)
value2 = function2(i+delta)
value3 = function2(i+2*delta)
value = (value1 + 4*value2 + value3)
total += value
i += 2*delta
return (total * delta) / 3
Confronti
Ho eseguito un breve test sui tre metodi eseguendo il seguente integrale:
Test con 10 intervalli
- Rettangoli : 2.124518813379953
- Trapezi: 1.6517198689244215
- Simpson: 1.658728115119567
Già con soli 10 intervalli notiamo che il metodo di Simpson ha un errore inferiore al millesimo.
Test con 100 intervalli
- Rettangoli : 1.7055620908116773
- Trapezi: 1.6582821963661236
- Simpson: 1.6583476291805812
Test con 1000 intervalli
- Rettangoli : 1.6630749297713836
- Trapezi: 1.6583469403268327
- Simpson: 1.6583475942223707
Per un test più preciso bisognerebe eseguire il test anche per molti altri tipi di funzione. In ogni caso, se vi trovate a dover integrare una funzione in un vostro programma, ora, sapete come fare.
Non conoscevo il Metodo di Cavalieri-Simpson. Molto interessante.
è meraviglioso nella sua precisione e semplicità implementativa 🙂
Se ti capiterà di calcolare un integrale finito numericamente usa il metodo di Gauss (vedi ad esempio http://it.wikipedia.org/wiki/Quadratura_di_Gauss) e non quei metodi naif.
Inoltre usare il python (o qualunque altro linguaggio interpretato) per il calcolo numerico è una scelta perdente, rispetto a linguaggi compilati i tempi sono enormemente più lunghi
Conosco il metodo di Gauss ma, come ho specificato nell’articolo è inutilmente complesso e non estendibile direttamente ad un metodo iterativo utile, ad esempio, ad integrare in real-time sistemi dinamici (o almeno io non ho trovato nulla sull’argomento).
Al massimo si potrebbe utilizzare il metodo di Runge-Kutta che dovrebbe dare risultati migliori ma sto ancora studiandomelo. 🙂
Inoltre l’uso del Python era ovviamente solo a scopo dimostrativo e didattico.
Non lasciarti trarre in inganno dalle formule per il calcolo dei coefficienti. Una volta deciso il polinomio, i coefficienti sono calcolati una volta per tutte e stanno in una tabella, assieme alle ascisse ottimali. A parità di punti, Gauss dà circa il doppio di cifre decimali corrette rispetto a Simpson.
Il problema nell’uso in real time è che la cadenza con cui ti arrivano i dati non corrisponde alle ascisse necessarie per Gauss. Di solito i valori campionati vengono filtrati; se si usa un filtraggio polinomiale (es. Savitzky-Golay) l’integrazione può essere fatta sul polinomio.
Allora vedrò di approfondire di più quest’aspetto 🙂 Tuttavia in pratica una precisione così elevata solitamente non è mai necessaria.
Questi aspetti sono sempre fonte di discussione fra chi ha una visione matematica e chi ha una visione ingegneristica xD
Che io sappia Runge-Kutta serve per integrare equazioni differenziali; non sono al corrente di un metodo con questo nome per l’integrazione di funzioni date per punti.
Si, integra equazioni differenziali ed è proprio quello che serve nel caso di applicazioni “simulative”. Ovviamente è un approssimazione ma più che sufficiente se “filtrata” (ad esempio con il favoloso filtro di Kalman).
A robotica la usiamo per passare da un modello dinamico espresso in forma di equazione differenziale ad un modello discreto utile per la computazione in realtime. Io l’ho usato con soddisfazione anche per piccole simulazioni fisiche in alcuni giochi. 🙂