INTEGRAZIONE NUMERICA

Integrazione
Corso: Analisi Numerica
Anno Accademico: 2004-2005
INTEGRAZIONE NUMERICA
Le procedure numeriche per approssimare l’integrale
definito:
b
I (f )   f (x )dx
a
Date da:
Qn f  
n
 a f x 
i 0
i
i
Sono note come formule di quadratura numerica.
[a,b] è un intervallo chiuso e limitato.
Gli n+1 punti distinti xi sono i nodi e gli ai sono i
pesi della quadratura.
Il problema è determinare xi ed ai in modo che
Q ( f ) approssimi I ( f ) un’ ampia classe di funzioni.
Se pn (x) Pn è un polinomio interpolante la f(x)
negli xi
la formula:
Qn f
n
   ai f x i   a pn x  dx
i o
b
si dice formula di quadratura interpolatoria.
I nodi e i pesi sono scelti in modo da minimizzare
l’errore:
E n (f)  I(f)  Qn (f)
Una misura di tale errore è dato dal grado di
precisione.
Un modo pratico di calcolarlo è determinare una classe
di funzioni per la quale la formula risulti esatta.
Generalmente tale classe è quella dei polinomi per cui
una formula si dice esatta di grado k se risulta esatta
per p P .
k
Un modo generale per costruire formule di quadratura
con grado di precisione fissato è il metodo dei
coefficienti
indeterminati,
che
consiste
nel
determinare i nodi e i pesi imponendo che la formula
sia esatta per polinomi del grado dato dalla precisione.
Se i nodi sono fissati, i pesi si trovano
risolvendo il sistema lineare:
n
 ai xi 
i0
r

b
a
x r dx
0  r  n
Se i nodi non sono fissati, il sistema è non
lineare, e ciò vedremo che darà luogo alle
formule col più alto grado di precisione possibile.
FORMULE DI QUADRATURA INTERPOLATORIE
Siano xi  x j , i  j,
i,j  0,..,n,
n  1 punti di interpolazione e costruiamo
il polinomio di interpolazione
pn (x)  Pn
per f(x)  C a,b  ovvero tale che :
f(x j )  pn (x j ),
j  0,..,n.
Qn (f) 
pn (x) 

b
a
n
l
j 0
nj
pn  x dx
( x )f ( x j )
Qn ( f ) 
anj 
lnj 

b
a
n
 a f x 
nj
j 0
j
lnj  x  dx
Wn ( x )
 x  x W  x 
'
n
j
n

Wn (x)   x  x j
i 0
j

j  0,...,n
Ogni formula di quadratura interpolatoria che usi n+1
nodi ha, per costruzione, grado di precisione almeno n.
Le formule più naturali sono quelle con i nodi
ugualmente spaziati in [a,b].
Tali formule sono le formule di NEWTON-COTES.
Sia:
b a
h
,
x j  x 0  jh,
n
x 0  a,
xn  b.
FORMULA DEL TRAPEZIO
La formula di NEWTON - COTES a due punti in cui :
x0  a, x1  b è detta formula del trapezio.
Ricaviamola per il generico intervallo -h, h  e poi per a, b 
x0  h, x1  h.
Il polinomio p1 (x)  P1 tale che : p(xi )  f(xi ),
i  0,1
h-x
h+x
è dato da : p1 (x)  f(-h)
 f(h)
2h
2h

h
QT (f)  I(p1 )   p1  x dx  hf  h   hf h 
h
Per ricavarlo per a,b  usiamo il
metodo dei coefficienti indeterminati.
1
QT ( f )   ai f ( xi )  a0f ( x 0 )  a1f ( x1 )
i o
imponiamo che il grado di precisione sia 1 e sia :
f ( x )  1, x
1
a
i 0
 a0  a1 
i

a
1
 a f ( x )  aa
i 0
da cui:
pertanto:
i
b
0
xdx  b  a
 ba1 

b
a
b2  a 2
xdx 
2
b a
a0  a1 
2
b-a
QT (f ) 
f ( a )  f ( b ) 
2
Geometricamente:
errore
f(b)
f(a)
a
b
Per ricavare l'errore ricordiamo che se f  C n 1 [a,b]
l'errore dell'interpolazione è :
f ( n 1 )  
e(x)  f(x)-pn ( x ) 
W( x)
n  1 !
n
W(x)  ( x  xi )
i 0
b
b
a
a
per cui : en   f ( x )  pn ( x ) dx  
ponendo :
si ha :
f ( n 1 )  
W ( x )dx
n  1  !
Mn 1  max f ( n 1 ) ( x )
en  Mn 1 
b
a
b
W( x)
dx
( n  1) !
eT   f  
a
''
x  a  x  b 
2
dx
Applichiamo ora il teorema del valor medio
sugli integrali per il quale , se g, h  C a,b  e g(x)
non cambia segno in a,b 
 n  a,b :

