Integrazione numerica

Regola dei trapezi - Metodo di Simpson
ed estensioni di ordine superiore
Torna
all'indice
generale
Supponiamo di voler calcolare l'integrale di una funzione f su un certo intervallo xa , xb e di conoscere la funzione
solo agli estremi f xa , f xb . Senza ulteriori ipotesi abbiamo solo due condizioni e possiamo calcolare l'integrale di f
approssimandola con una retta che passa per i punti f x a , x a e f x b , x b . La retta ha equazione
r( x) a 0
a 1.x
dove a 0 , a 1 sono le soluzioni del sistema
r xa
f xa
r xb
f xb
1 xa
ovvero
.
1 xb
a0
f xa
a1
f xb
che ha per soluzione
1 xa
1
.
1 xb
f xa
x b. f x a
xb
simplify
f xb
f xa
xb
x a. f x b
xa
f xb
xa
L'equazione della retta che passa per x a e x b è quindi
x b.f x a
xb
x a. f x b
xa
f xa
xb
.
f xb
xa
1
x b. f x a
xb
collect , x
x
x a. f x b
xa
f xa
xb
f xb
.x
xa
Integrando tra x a e x b
xb
x b. f x a
xb
x a. f x b
xa
f xa
xb
f xb
. x d x factor
xa
1.
f xb
2
f xa . xb
xa
xa
nota come regola dei trapezi. Con una procedura analoga si calcola anche l'ordine successivo. Conoscendo la
a b
funzione anche in un punto intermedio dell'intervallo avremo infatti (f ( a ) f 0 , f
f 1 ,f ( b ) f 2 )
2
x0
2. h
1
1
x0
1 x0
1
2
x0
x0
h
2. h
x0
x0
f0
h
2
2.h
.
2
f1
f2
1
.
x d x factor
2
1. . .
h 4f1
3
f2
f0
x
x0
questa espressione è come regola di Simpson
Con la stessa procedura possiamo calcolare anche approssimazioni di ordine superiore
13
Suddividendo l'intervallo in 3 parti si ottiene la regola 3/8 di Simpson
x0
3. h
1
1
2
x0
x0
h
x0
1
3
x0
x0
2
h
x0
1 x0
2. h
x0
2
2.h
1 x0
3. h
x0
3.h
2
f0
3
h
.
1
f1
x
.
2
x0
3
2.h
f2
x
x0
3.h
3
f3
x
3. . .
h 3f1
8
d x factor
3.f 2
f0
f3
3
x0
Suddividendo l'intervallo in 4 parti si ottiene la regola di Bode
x0
4. h
1
1
x0
3
2
x0
h
x0
2
h
x0
1 x0
2. h
x0
2. h
2
1 x0
3. h
x0
3. h
2
1 x0
4. h
x0
2
4. h
1
4
x0
x0
x0
3
h
x0
x0
2. h
3
x0
3. h
3
x0
3
4. h
4
h
x0
2. h
4
x0
3. h
4
x0
4
4. h
.
x0
=
f0
1
f1
x
f2
.
2
x
d x=
3
f3
x
f4
x
4
2. . .
h 7 f0
45
32 . f 1
12 . f 2
32 . f 3
7 . f4
Analizziamo ora che accuratezza ci dobbiamo aspettare da un generico ordine di sviluppo. Le formule indicate
all'ordine n sarebbero corrette se f fosse esattamente un polinomio di grado n. La funzione f tuttavia si differenzia in
generale da un polinomio (in effetti dal suo polinomio di Taylor di grado n) in quanto ha le derivate di ordine n 1
diverse da zero. L'errore dovuto ad un troncamento del polinomio all'ordine n 1 su un intervallo h è di ordine
n (n 1)
n 1
O h .f
. Integrando la funzione approssimata di conseguenza il risultato sarà di ordine O h
se n è dispari.
n
2
Sarà invece di ordine O h
se n è pari perchè il termine dispari relativo all'errore grazie al fatto di aver scelto le
posizioni su cui si calcola la funzione equispaziate, e quindi simmetriche sull'intervallo di integrazione, si annulla.
Riassumendo avremo quindi che
Trapezi
1
Simpson
2
h.
f xb
2
f1
Simpson 3/8 3
Bode
4
3
4. f 2
3. .
h f3
8
3.f 1
2. . .
h 7 f0
45
32 . f 1
3 (2)
O h .f
f xa
f3
3. f 2
12 . f 2
5 (4)
O h .f
.h
3
5 (4)
O h .f
f0
32 . f 3
7 . f4
7 (6)
O h .f
In generale tuttavia non siamo interessati al calcolo di un integrale su un dominio h ma su un intervallo esteso. Le
formule elencate devono essere composte in modo da coprire tutto l'intervallo di integrazione. Con la regola dei
b a
trapezi avremo ad esempio, con h
N
N
b
1
xi h
f( x ) dx
a
f( x ) dx
i=0
xi
h.
f
2 0
2. f 1
2 . ......
2 . fN 1
fN
O
b
a
N
3
.N
Visto che l'intervallo di integrazione è fissato dal problema, riducendo h cresce il numero di intervalli su cui è
1
approssimato l'integrale. La regola dei trapezi converge quindi con l'ordine O
. La stessa procedura può essere
2
N
14
applicata alla regola di Simpson per calcolare una sequenza che converga con l'ordine O
1
N
4
. Applicando la regola
ripetutamente su intervalli consecutivi si ottiene
xi h
N
b
f( x ) dx
a
h.
f( x ) dx
xi
i=1
1.
f
3 1
4.
f
3 2
2.
f
3 3
....
2.
f
3 N 2
4.
f
3 N 1
1.
f
3 N
O
(b
a)
N
5
4
Possiamo scrivere a questo punto una routine che data la funzione f ( x ) , ne calcoli l'integrale. Il diagramma di flusso è
in questo caso molto semplice
TrapeziN ( f , x , N )
h
x1
x0
N
1
S
0.5 .
v
x0
N è il numero di punti su cui è calcolato l'integrale, divido quinti l'intervallo
totale per il numero di intervalli risultanti
f x1 Preparo l'integrale sommando nella variabile S il valore della funzione
agli estremi ed il valore dell'ascissa v
f x0
for i ∈ 2 .. N
Calcolo la funzione sulle ascisse nell'intervallo, distanziate di h
1
v
v
h
S
S
f( v )
Restituisco il risultato dell'integrale
S.h
Nella figura seguente è rappresentato l'errore risultante dal calcolo dell'integrale della funzione f ( x ) sin ( x )
sull'intervallo ( 0 , π ) in funzione del numero di suddivisioni dell'intervallo N 4 , 6 .. 40 . La curva continua è la funzione
2
Errore (1/N^2 è la linea continua)
N
0.2
0.1
0
0
10
20
30
40
Numero di valutazioni di f, N
Questa routine assume che il numero di suddivisioni sia scelto a priori. In pratica tuttavia quello che si desidera è di
calcolare un integrale con un certo grado di accuratezza, e di scegliere il numero di suddivisioni necessario in funzione
di questa accuratezza. A tal fine risulta utile una routine che permette di valutare l'integrale e di incrementare il
numero N senza per questo dover ricalcolare la funzione integranda in posizioni in cui il calcolo è gia stato effettuato.
La seguente routine esegue questo compito suddividendo l'intervallo in potenze di due
La routine che calcola l'integrale parziale è la seguente
15
TrapN ( f , a , N , ∆ )
v
a
j
0
0.5 . ∆
Parziale
Trasla l'inizio dell'intervallo di ∆/2
ed inizializza le variabili
0
Calcola la funzione su intervalli equispaziati ∆
while j< N
Parziale
xj
v
v
v
j
j
Parziale
f( v )
∆
1
Restituisce la somma parziale
Parziale
La seguente routine Trapezi, utilizzando la routine descritta, calcola l'integrale con la regola dei trapezi con
accuratezza ε
Trapezi ( f , x , ε )
0.5 . f x0
S
∆
x1
I
S.∆
it
1
f x1
x0
for i ∈ 0 .. 40
S
S
it
it . 2
TrapN f , x0 , it , ∆
∆
∆
2
Iv
I
S.∆
I
I
return
if I
i
2
1
Iv <ε
Piuttosto che scrivere un algoritmo analogo per la regola di Simpson, consideriamo la seguente regola di somma di
Eulero - Mclauren
b
f( x ) dx
a
h.
f
2 0
+
B2
2. f 1
. h2
2!
.
2 . ......
(1)
f
2 . fN 1
(1)
N
f
1
....
f N ...
B2 . k
( 1)
.
. h2 k
( 2 . k) !
.
( 2 .k
f
( 2 .k
1)
N
f
1)
1
I numeri Bk sono noti come numeri di Bernoulli e sono definiti in base alla seguente funzione generatrice
∞
t
t
e
1
Bn .
n=0
n
t
n!
t
moltiplicati per n ! . L'espansione ( 1 )
e 1
rappresenta l'espressione dell'errore commesso con la regola dei trapezi, in termini delle derivate successive della
funzione calcolate agli estremi dell'intervallo ( a , b ) . Osserviamo che, oltre al fatto che la serie contiene solo potenze
ovvero sono i coefficienti dello sviluppo in serie di Taylor della funzione
t
2
pari di h, all'ordine O h , il valore dell'errore ottenuto suddividendo l'intervallo in N sottointervalli, è esattamente 4
volte l'errore che si ottiene suddividendolo in 2 . N sottointervalli. Potremo quindi scrivere
16
b
b
f( x ) dx
a
SN
2
E 2.h
4
E 4. h
ed anche
....
f( x ) dx
E4
E. 2
h
4
S 2 .N
a
16
. h4
.....
dove h rappresenta il passo della prima suddivisione. Combinando le due espressioni possiamo eliminare l'errore
2
all'ordine O h
b
f( x ) dx
a
4.
S .
3 2N
1.
S
3 N
13 . E 4
12
. h4
.....
1/|Errore|
e questa non è altro che la regola di Simpson (dei 4/3) ricavata precedentemente. Nella figura seguente è
rappresentato in scala doppio-logaritmica il grafico dell'inverso dell'errore commesso sulla funzione f ( x ) sin ( x )
integrata sull'intervallo ( 0 , π ) per N 6 , 8 .. 60
1 .10
5
1 .10
4
1 .10
3
100
10
1
10
100
Numero di valutazioni di f, N
Simpson
Trapezi
Una routine che implementi il metodo di
Simpson con un controllo dell'errore è
quindi la seguente:
Simpson ( f , x , ε )
S
0.5 . f x0
∆
x1
Trap
f x1
x0
S.∆
I
Trap
it
1
for i ∈ 0 .. 40
TrapN f , x0 , it , ∆
S
S
it
it . 2
∆
∆
2
Trap v
Iv
Trap
I
S.∆
Trap
Trap . 4
I
Trap v
3
I
return
i
2
if I
1
Iv <ε
17
v
Verifichiamo questa routine calcolando
0 , 0.1 .. π . Nel grafico seguente (figura a sinistra) è
sin ( x ) d x per v
0
v
rappresentato in funzione di v il modulo della differenza
sin ( x ) d x ( 1 cos ( x ) ) dove l'integrale è calcolato
0
rispettivamente con il metodo dei Trapezi e con la routine suindicata per il metodo di Simpson, entrambe con
5
accuratezza ε 10 . Nella figura di destra è invece rappresentato il numero di valutazioni della funzione integranda
utilizzate effettivamente dalle routine che implementano il metodo dei Trapezi e di Simpson in funzione della variabile
v
4 .10
6
3 .10
6
600
2 .10
6
1 .10
6
N
Accuratezza
400
200
0
0
0
1
2
3
4
0
1
2
3
4
v
v
Trapezi
Simpson
Trapezi
Simpson
E' interessante la possibilità di estendere questo metodo anche agli ordini superiori. Dalla regola di somma di
Eulero-McLaurin abbiamo
f ( x ) d x S0
2
∆ . E1
4
∆ . E2
....
Le quantità incognite sono per noi
(N
in 2
1)
∆
2 .N .
B2 . n
.
. f( 2 n 1 )
En 1
b
.
( 2 n) !
EN
f ( x ) d x E1 E3
.. EN
e cioè N
( 2 .n
f
1)
a
1. E' necessario quindi valutare la funzione
punti per determinare queste quantità. In questo modo si determina automaticamente il valore dell'integrale
eliminando l'errore sino all'ordine 2 . N . Introduciamo u ∆2 , in forma matriciale abbiamo
1
1
2
u
u
u
u
2
2
u
2
N
2 .2
2 .N
2
.
2 .2 . 2
u
....
2
u
1
N
....
.
.
.
2
.
u
.
.
2 . n. m
.
2
2
u
1
2 .N
2
u
2 .2 . N
2
E1
.
n
.
f( x ) dx
EN
N
u
.
2 .N
S0
S2
.
.
S 2 .N
2
2
In effetti a noi non interessa calcolare i termini E1 , E2 .. EN , ma piuttosto solo l'integrale di f , ovvero il coefficiente r0 .
Possiamo quindi eliminare la dipendenza dal parametro u, ovvero dall'internallo di integrazione, ridefinendo gli errori
E1 , E2 .. EN
2
N
u . E1 , u . E 2 .. u . EN . La matrice diventa in questo caso
18
1
1
1
1
1
1
2
2
2 .N
2
2
.
2 .2 . 2
1
....
2 .2
1
1
1
....
.
.
.
.
.
1
1
1
2 .N
2
2 .2 . N
2 .N
S 2 .N
EN
1
.
2
.
.
2
O utilizzando una notazione sintetica, .U . r s dove
1
avremo per l'elemento generico di Un, m
.
( 2 ( n .m ) )
2
.
.
.
2 . n. m
S2
E1
.
2
1
S0
f( x ) dx
2
2
Per estrarre il primo elemento dalla soluzione scriviamo
1
f( x ) dx ( 1 0 . . 0 ) . U . s
r0
Supponiamo quindi che questo integrale si possa scrivere come
S0
N
w i . S 2 .i
f( x) d x
w 0 . . wN .
i=0
S2
.
S2 . N
Il vettore w deve quindi soddisfare l'identità
w 0 . . wN
(1 0 . . 0 ) .U
1
Moltiplicando a destra per U e trasponendo (U è simmetrica) si ottiene
U.
w0
1
w1
0
.
.
wN
0
e cioè il vettore dei pesi è la soluzione di un sistema lineare che può essere determinata una volta per tutte prima
dell'inizio dell'integrazione.
U è una matrice di tipo Vandermond
1 x0
x0
1 x1
2
....
x0
.
.
x1
Matrice di Vandermond: 1
.
.
.
.
1
.
.
.
.
1 xN
xN
2
....
xN
N
N
N
ed il determinante di una matrice di Vandermond è dato dalla produttoria
xi
xj . Nel nostro caso particolare,
( i< j )
xn
1
n
e la soluzione si scrive nella forma
4
19
xi
xj
( i< j )
x 0
k
wk
xi
xj
( i< j )
Per ora ci limitiamo ad utilizzare le routine di Mathcad per la risoluzione dei sistemi lineari ed a fornire il valore
1
numerico dei pesi w per N 5. Definiamo gli indici m 0 .. N ed n 0 .. N , la matrice Un, m
ed il vettore
2 . ( n. m )
2
vn
0 eccetto v0
1. Il vettore w è semplicemente w
T
w =
1.352 . 10
9
1.844 . 10
6
1
U . v . Nel nostro caso
5.017 . 10
4
0.032
0.483 1.452
La seguente una routine che valuta il vettore w che utilizzeremo per implementare l'integrazione con il metodo di
Simpson esteso.
w(N)
for i ∈ 0 .. N
vi
0
Definisco il vettore v
for j ∈ 0 .. N
Ui , j
v0
U
i .j
Definisco la matrice U
0.25
1
1
1
2 .i .j
2
v
Calcolo il vettore w
Ed ecco infine la routine che, sfruttando la precedente ed il mattone fondamentale che implementa la regola dei
trapezi, calcola la regola di Simpson all'ordine più elevato possibile con gli N punti a disposizione. Questo metodo è
noto come metodo di Romberg
Romberg ( f , x , ε )
S
0.5 . f x0
∆
x1
I
S.∆
it
1
f x1
Inizializza le variabili
x0
for i ∈ 0 .. 20
TrapN f , x0 , it , ∆
S
S
it
2 . it
∆
0.5 . ∆
S.∆
Trapi
Iv
I
Calcola la somma parziale dei valori di f
Conserva nel vettore Trap gli integrali parziali valutati con il
metodo dei Trapezi
I
Conserva il valore precedente dell'integrale per una stima
dell'errore e valuta l'integrale come prodotto scalare con il
vettore w
w ( i) . Trap
I
return
i
2
if I
1
Iv <ε
Restituisce il risultato in modo analogo alle altre routine di
esempio riortate
20
v
Verifichiamo anche questa routine calcolando
0 , 0.1 .. π . Nel grafico seguente (figura a sinistra)
sin ( x ) d x per v
0
v
è rappresentato in funzione di v il modulo della differenza
sin ( x ) d x ( 1 cos ( x ) ) dove l'integrale è calcolato
0
rispettivamente con il metodo di Simpson e con la routine suindicata per il metodo di Romberg, entrambe con
8
accuratezza ε 10 . Nella figura di destra è invece rappresentato il numero di valutazioni della funzione integranda
utilizzate effettivamente dalle routine che implementano i due metodi in funzione della variabile v
0 , 0.1 .. π
Accuratezza
1 . 10
9
1 .10
10
1 .10
11
1 .10
12
1 .10
13
1 .10
14
1 .10
15
1 .10
16
1 .10
17
150
100
N
v
50
0
1
2
3
4
0
0
1
v
2
3
4
v
Simpson
Romberg
Simpson
Romberg
Per confrontare la velocità di convergenza dei vai metodi, nella figura seguente è riportato il grafico dell'inverso
M
1/|Errore|
dell'errore commesso in funzione del numero di valutazioni della funzione N 2 , con M 1 , 2 .. 5 . Si osservi il
comportamento non lineare della curva corrispondente al metodo di Simpson Esteso, dove l'ordine di convergenza
cresce con il numero di punti su cui è valutata la funzione integranda.
1 .10
12
1 .10
11
1 .10
10
1 . 10
9
1 . 10
8
1 . 10
7
1 . 10
6
1 . 10
5
1 . 10
4
1 . 10
3
100
10
1
1
10
100
Numero di valutazioni di f, N
Simpson
Trapezi
Romberg
21
L'implementazione fornita del metodo di Romberg non è ne l'unica ne probabilmente la migliore. Se da un lato la
scelta di utilizzare il vettore dei pesi tabulati w consente di ridurre al minimo le operazioni da compiere, dall'altro ci
sono due controindicazioni intrinseche che sono, da una parte il fatto che una volta tabulati i valori w , siamo limitati nel
D
numero di suddivisioni dell'intervallo al numero 2
1 dove D è la dimensione massima di w . Dall'altra, osserviamo
che il vettore w contiene valori sempre più piccoli a mano a mano che la sua dimensione cresce, ed è moltiplicato per
numeri sempre crescenti, relativi alle somme parziali che contengono un numero sempre crescente di termini.
L'alternativa è la seguente. Riprendiamo la regola di Eulero Mclaurin
Stiamo approssimando la nostra funzione con un polinomio, e poi stiamo integrando. Il polinomio sarà in generale di
grado pari al numero di punti in cui è valutata la funzione (meno uno).
2
∆ . E1
f ( x ) d x S0
4
∆ . E2
....
∆
2 .N .
EN
ed invertiamola nel modo seguente.
S( ∆ )
2
∆ .
f( x) d x
E1
4
∆ .
E2
....
∆
2 .N .
EN
ovvero il termini relativo alle somma ottenuta con il metodo dei trapezi, è esprimibile come un polinomio in ∆2 , il cui
termine noto è l'integrale da calcolare. Conosco le somme parziali N 1 valori di ∆ ed ovviamente i valori di ∆2 .
Posso quindi utilizzare la routine che calcola il polinomio di Legendre che passa per questi punti. Se valuto questo
polinomio in ∆ 0 ottengo una routine che stima l'integrale allo stesso ordine di quella descritta precedentemente (N).
La seguente è una routine che calcola il polinomio interpolante, analoga a quella riportata nella sezione relativa
all'interpolazione, ottimizzata per il calcolo del polinomio in ∆ 0
polint ( Dx , Dy )
n
last( Dx )
P
Dyn
ns
n
C
Dy
D
Dy
1
for m ∈ 1 .. n
for i ∈ 0 .. n
dy
m
T
Ci 1 Di
Dxi Dxi m
Ci
T. Dxi
Di
T. Dxi m
C ns 1 if ( 2 . ns ) < ( n
m)
otherwise
dy Dns
P
ns
ns
P
dy
1
P
dy
22
La seguente è invece una routine che implementa il metodo di Romberg con il metodo esposto
Romberg2 ( f , x , ε )
∆
x1
St
0.5 . x1
it
1
H0
x0
x0 . f x0
Inizializza le variabili
f x1
1
for j ∈ 0 .. 20
St
0.5 . St
Sj
St
it
2 . it
∆
0.5 . ∆
∆ . TrapN
f , x0 , it , ∆
Calcola la somma parziale dei valori di f
if j> 0
SS
polint ( H , S )
return
SS0
j
2
Sj 1
Sj
Hj 1
0.25 . Hj
Calcola il polinomio interpolante
if SS1 < ε
Restituisce il risultato in modo analogo
alle altre routine di esempio riportate
1
Conserva in S e H le somme parziali ed i
valori dell'intervallo
v
Verifichiamo anche questa routine calcolando
0 , 0.1 .. π . Nel grafico seguente (figura a sinistra)
sin ( x ) d x per v
0
v
è rappresentato in funzione di v il modulo della differenza
sin ( x ) d x
(1
cos ( x ) ) dove l'integrale è calcolato
0
8
9
1 . 10
10
1 .10
11
1 .10
12
1 .10
13
1 .10
14
1 .10
15
1 .10
16
1 .10
17
1 .10
40
30
N
Errore
rispettivamente con le routines Romberg e Romberg2, entrambe con accuratezza ε 10 . Nella figura di destra è
invece rappresentato il numero di valutazioni della funzione integranda utilizzate effettivamente dalle routine che
implementano i due metodi in funzione della variabile v
20
10
0
1
2
3
4
0
0
1
v
Romberg 2
Romberg
2
3
4
v
Romberg 2
Romberg
Le due routines in questo caso presentano risultati analoghi. In generale, per una routine che deve calcolare un gran
numero di volte lo stesso integrale o integrali simili tra loro, è preferibile la prima, per una routine di integrazione di
carattere generale è invece preferibile la seconda.
Un'ultima considerazione. Ordine elevato non significa necessariamente precisione elevata. E' necessario
ricordare come nell'espansione in termini della regola di Eulero McLaurain (come del resto nello sviluppo d Taylor),
compaiano le derivate di ordine sempre più elevato con il crescere dell'ordine, calcolate agli estremi. Questo può
portare al crescere dll'ordine ad un crescita del coefficiente che moltiplica l'errore. Ordine elevato significa
accuratezza solo per quelle funzioni che possono essere approssimate efficientemente con dei polinomi.
Torna all'indice generale
23