Integrazione Numerica
•Data una funzione f integrabile su [a,b] consideriamo
b
 f (x) dx  I[f ]
a
E consideriamo il problema di valutare numericamente tale integrale.
sua approssimazione. Perché?
Non sempre è possibile esprimere la primitiva della funzione integranda in
termini di funzioni elementari
In alcuni casi non si conosce la primitiva
La funzione da integrare può essere data non in forma analitica, ma per punti.
Si cercano metodi numerici in grado di fornire una approssimazione di un
integrale in termini di un numero finito di valori della funzione integranda
FORMULE DI QUADRATURA
1
Supponiamo di conoscere (o di poter valutare) la funzione integranda f(x)
in punti {xi } (scelti o prefissati), distinti in [a,b]
Costruiamo formule del tipo
b
n
 f (x) dx   w f (x )
i
i 0
a
i
{xi } = nodi della formula di quadratura
{wi } = pesi della formula di quadratura
Dato che l’operatore integrazione è un funzionale lineare, tale formula ne
preserva, tra l’altro, questa proprieta’
En 1[f ]  I[f ]  In 1[f ]
errore di quadratura
• IDEA IMMEDIATA:
• Approssimare f(x) con il polinomio interpolante la funzione nei nodi {xi } di
grado n (unico perche’ i nodi sono distinti)
b
b
a
a
b
b
 f (x) dx   [L (f ; x)  e (f ; x)] dx   L (f ; x) dx   e (f ; x) dx
n
dove
n
n
a
n
a
Ln (f ; xi )  f (xi )
2
Formule interpolatorie
•
Se rappresentiamo Ln (f ; x) nella forma di Lagrange
Ln (f ; x) 
n
 l (x)f (x )
i 0
i

i
b
b
b
a
a i 0
n
b
n
b
 f (x) dx    l (x)f (x ) dx   e (f ; x) dx  f (x )  l (x) dx   e (f ; x) dx
i
i
n
i 0
a
i
i
n
a
a
wi
•
En 1
Le formule costruite in questo modo si chiamano formule di quadratura
interpolatorie
• ESEMPIO
Consideriamo i punti (a,f(a)) e (b,f(b)) e sostituiamo alla
funzione la retta che passa per i punti dati
b
 f (x) dx 
a
xb
l0 (x ) 

ab
xa
l1( x ) 

ba
ba
[f (b)  f (a)]
2
Regola dei Trapezi
b
1(a  b)2
ba
w0   l0 (x ) dx  

2 ab
2
a
b
w1 
 l1(x) dx 
a
ba
2
3
Grado di precisione
• La precisione di una formula di quadratura è legata alla bontà con cui In 1[f ]
approssima I[f ] , pertanto in generale è dipendente dalla funzione
integranda.
• Si esamina per quale classe di funzioni e’ esatta (cioè In 1[f ]  I[f ] )
Definizione: Una formula di quadratura ha grado di precisione k se è
esatta quando la funzione integranda è un polinomio di grado k, ed
esiste almeno un polinomio di grado k+1 per cui l’errore En 1[f ] risulti
non nullo
• Tale definizione è giustificata dal teorema di Weierstrass
• Vale il teorema seguente
Teorema: Le formule di quadratura interpolatorie costruite su n+1 nodi,
hanno grado di precisione almeno n .
en  f (x )  pn (x )   n 1(x ) f [x0 , x1,..., xn 1, x]
b
 En 1[f ]  I[f ]  In 1[f ] 
b
b
 f (x)dx   p (x)dx   
n
a
a
n 1
(x ) f [x0 , x1,..., xn 1, x]
a
tenendo presente f [x0 , x1,..., xn 1, x]  0 per f  Pn
4
Convergenza
• Dal teorema di Weierstrass discende anche il seguente
Teorema: Sia In 1[f ] una successione di formule di quadratura tali che
In 1[f ] abbia grado di precisione almeno n, ed equilimitate
(i.e.  c : In 1  c n ). Allora si ha
lim In 1[f ]  I[f ]
n 
Dim.
 p  Pn : In 1[p]  I[p] 