b
a
b
g(x)h  x   h  x   g  x dx
a
ponendo : g(x)   x-a   x  b  ,
h(x)  f '' x 
si ha :
f '' (n)
eT 
2

b
a
f '' (n)
3
x  a  x  b   
b  a 
12
Si può verificare che il grado di precisione è 1.
REGOLA DI SIMPSON
La formula di Newton-Cotes a 3 punti è detta
regola diSimpson.
Poniamo : x0  h, x1  0, x2  h
ed imponiamo che :
Per:
2
 a f(x ) 
Q2 (f) 
i 0
i
i
f(x)  1, x, x 2
a0  a1  a2 

h

h
h
-a0 h  a2 h 
-h
1dx  2h
xdx  0
2 3
a0 h  a2 h   x dx  h
h
3
2
2
h
2
h
h
f ( x ) dx
h
4
 a0  a2  , a1  h
3
3
h
Q2 (f)  f (-h)  4 f(0)  f(h) 
3

b-a 
a  b 
e per a, b  : Q2 (f) 
f (a)  4 f 
  f  b 

6 
 2 

Si può facilmente verificare che il grado di
precisione è 3 e ciò è sfruttato per determinare
l’errore.
Infatti, poiché  x-x0   x-x1  x-x2  cambia segno
in [a,b] non si può procedere come prima.
Si definisce invece il polinomio hermitiano p3 (x)  P3
con le seguenti condizioni:
 a+b2    a+b2 =f  a+b2  ,
' a+b
p'  a+b

f
2 
2 
p3 (a)  f (a), p3
p3 (b)  f(b)
e poichè il grado di precisione è 3 :
f IV  
f  x   p3  x  
4!
 x-x0   x-x1   x-x2 
2
a cui può applicarsi il teorema del valore medio.
Da cui:
eS
-f IV ( x )

90
 
b-a 5
2
L’errore dell’integrazione delle formule di Newton-Cotes
ha ordine 2n+1 se i nodi sono n+1 , mentre si può fare
vedere che la precisione dipende da n.
In particolare se:

n dispari
precisione n

n pari
precisione n+1
Esempi:
Trapezio :
Simpson:
Generale:
2 nodi  n  1, eh 3 , prec. = 1
3 nodi  n  2,
n + 1 nodi ,
e h
2n 1
eh 5 , prec. = 3
, prec.
n
n +1
dispari
pari
Per aumentare la precisione si hanno 2 alternative:
i)Aumentare il numero di nodi in modo che Qn (f)
sia integrale di un polinomio interpolante di alto grado:
Quadrature Gaussiane
ii)Si divide [a,b] in sottointervalli, in essi si usano
formule di bassa precisione, si sommano i risultati:
Regole di Quadratura Composte.
Esaminiamo prima le quadrature composte
Regole di Quadratura Composte
Suddividiamo [a,b] in n intervallini:
Qn (f) 
n-1

j 0
x j 1
xj
f(x) dx
e usiamo in [x j ,x j 1 ] la regola del trapezio.
b-a
Sia : h 
, x j  a  jh, j  0,...., n
n
f (n j ) 3
x j 1
h
x j f(x)dx  2 f(x j )  f(x j 1 )   12 h
Sommando si ha :
n f (η )
h
j
Tn (f)  h f(x j )  f(x0 )  f(xn )   
h3
2
12
j 1
j 0
n 1
Per semplificare l’espressione dell’errore usiamo il
lemma:
Sia g(x)  C a,b  e
a 
n 1
j
j 0
a j   tutte dello stesso segno
x j  a,b  , j  0,...n-1
n 1
a
j 0
n-1
j
g(x j )  g    a j
i 0
  η  a,b :
h3
Identificando f (η ) con g(x) e a j con si ha :
12
n
h3
h3
et  -f     j 12  f    n 12  f    h 2  b12a 
 
