Dispense sui vari metodi (rettangoli, trapezi, Cavalieri

Integrazione numerica
Ovvero come calcolare numericamente l'area sottesa al grafico di una funzione.
Supponiamo di avere una funzione
f (x ) integrabile nell'intervallo (a,b). L'integrale può essere stimato
a) suddividendo la regione sottesa alla funzione in tanti rettangolini (metodo dei rettangoli)
In questo modo, suddividendo l'intervallo (a,b) in N parti, si ha che la stima dell'area risulta essere pari alla somma delle
aree di tutti i rettangolini, ovvero:
b
∫a
f ( x) dx ≃
(b−a) N
∑i=0 f ( x i)
N
b−a
la base dei rettangolini, l'errore commesso nell'approssimazione andrà linearmente con dx ,
N
ovvero ε ∝ dx . L'integrale di una funzione costante ( f ( x )=q )sarà calcolato in modo esatto.
Detta dx=
b) suddividendo la regione sottesa alla funzione in tanti trapezi (metodo dei trapezi), migliorando in questo modo il
precedente metodo
In questo modo, suddividendo l'intervallo (a,b) in N parti, si ha che la stima dell'area risulta essere pari alla somma delle
aree di tutti i trapezi, ovvero:
b
∫a
f ( x) dx ≃
(b−a)
⋅
N
(
f (a)+ f (b)
+
2
N−1
∑i=1
f ( x i)
)
b−a
la base dei rettangolini (ora trapezi), l'errore commesso nell'approssimazione andrà quadraticamente
N
con dx , ovvero ε ∝ dx 2 . L'integrale di una funzione lineare ( f (x )=m x +q ) sarà calcolato in modo esatto.
Detta dx=
c) suddividendo la regione sottesa alla funzione in tanti archi di parabola (metodo dei Cavalieri-Simpson) migliorando
così entrambi i precedenti metodi.
Si suddivide l'intervallo (a,b) in N parti. Si considerano 3 punti alla volta e si determina l'equazione della parabola
passante per quei tre punti. A questo punto si determina l'area sottesa a quell'arco di parabola, pari (dopo un po' di
dx
b−a
f (x 0 )+4 f ( x 1)+ f ( x 2 ) ) con dx=
passaggi) a
. Sommando su tutti gli N intervalli si ottiene
(
3
N
b
∫a
f ( x) dx ≃
(b−a)
⋅ ( f (x 0 )+ 4 f ( x 1 )+ 2 f (x 2 )+4 f ( x 3) + .. . + 4 f ( x N −3 )+ 2 f (x N −2)+ 4 f ( x N −1 )+ f ( x N ) )
N
L'errore commesso nell'approssimazione andrà come ε ∝ dx 3 . L'integrale di una funzione quadratica (
2
f ( x )=a x + b x +c ) sarà calcolato in modo esatto.
d) utilizzando un metodo stocastico (o Monte Carlo).
In questo caso si calcola l'area come nel caso a) prendendo però i punti x i non più a distanze regolari uno dall'altro
ma prendendo N valori x i a caso nell'intervallo (a,b) . È evidente che se il campionamento di x i è uniforme e che
se N è sufficientemente grande otteniamo un valore ragionevole per la stima dell'integrale:
b
(b−a) N
∫a f ( x) dx ≃ N ∑i=0 f ( x i)
L'errore va come ε ∝ √ (dx ) , quindi “peggio” che nel caso a). L'utilità e l'efficienza dell'algoritmo risulta evidente
nell'integrare una funzione a più variabili (integrando su 2 o più dimensioni).
d) Metodo stocastico (Monte Carlo)
Stima di pi greco con l'uso di metodi stocastici.
Parte prima
Con la funzione rand() otteniamo un numero intero casuale compreso fra 0 e la costante predefinita RAND_MAX. Per ottenere
un numero x reale compreso fra 0 e 2 dovremo quindi scrivere:
x=2.*(float)rand()/(float)RAND_MAX;
(1)
Le due istruzioni (float) nella precedente servono per convertire quanto sta a destra ad una variabile di tipo float e vengono
dette “cast”, ovvero conversione di tipo. Senza tali istruzioni la divisione dell'istruzione precedente risulterebbe essere una
divisione fra interi (quali sono appunto l'output della funzione “rand” e la costante predefinita RAND_MAX), e quindi il risultato
verrebbe troncato all'intero, anche se la variabile x è stata definita come float.
Proviamo quindi a calcolare l'area di un cerchio sparando dei punti a caso, compresi nel quadrato
circoscritto al cerchio. È chiaro che la probabilità di trovare un punto all'interno del cerchio è pari al
rapporto fra l'area del cerchio e quella del quadrato.
2
Campioneremo quindi un numero di punti N TOT sufficientemente grande, di coordinate (x,y), con x e y
0
calcolati come dalla precedente formula.
2
Per ogni punto verificheremo se è interno alla circonferenza, ovvero se la distanza del punto (x,y) dal centro del cerchio è
minore uguale al raggio.
Il rapporto del numero dei punti interni NINT alla circonferenza sul numero dei punti totali NTOT soddisferà la seguente:
lim
N TOT  ∞
N lNT  r 2 
=
=
N TOT 4 r 2 4
(2)
Limitandoci ad un NTOT finito otterremo una stima approssimata di π
≃4
N lNT
N TOT
(3)
Parte seconda
Con una piccola modifica vediamo come migliorare l'algoritmo. Considera la figura seguente y
corrispondente ad un quarto di cerchio.
L'arco di circonferenza rappresentato può essere visto come una funzione di equazione
y= f  x= r 2− x 2
per 0 < x < r
(4)
Il rapporto fra le due aree sarà sempre pari a π/4. Come sopra potremmo campionare un insieme di
punti (x,y) distribuiti uniformemente nel quadrato e proseguire allo stesso modo.
a
b
Possiamo notare però che in questo caso non è più necessario campionare la y, e si può migliorare
l'efficienza dell'algoritmo. Per una data x la probabilità che il punto sia interno alla circonferenza (ovvero sotto la curva di
equazione (4)) sarà proporzionale al valore della funzione in quel punto. Ovvero, indicando con {xi} l'insieme delle x
campionate, con 1 < i < NTOT , avremo che
N TOT
f  xi b−a 

= lim ∑
4 N  ∞ i =1
N TOT
(5)
TOT
Parte terza
Siamo partiti con l'intenzione di stimare semplicemente il valore di π. Se osserviamo però la formula (5) ci accorgiamo come
questa si possa usare per il calcolo dell'area sottesa ad una generica funzione nell'intervallo (a,b).
Volendo, per esempio, calcolare l'area sottesa al grafico della funzione
g  x=exp x −1
(6)
nell'intervallo (a=2,b=5) posiamo stimarla semplicemente campionando un numero NTOT di variabili stocastiche xi
uniformemente distribuite nell'intervallo (a,b), ovvero:
x
xi=(float)rand()/(float)RAND_MAX * (b-a) + a;
(7)
e, allo stesso modo dell'equazione (5), stimare l'area con la seguente:
N TOT
A≃ ∑
i=1
g  x i b−a
N TOT
(8)
che in questo caso risulta uguale 138.024.
Curiosità
Come probabilmente avrai notato, il metodo dei “rettangoli” risulta essere più accurato a parità di costo computazionale
(tempo impiegato) per calcolare l'area. L'approccio appena visto sembrerebbe per tanto una inutile complicazione. Così però
non è se si considera un problema in più dimensioni, ovvero se si considera una funzione di più variabili, f(x,y,z....).
Facciamo un esempio.
Per il calcolo con una certa accuratezza dell'area sottesa ad una funzione f(x) si
necessita di un campionamento di almeno NTOT punti, ovvero del calcolo
dell'area di NTOT rettangolini.
Per analogo calcolo con analoga accuratezza del volume(e non l'area, ora c'è
una dimensione in più!) sotteso alla curva di una funzione f(x,y) si necessiterà di
un campionamento di almeno (NTOT * NTOT)= NTOT2 punti, ovvero del calcolo del
volume di NTOT 2 parallelepipedi.
Ti accorgi così che il numero di parallelepipedi da considerare va
esponenzialmente con il numero di variabili della funzione f. Ovvero calcolare l'iper-volume sotteso ad una funzione
f(x1,x2,...xK) costerà il calcolo di NK iper-volumetti. Al contrario col metodo stocastico appena visto questo non è altrettanto
vero, e per problemi a più dimensioni risulta pertanto essere molto più efficiente.
a) Regola del rettangolo
Da Wikipedia, l'enciclopedia libera.
La regola del rettangolo o regola del punto medio, è il più semplice procedimento di integrazione numerica per approssimare
un integrale definito nella forma :
.
Tale formula approssima l'integrale (e quindi l'area sottesa dalla funzione) come un rettangolo di base
e di altezza
dove a e b sono gli estremi di integrazione e c è il punto medio dell'intervallo, ottenendo un'espressione finale per l'integrale
,
pari a:
Formula composita
Illustrazione del metodo del punto medio composito
Per calcolare con più accuratezza l'integrale, si divide l'intervallo di integrazione
in
sottointervalli di ampiezza
uniforme pari a
La formula del punto medio diventerà dunque
sottointervallo.
, dove
rappresenta il punto medio del k-esimo
Analisi dell'errore
L'errore sviluppato con il metodo del rettangolo, assumerà la seguente espressione:
, dove
è un opportuno punto compreso nell'intervallo
.
Nel caso si usi il metodo composito l'errore sarà
. Dalla formula dell'errore si deduce che il metodo
integra esattamente polinomi di primo grado, e che l'errore diminuisce quadraticamente rispetto all'ampiezza dei sottointervalli
.
b) Regola del trapezio
Da Wikipedia, l'enciclopedia libera.
Regola del trapezio.
In analisi numerica, la regola del trapezio fornisce un procedimento per il calcolo approssimato di un integrale definito della
forma
.
La regola del trapezio o di Stevin, nella sua formulazione elementare, propone di approssimare l'integrale, cioè l'area della
regione piana compresa fra il grafo della funzione
(b,0) e (a,0). Di conseguenza
e l'asse delle ascisse con l'area del trapezio di vertici (a,f(a)), (b,f(b)),
.
Come è visivamente intuibile, questa approssimazione è accettabile se nell'intervallo di integrazione la funzione ha un
andamento che si scosta poco dal lineare. Se questo non accade si può suddividere l'intervallo complessivo
in un numero
n opportuno di sottointervalli in ciascuno dei quali in genere accade che la funzione ha andamento poco lontano dal lineare;
quindi la regola del trapezio nella forma composta dice di applicare l'approssimazione precedente a tutti i sottointervalli. Si
ottiene quindi la formula
La regola del trapezio fa parte della famiglia di formule per l'integrazione numerica chiamate formule di Newton-Cotes. Di
questa famiglia fa parte anche la regola di Cavalieri-Simpson che dà spesso risultati più accurati. La regola di CavalieriSimpson e altri metodi simili migliorano la regola dei trapezi per gran parte delle funzioni dotate di due derivate continue.
Accade però per certe funzioni dal comportamento irregolare che la più semplice regola del trapezio risulti preferibile. Inoltre la
regola del trapezio tende a diventare molto accurata per gli integrali di funzioni periodiche nei rispettivi intervalli di periodicità;
questo comportamento viene chiarito in relazione alla formula di Eulero-Maclaurin.
c) Regola di Cavalieri-Simpson
Da Wikipedia, l'enciclopedia libera.
Per regola di Cavalieri-Simpson o regola di Simpson si intende un metodo per il calcolo numerico approssimato di integrali
definiti della forma:
.
Come tutti i procedimenti per il calcolo approssimato di integrali definiti e per altri calcoli approssimati a partire da funzioni di
variabile reale, tale metodo si utilizza per funzioni
delle quali non si conosce la funzione primitiva, oppure della cui
primitiva si conoscono solo caratteristiche dalle quali non si riesce a ricavare una espressione tramite funzioni elementari che
possa essere ragionevolmente utilizzata per i calcoli richiesti. Questi metodi approssimati si utilizzano inoltre nei casi in cui non
è nota una espressione analitica della funzione da integrare, ma si conoscono soltanto alcuni suoi valori (ottenuti
sperimentalmente o ricavati da altre fonti), oppure quando è noto soltanto il suo diagramma (tracciato con l'ausilio di appositi
strumenti o ricavato dalla letteratura).
La formula di quadratura o metodo delle parabole
La regola di Cavalieri-Simpson approssima l'integrale della funzione richiesta (in blu) con quello della parabola che la interpola
nei nodi (in rosso)
La regola di Cavalieri-Simpson prevede la suddivisione dell'intervallo di integrazione in sottointervalli e la sostituzione in
questi sottointervalli della funzione integranda mediante archi di parabola, cioè mediante polinomi quadratici.
Consideriamo dunque
integrazione
.
Suddiviso
; per semplicità di raffigurazione supponiamo sia
in un numero pari n=2m di sottointervalli, ciascuno di ampiezza
in tutto l'intervallo di
. Introduciamo poi le notazioni
per gli estremi dei successivi sottointervalli e per i valori che la funzione assume in loro corrispondenza.
Consideriamo anche l'intervallo parziale formato da due sottointervalli consecutivi avente come estremi
consideriamo anche i successivi m-1 intervalli parziali aventi come estremi rispettivamente
e , ... ,
e
; oltre a questo
e .
In ciascuno di questi intervalli parziali ci proponiamo di sostituire
con una funzione razionale intera di secondo grado.
Cominciamo dal primo intervallo parziale e scegliamo un polinomio della forma
in modo che il suo integrale tra
trascurabile.
e
differisca da quello della funzione originale di una quantità che possa risultare
La espressione polinomiale sostitutiva rappresenta una generica parabola con asse di simmetria verticale. Per determinare il
valore delle costanti A, B e C si impone il passaggio della parabola per i punti di coordinate:
.
In questo modo la parabola è univocamente determinata dalla risoluzione del seguente sistema di equazioni lineari:
,
da cui risulta:
Per il valore dell'integrale di questo polinomio che scriviamo
si trova:
=
.
Sostituendo i valori di A, B e C ricavati dal sistema, si ottiene il valore approssimato
.
Operiamo in modo analogo per il calcolo degli integrali dei polinomi nei successivi m-1 intervalli parziali; successivamente si
sommano i valori ottenuti sugli m intervalli parziali e per l'intero intervallo d'integrazione si ottiene un valore approssimato che
denotiamo J per l'integrale da valutare:
.
Dunque:
.
Tale formula prende il nome di formula di Cavalieri-Simpson o formula di quadratura delle parabole.
Gli errori
Il metodo di quadratura di Cavalieri-Simpson, come ogni metodo di approssimazione numerica, è suscettibile di errori. Oltre
all'errore dovuto alla sostituzione della funzione integranda con una sequenza di funzioni approssimanti, i polinomi quadratici,
intrinseco al metodo utilizzato, si riscontrano anche errori dovuti all'arrotondamento dei valori che vengono concretamente
calcolati con strumenti che inevitabilmente operano con precisione limitata.
Per ridurre al minimo questi ultimi, è consigliabile:
• scegliere un passo di integrazione
con un numero finito di cifre decimali;
• eseguire i calcoli con un numero di cifre decimali almeno doppio di quello delle cifre che si desiderano esatte nel
risultato.
Indicando l'errore intrinseco al metodo
, si può dimostrare che:
,
dove k è una costante che dipende dalla funzione integranda e dall'intervallo d'integrazione. La regola di Cavalieri-Simpson è
dunque un metodo del quarto ordine.
Di questo errore può essere molto utile conoscere una maggiorazione; la valutazione accurata di tale maggiorazione non è
semplice, poiché richiede di calcolare la derivata quarta della funzione integranda. Per molte funzioni integrande date
analiticamente il calcolo della derivata quarta risulta molto oneroso; per funzioni note empiricamente la stessa valutazione della
derivata quarta costituisce di per se un problema di calcolo approssimato tendenzialmente oneroso. Di conseguenza in genere
per la valutazione dell'errore si preferisce ricorrere a metodi empirici: il più noto e utilizzato è il metodo del dimezzamento del
passo.
Da quanto osservato precedentemente segue che, applicando il metodo di Cavalieri-Simpson con un passo di integrazione h, si
ottiene l'approssimazione che ora indichiamo con
Utilizzando il passo di integrazione
affetta da un errore che scriviamo
, si otterrà il valore approssimato
con errore:
.
.
Da tali relazioni segue:
,
da cui risulta:
.
Poiché, trascurando le approssimazioni da arrotondamento, l'approssimazione migliore è data da
in
, sostituendo il valore di k
, si ottiene:
.
Si può quindi assumere come valore assoluto di
:
come maggiorazione dell'errore assoluto
È interessante osservare che, se le approssimazioni
.
e
coincidono per le prime r cifre decimali, risulta:
il che equivale a dire che le prime r cifre decimali non sono affette da errore.
Si può quindi concludere che se due approssimazioni di un integrale, di cui la seconda ottenuta dimezzando il passo di
integrazione utilizzato per calcolare la prima, coincidono per le prime r cifre decimali, tali cifre si possono ritenere esatte.
Più in generale se si vuole conoscere un'approssimazione di un integrale con la garanzia dell'esattezza per un determinato
numero s di cifre decimali, si deve calcolare un certo numero di approssimazioni successive, dimezzando di volta in volta il
passo, fino ad ottenerne due che coincidono per s cifre.
Osserviamo che può accadere di arrivare alla prima coppia di valori approssimati soddisfacenti che presentano più di s cifre
coincidenti. Osserviamo anche che procedendo con la riduzione dell'ampiezza dei sotto-intervalli, oltre al maggior tempo di
calcolo richiesto, si possono avere errori di arrotondamento tutt'altro che trascurabili a causa dell'aumento del numero di
operazioni richieste; si può anche arrivare a situazioni che vedono aumentare l'errore complessivo con il ridursi del passo di
integrazione. Per tale motivo sono stati concepiti dei metodi di integrazione numerica "adattivi" (o adattativi) che aumentano il
numero dei sotto-intervalli solo nelle zone indicate da un apposito test di errore.
Nota storica
Thomas Simpson è conosciuto al giorno d'oggi soprattutto per la regola di Simpson, il termine prevalente a livello
internazionale. In realtà Simpson utilizzò nelle sue opere una formulazione comunque già diffusa, ma ancora non formalizzata
all'epoca in cui scriveva. Il procedimento era stato trovato 200 anni prima da Keplero e molti testi tedeschi la chiamano
Keplersche Fassregel. Inoltre la formula era stata introdotta e dimostrata dal matematico italiano Bonaventura Cavalieri nel
1635 nella sua formulazione geometrica; per questo motivo in molti testi italiani si vuole rendere omaggio anche a Cavalieri.
Occorre aggiungere che questa regola era nota anche a Evangelista Torricelli. Per quanto riguarda il mondo britannico, sembra
che la regola che porta il nome di Simpson fosse conosciuta e utilizzata in precedenza; essa è stata utilizzata anche da James
Gregory. Inoltre si pensa che nella nota serie televisiva "I Simpsons" il cognome sia stato scelto in base alla storia di Thomas,
in realtà non meritevole della scoperta a lui associata.