En 1[f ]  I[f ]  In 1[f ]  I[f ]  In 1[p]  I[p]  In 1[f ]  I[f  p]  In 1[f  p]

En 1[f ]  ( I  In 1 ) f  p  ( I  c) f  p
Per Teorema di W.
p
convergente a f  En 1[f ]  0
Teorema: Data unanfamiglia di formule di quadratura interpolatorie In 1[f ] tali che
 H :  wi  H Allora si ha
i 0
lim En 1[f ]  0
n 
5
Osservazione: In una formula di quadratura interpolatoria si ha
(dim. per esercizio)
n
w
i 0
i
ba
Corollario: Data f  C [a, b], se i pesi di una formula di quadratura di grado
di precisione n, sono tutti positivi allora
lim In 1[f ]  I[f ]
n 
6
•
Tra le formule interpolatorie piu’ usate si
possono mettere in evidenza due classi
importanti:
1.
Formule di Newton-Cotes: i nodi x0 , x1,..., xn sono prefissati
nell’intervallo [a,b] e sono equispaziati. Queste formule hanno
grado di precisione n o n+1 ed hanno i pesi wi facilmente
ricavabili ed espressi con semplici numeri razionali. Hanno però lo
svantaggio che per n>9 i
non sono tutti dello stesso segno.
wi
2.
Formule gaussiane: i nodi x0 , x1,..., xn non sono prefissati a priori,
ma assime ai pesi vengono ricavati in modo da massimizzare il
grado di precisione che risulta di 2n+1. Queste formule, rispetto a
quelle di Newton-cotes, hanno il vantaggio, oltre all’elevato grado di
precisione, di avere i pesi sempre positivi ,al prezzo pero’ che
l’espressione dei nodi e dei pesi e’ spesso non razionale.
7
Formule di Newton-cotes
• Basate sul metodo di interpolazione di Lagrange con nodi equispaziati in
[a,b]:
xi  x0  i h ,
i  1,..., n ;
ba
n
se
x0  a e xn  b  h 
FORMULE DI TIPO CHIUSO
se
x0  a  h e xn  b  h  h 
ba
n2
FORMULE DI TIPO APERTO
• I pesi, che dipendono solo da n e h (per tutte le formule interpolatorie),
non dipendono dall’intervallo di integrazione [a,b] e sono esprimibili
come
b
wi   l dx 
a
( n)
i
xn ( xn  h)

n

x0 ( x0  h) j  0, j  i
x  xj
xi  x j
n( n 1)
dx 

0( 1)
t j
h dt  h

j  0, j  i i  j
n
n( n 1)

0( 1)
t j
dt  h  i
j  0, j  i i  j
n

i
8
Esempi
• n=2 chiusa: SIMPSON
ba ab
ba

, x2  b  h 
2
2
2
h
4h
h
w0  , w1 
, w2  ,
3
3
3
x0  a , x1  a 
b
 f (x) dx I3[f ] 
a
h
ab

f
(
a
)

4
f
(
)

f
(
b
)

3 
2

• n=0 aperta: MIDPOINT
x0  a  h ,  h 
ba
2
w0  b  a
b
 f (x) dx I1[f ]  (b  a)f (
a
ab
)
2
9
Errore per formule di N.-C.
• Teorema: Sia f  C n 2[a, b]con n pari, allora si ha
f (n  2)( ) n  3
En 1[f ] 
h
(n  2)!
Se invece n e’ dispari,
n 1  i
n
i 1
j 0
 t  (t  j) dt
f  C n 1[a, b]
f (n 1)( ) n  2
En 1[f ] 
h
(n  1)!
n 1  i
n
i 1
j 0
 t  (t  j) dt