e indicando con Σ  la sommatoria dimezzata agli estremi si ha :
n
Tn (f)  h ''f(x j )  f   
j 0
h 2 ( b a )
12
Nelle formule di Newton – Cotes il calcolo dei pesi è
indipendente dalla spaziatura h ed essi possono essere
quindi tabulati.
Si può vedere che, per n grande, i pesi aumentano di
modulo mentre il segno varia. Ciò rende instabili tali
formule dal punto di vista della propagazione degli
errori, inoltre un aumento del grado di precisione,
ovvero dei nodi della quadratura, non implica
necessariamente la convergenza della quadratura
all’integrale quando la funzione non è polinomiale.
Il seguente teorema mostra sotto quali condizioni
l’aumento dei punti di interpolazione porti alla
convergenza della quadratura all’integrale.
Teorema
Sia f  C a, b , Qn (f) 
n
 
(n)
(n)
a
f
x
 j
j
j 0
dove a j(n) , x j(n) sono i pesi e i nodi della quadratura
interpolat a dipendenti da n.
Se  K  0 :
n

j 0
a j(n)  K n  N
 lim Qn f   I(f) f  C a, b 
n 
Dim.:
Per Weierstrass ε  0
f-qN
 qN (x)  PN (N  f    ) :

ε
Poichè la quadratura è interpolatoria :
Qn (qN )  I(qN ) n  N
Scegliendo n  N si ha :
I(f) - Qn (f)  I(f)  I(qN )  Qn (qN )  Qn (f)  I(f)  I(qN )  Qn (qN )  Qn (f)
I(f)  I(qN )  f  qN
Qn (qN )  Qn (f) 
n
a
j 0
(n)
j

(b  a)
qN (x

(n)
j
)  f(x
_
 I(f)  Qn (f)  k  b  a  ε   ε
Si può provare che è vero il viceversa
(n)
j
)   f  qN
n


j 0
a j(n)  Kε
Se gli a j(n) sono tutti  0 la convergenza è garantita.
Infatti , poichè il polinomio p0 (x)  1 è integrato esattamente
si ha :
b
0  I 1  =  dx  Qn ( 1 ) 
a
n
(n)
a
 j
j 0
Quindi se i pesi sono tutti positivi :
n
a
j 0
(n)
j

n
a
j 0
(n)
j


b
a
dx
Un vantaggio delle formule con pesi positivi è che hanno
buone proprietà di arrotondamento poiché gli errori
tendono a cancellarsi. Inoltre l’errore è minimizzato se
i pesi sono quasi uguali. Un’idea è allora di determinare
formule con pesi uguali e nodi determinati imponendo
che la formula abbia grado di precisione n.
Si ha :
n
Qn (f)  an  f(x j )
j 0
Perchè integri esattamente f(x)  1
si deve imporre an 
Per n  1: a1  b  a
Q1 (f)  (b  a) f
 a+b2 
b a
n
x1  21 (b  a) da cui :
MIDPOINT
RULE
Metodo Midpoint
Integrazione esatta di un’approsimazione lineare di
Taylor dell’integranda.
I f



b
a
f x dx
Approssimazione lineare di Taylor ad f(x) in
a  b
c 
2
p1  x   f c    x  c f ' c 
QMP f  
poichè :


b
a
b
a
p1  x dx  b  a f c 
 x  c   f c dx  0
'
Ricaviamo l’errore
2
f  x   f c    x  c f c   x  c  f ''  
'
EMP f   
b
1
a 2
2 ''
 x  c  f   
1
2
1
24
b  a 
3
che è la metà dell’errore del metodo dei trapezi.
f  
''
Formule di questo tipo hanno lo svantaggio di dover
trovare le radici di polinomi di grado crescente.
Vediamo ora le formule di Quadratura Gaussiana in cui
i nodi che i pesi sono indeterminati.
Formule di QUADRATURA GAUßIANA
Risolvendo il SISTEMA NON LINEARE:
Qn (f) 
n
 a f(x
i 0
i
i
)
in cui sia ai che xi siano INDETERMINATI e
imponendo che la formula abbia precisione 2n+1 se
n+1 sono i nodi della Quadratura, si ottiene la
quadratura di tipo Gaussiano.
Il sistema risultante avrà 2n+2 incognite.
Per n =0 e [a,b] =[-1,1] :
I(f ) 

1
1
f ( x )dx
Q0 ( f )  a0f ( x 0 )
I ( f )  Q0 ( f )  E 0 (  )
imponendo che


1
1
1
1
E0(f) = 0
1dx  2  a0
xdx  0
per
f(x) =1, x
a0  2, x 0  0