dove i=1 per le formule di tipo chiuso, mentre i=0 per quelle
di tipo aperto.
•
•
A parita’ di numero di nodi formule aperte e formule chiuse hanno lo stesso
grado di precisione.
Per n pari (numero dispari di nodi) il grado di precisione e’ n+1 e non n
(come invece si ha per n dispari)
E’ piu’ conveniente usare formule con n pari.
10
Convergenza formule di N.-C.
• Contrariamente a quanto potrebbe sembrare “a prima vista” non
conviene usare formule di Newton-Cotes di grado di precisione via
via crescente.
• Teorema (Kusmin): Per ogni successione In 1[f ] di formule di
quadratura interpolatorie costruite su un intervallo chiuso con nodi
equidistanti si ha
n
lim  wi(n 1)  
n 
i 0
• I pesi tendono a crescere in modulo e ad essere di segno alterno,
dando luogo a rilevanti errori di arrotondamento (per es. errori di
cancellazione)
11
Condizionamento
• In generale la quantita’  wi da’ una misura di quanto si
amplifichino gli errori sui dati iniziali e quindi puo’ essere messa
anche in relazione con il condizionamento del problema.
• Supponiamo che a causa degli errori di arrodontamento si abbia:
~
~
f j  f (x j )  fl (f (x j ))  f (x j )   j
~
I n 1[f ]  In 1[f ] 
n
n
 w (f (x )   )   w f (x )
i 0
i
i
i
i
i 0
i

n
w 
i 0
i
n
i

   wi ,
i 0
• L’errore effettivamente commesso e’ dato da
~
n
~
E n 1  I[f ]  I n 1[f ]  I[f ]   wi f (xi ) 
i 0
n
• Se I pesi sono positivi:
w 
i 0
i
n
w
i 0
i
i
n
w
i 0
i
i
  max  i
 En 1[f ] 
n
w
i 0
i
i
n
i
   wi   (b  a)
i 0
12
Formule composite
• Idea: si suddivide l’intervallo di integrazione [a,b] in N sottointervalli
di ampiezza uguale h 
formula di grado basso
N
b
 f (x)dx 
a
N  1 z k 1
 
f (x )dx
k 0 z
k
ba
N
N 1
I
k 0
e su ciascuno di essi applicare una
(k )
n 1
[f ],
z k : k  0,..., N; z0  a, z N  b
• Le formule piu’ usate:
N 1
z k 1  z k
ba
I [f ]  
(f (zk )  f  zk 1 ) 
[f (a)   2f (zk )  f (b)] TRAPEZI
2
2
N
k 0
k 1
N 1
N 1
N 1
z k 1  z k 1
z  zk
4 z k 1  z k
1
ba
N
I3 [f ]  
( f (zk )  f (
)  f  zk 1 ) 
[f (a)  2  f (zk )  4  f ( k 1
)  f (b)]
2
3
3
2
3
6N
2
k 0
k 1
k 0
N
2
N 1
SIMPSON
Il grado di precisione delle formule composite e’ lo stesso delle
corrispondenti formule di Newton-Cotes “semplici”.
13
Errore
Si puo’ facilmente dimostrare:
N
2
E
(b  a)3 (2)

f
2
N
[a, b]
1
12
5
N
3
E
 b  a
(4)

N
f

2
N


[a, b]
1
90
Per funzioni sufficientemente regolari si ha quindi:
lim EnN1[f ]  0
N 
14
Implementazione
• In pratica e’ importante determinare un valore adeguato del numero
di suddivisioni dell’intervallo che bisogna fare.
• Stima dell’errore in modo automatico valutando:
InN11[f ]  InN21[f ]
• Di solito e’ conveniente considerare N2=2N1
• Nell’ipotesi in cui la derivata s-esima vari lentamente e’ possibile
dare una stima dell’errore a partire da due formule composite con
valori diversi di N
s 1
s
(b  a)
1 f ( )  N
ns 1
Ns s!
Ns
 2N