Q0 ( f )  2f ( 0 )
si ha:
che per [a,b] generico dà:
Q0( f )  ( b  a )f (
a b
2
)
che è la regola del punto di mezzo, che quindi è di
tipo Gaussiano.
Per n=1:
a0  a1  2
a x  a x  0
a0  a1  1
1 1
 0 0


2
2
2
3
a
x

a
x

x


1 1
3
0
3 , x1 
 0 0
a x 3  a x 3  0
1 1
 0 0
3
3
Notiamo che tale formula ha grado di precisione 3 e
usa 2 punti mentre la regola di Simpson per avere la
stessa precisione usa 3 punti.
Quindi, in generale, si deve risolvere il sistema non
lineare:
n
 ai x i
i 0
r


b
a
x r dx
r  0 ,...,2 n  1
nelle 2n+2 incognite a0, …, an, x0, …, xn
Però, nell’ambito delle formule di Quadratura
Interporlatorie si può trovare un’opportuna formula
per Qn(f) con grado di precisione 2n+1, che, per n+1
nodi, è il max possibile, quando si conoscono gli n+1 nodi
senza dover risolvere il sistema non lineare.
A tale scopo si ha: TEOREMA
Se Qn (f) 
n
a
i 0
i
f ( xi )
è una Formula di Quadratura di tipo
INTERPOLATORIO, ovvero:
b
Qn ( f )   pn ( x )dx
a
dove pn(x) n è in un polinomio interpolante f(x) negli n+1
nodi: x0,…,xn e tali nodi sono gli zeri di un polinomio
pn+1Tn+1 insieme dei polinomi ortogonali su [a,b], allora il
grado di precisione della formula è 2n+1.
Dimostrazione
Sia:

b
a
f ( x ) dx 
n
af( x
i 0
i
i
)  En ( f )
Sia f(x) P2n+1 e dividiamolo per pn+1(x) dell’enunciato:
f(x)=pn+1 (x) q(x)+r(x)
dove q(x) ed r(x) sono polinomi al più di grado n.
Poiché gli xi sono gli zeri di pn+1(x) si ha:
f(xi)=r(xi)
i=0,…,n
pertanto:

b
a
b
b
a
a
f ( x ) dx   pn  1 ( x ) q ( x ) dx   r ( x ) dx 
n
 a r( x
i 0
i
i
)  En ( f )
essendo la formula di tipo interpolatorio, essa ha almeno
precisione n.

b
a
r ( x ) dx 
n
 a r( x
i 0
i
i
)
e poiché q(x)pn+1(x):

b
a
pn  1 ( x ) q ( x ) dx  0  E n ( f )  0
avendo imposto f P2n+1 la formula ha precisione 2n+1
Mostriamo ora che le formule Gaussiane hanno i pesi
positivi.
Se:
n
af( x
Qn ( f ) 
i 0
i
i
)
è Gaussiana, ha precisione 2n+1 e come f(x) prendiamo il
quadrato dei polinomi di Lagrange:
lk ( x ) 
n

i  0 ,i k
x xi
xk xi

, 0kn
(lk(x))2P2n e poiché: lk(xi)=ik si ha:
0 

b
a
(lk (x)) 2 dx  ak
0kn, c.v.d.
Calcolo dei Nodi e dei Pesi ( QUADRATURA)
Per calcolare i nodi di una quadratura Gaussiana si
procede nel seguente modo:
 Si generano prima i polinomi ortogonali usando le
formule di ricorrenza. Poiché gli zeri di tali polinomi
sono semplici reali ed interni all’intervallo di
ortogonalità si può usare il metodo di NEWTON per
determinarli.

Per calcolare i pesi invece si possono usare:
1. IL metodo dei Coefficienti Indeterminati oppure
2. Si ricavano da a j 

b
a
lnj2( x )dx con
0  j  n dove
lnj sono i polinomi di Lagrange di grado n. Se l’integrale
da calcolare è del tipo: I ( f ) 

b
a
f ( x ) ( x ) dx
b
 ( x)dx  0
dove (x) è una funzione peso tale che:
a
allora la QUADRATURA cioè i nodi e i pesi dipendono da (x).
In tal caso si scelgono i polinomi ortonormali in [a,b] rispetto
ad (x).
Se ( x )  1  x  p 1  x  q in [-1,1] con p<1, q<1 i polinomi
sono quelli di JACOBI.
Se invece
1
p  q 
2
ovvero
( x ) 
1
(1  x2 )
i polinomi sono
quelli di CHEBICHEV. Con tali polinomi i coefficienti sono