(b  a)s 1 1 f s ( )
2N
I[f ]  I [f ]  kn
n s 1
(2N )s s !
(2N )s
I[f ]  I N [f ]  kn

2N
N
I [f ]  I [f ]

s
2 N
s
s
(2  1)

2N
E [f ]

I 2N [f ]  I N [f ]

s
s
2 N
2s  1
Estrapolazione di Richardson
15
Formule adattive
• Le formule composite con suddivisione uniforme dell’intervallo di
integrazione sono ormai superate, tranne in casi particolari (funzioni
periodiche, regola dei trapezi)
• Si usano formule di tipo adattivo
• Quando la funzione integranda presenta delle irregolarita’ c’e’ la
necessita’ di addensare I nodi nelle vicinanze delle irregolarita’
• L’intervallo viene suddiviso in sottointervalli di ampiezza diversa
• Viene applicata una formula base (non con molti nodi) addensando I
nodi la’ dove e’ necessario.
16
Polinomi ortogonali
•
Uno spazio vettoriale G sul campo R dei reali, su cui è definito un prodotto
scalare, è detto spazio di Hilbert se ogni successione di Cauchy di elementi di G
è convergente ( spazio di Banach + prodotto scalare)
• Data una funzione peso w(x) non negativa nell’intervallo (finito o infinito)
(a,b) e non identicamente nulla ed un insieme pi (x)iN di polinomi, in
cui pi (x) è di grado i, esso è uno spazio di Hilbert rispetto al prodotto
scalare
b
pi , p j   w(x)pi (x)p j (x)dx
a
• Un sistema di polinomi è detto ortogonale rispetto ad una funzione peso
w(x) ed al prodotto scalare sopra definito se
pi , p j  0
per i  j
• L’intervallo (a,b) e la funzione peso w(x) definiscono univocamente i
polinomi pi (x), a meno di fattori costanti non nulli.
17
• Teorema: Sia q(x) un polinomio di grado n. Allora
pi , q  0
per i  n
• Teorema: per ogni n>0 il polinomio ortogonale pn(x) possiede n zeri
reali distinti e tutti contenuti in (a,b). Inoltre gli zeri di pn(x) si
alternano con quelli di pn 1(x) ossia tra due zeri consecutivi di pn (x)
esiste un solo zero di pn 1(x)
• Teorema: Ogni sistema di polinomi ortogonali pi (x)iN soddisfa
una relazione per ricorrenza a tre termini del tipo:
pi 1(x)  ( Ai x  Bi )pi (x)  Ci pi 1(x)
18
Formule Gaussiane
• Nelle formule gaussiane i nodi non sono determinati a priori, ma
sono scelti in modo da massimizzare il grado di precisione della
formula
 Teorema: Una formula di quadratura di tipo interpolatorio costruita
su n+1 punti ha grado di precisione almeno n ed al massimo 2n+1.
Il grado massimo viene raggiunto se e solo se i nodi sono gli zeri
dell’ (n+1)-esimo polinomio ortogonale rispetto alla funzione peso
w(x)=1
b
Si considera il resto nella forma di Newton En 1[f ]   n 1(x) f [x0 , x1,...xn , x] dx
a
Per costruzione se f  Pk con k  n  f [x0 , x1,...xn , x]  0
Vediamo ora f  Pn  r con r  0  f [x0 , x1,...xn , x] Pn  r (n 1) r 1
b

a
b
n 1
(x) pr 1 dx  0   n 1(x) x s dx  0 per s  0,1,...r  1
a
se e solo se n 1e’ ortogonale allo spazio Pr 1 rispetto a
b
f , g   f (x)g(x) dx
a
E’ sufficiente prendere come nodi gli zeri dell’(n+1)-esimo polinomio
ortogonale r-1=n
r=n+1. Il massimo grado di precisione ottenibile e’
2n+1
19
 Tutti i coefficienti di una formula interpolatoria In 1 di grado di