uniformi e per n nodi sono dati da: n cioè:
n 1

b
a
Se
 ( x)  1
f( x)
1 x
2
dx 

n
f ( x
i 0
i
)
i polinomi sono quelli di LEGENDRE.
Metodi di Estrapolazione
Servono per prendere informazioni da poche
Approssimazioni e usarle sia per stimare l’errore che per
avere un’approssimazione migliore.
Supponiamo che si abbia:
I f

I n f

 cn  p
P=2 Trapezi, Midpoint
P=4 Simpson
non valida per Gauss
Servono per stimare p , l’errore, e migliorare
l’approssimazione
Estrapolazione di Richardson

Con tale procedimento è possibile ottenere da una
formula di basso ordine di accuratezza una formula di
accuratezza maggiore.
Sia Q(h) una formula con accuratezza p ovvero:
I ( f )  Q(h)  Ch  o(h )
p

(1°)
p
Dove o(h p ) è un infinitesimo di ordine superiore a p
usando un passo qh si ha:
I ( f )  Q(qh)  C (qh)  o(h )
p
p
(2°)
Moltiplicando la (1°) per qp e sottraendo la (2°) si ottiene:
q Q(h)  Q(qh)
p
I( f ) 
 o( h )
p
q 1
p
Che per tanto ha un ordine più elevato. Se la formula di
partenza ammette uno sviluppo dell’errore del tipo:
I ( f )  Q(h)  c1h p1  c2 h p2  ....  ck h pk  ...
Si ha:
Con l’errore
q pk Qk (h)  Qk (qh)
Qk 1 (h) 
q pk  1
 O h
pk 1

Stima di p
Supponiamo di avere n, 2n, 4n punti e applichiamo la :
I f   I n f   cn  p

I f   I 2n f   c( 2n )  p
I f   I 4n f   c( 4n )  p
Consideriamo:
r4n 

r4n 
I n  I   I
I 2n  I   I
2n  p  n  p
4n  p  2n  p
I n  I 2n
I 2n  I 4n
 cn  p  c2n  p
 I 2n 


 I 4 n   c  2n   p  c ( 4 n )  p
2 p  1
p
 p

2
4  2 p

r4n  2 p
lg r4n
p 
lg2
Che può essere usata sia per verificare se il programma
lavora correttamente, sia per stimare la rapidità di
Convergenza quando l’integranda non è così regolare da
poter applicare la teoria dell’errore.
Per stimare l’errore si ha:
 
I  I 2n  c2n  p  2  p cn  p  2  p I  I n 
2 p I 2n  I n
 I 
2p  1
R2n
Inoltre
2 p I 2n  I n

2p  1
E 2n  R2n  I 2n
I n  I 2n

2p  1
Integrazione di Romberg

Tale metodo si ottiene applicando l’estrapolazione
di
Richardson al metodo dei trapezi. In tale metodo l’errore
è:
E h   h 2
L’integrale è:
I(f )  T( h)  E( h)

Allora per due valori h1e h2 si ha:
I ( f )  T ( h1 )  E ( h1 )  T ( h2 )  E ( h2 )
e poiché:
E h1 
E h2 
h

h
2
1
2
2
si ha:
h12
E h1   E h2  2
h2
Quindi:
h12
I ( h1 )  E ( h2 ) 2  I ( h2 )  E ( h2 )
h2
Da cui si ricava che:
E ( h2 ) 
Sostituendo:
I ( h1 )  I ( h2 )
2
 h1 
1

h
 2 
I ( f )  I ( h2 )  E ( h2 )  I ( h2 ) 
I ( h1 )  I ( h2 )
2
 h1 
1 
 h2 
Siano
h2
h1

2
:
I ( h2 )  I ( h1 )
4
1
I ( f )  I ( h2 ) 
 I h2   I ( h1 )
4 1
3
3
4
I 
I h2
3

1

I h1 
3
h2 
h1
2
 0
Quindi indicando con Tn
il metodo dei trapezi si
ha:
4T  0  T  0 
1 
T2n 
T2n
 j  1
Tn ( 0 )
T2n ( 0 )
2n
3
n
4 j  1T2n  j   Tn  j 

4 j 1  1
T2n (1)
Trapezi
T4n ( 1 )
T2k n ( 0 )
e l’errore :
T2 K n ( 1 )
T2 k n ( k )
E 2n ( k ) 
Tn ( k )  T2n ( k )
4k  1