precisione almeno 2n, sono tutti positivi
 Per quanto riguarda l’errore si ha che, se
f  C 2n 2(a, b)
b
f (2n 2)( ) 2
En 1[f ] 
qn 1(x)dx

(2n  2)! a
dove qn 1 e’ l’(n+1)-esimo polinomio ortogonale
 Le formule gaussiane classiche sono quelle associate ai classici
polinomi ortogonali.
 Essi sono in generale definiti in [-1,1]. Se [a, b]  [1,1]
x
1
[(b  a)  t (b  a)]
2
20
Formule di Gauss-Legendre
• I polinomi di Legendre sono definiti per ricorrenza

P0 (x )  0

P1(x )  x

2n  1
n
Pn 1(x ) 
x Pn (x ) 
Pn 1(x )

n 1
n 1
• Per queste formule si possono calcolare esplicitamente i nodi ed i
pesi, per ogni n, e tabularli (nell’intervallo di riferimento [-1,1])
• Difficile implementazione per loro utilizzo in modo iterativo
• I nodi sono tutti interni all’intervallo
21
Formule di Gauss-Lobatto
• Si considerano gli zeri del polinomio  n 1  Pn 1(x)  aPn(x)  bPn 1(x)
dove a e b sono determinate in modo tale che
 n 1(1)   n 1(1)  0
• Nei nodi sono inclusi gli estremi
G L
n 1
I
n 1
2
[f ] 
(f (1)  f (1))   w j f (x j )
n(n  1)
j 1
'
• I nodi sono gli zeri di Pn 1( x )ed i pesi sono dati da
wj 
2
(Pn (x j ))2
n(n  1)
• Le formule di Gauss-Lobatto hanno grado di precisione 2n-1
22
Formule gaussiane pesate
b
• Si studiano formule per approssimare
 w(x) f (x) dx
a
• Questi integrali sono utili quando si vuole integrare una funzione g(x)
che presenta delle singolarità o dei punti di discontinuità nell’intervallo
[a,b], ma può essere fattorizzata g(x)=w(x)f(x), in cui w(x) è una
funzione di forma semplice contenente la singolarità di g(x), mentre f(x)
è la parte regolare di g(x).
• Si costruiscono formule interpolatorie generalizzando il procedimento
visto finora (che corrisponde a w(x)  1)
b
b
n
In 1[f ]   w( x ) pn (x )dx   w(x ) l
a
a
i 0
( n)
i
(x )f (xi ) dx 
n
 w f (x )
i 0
i
i
b
wi   w(x )li(n)(x )dx
a
• La funzione peso deve essere tale da garantire l’esistenza degli
integrali coinvolti e permettere la costruzione dei pesi
23
• Formule gaussiane pesate in cui i nodi sono scelti coincidenti con gli
zeri dei polinomi ortogonali rispetto al prodotto scalare
b
f,g
w
  w(x)f (x)g(x)dx
a
• Particolarizzando l’intervallo [a,b] e la funzione peso w(x) si hanno
diverse formule:
 Gauss-Chebychev:
[-1,1]
w(x ) 
 Gauss-Laguerre: [0, ), w(x)  e
 Gauss-Hermite :
(, )
x
1
1 x
2
;
T0 ( x)  1

T1(x )  x
T (x )  2 x T ( x)  T ( x)
n
n 1
 n 1
L0 ( x)  1

L1(x )  x
(n  1)L (x )  (2n  1  x) L (x )  nL (x )
n 1
n
n 1

w(x)  e  x
2
H0 (x )  1

H1(x)  2 x
H (x )  2 x H (x )  2nH (x )
n
n 1
 n 1
 Tutti i nodi ed i pesi sono ottenibili da opportune tavole, tranne che per
le formule di Gauss-Chebychev:
x j  cos(
2j  1
),
2n  2
j  0,.., n;
wj 

n 1
24