Integrazione

Corso: Analisi NumericaAnno Accademico: 2004-2005

Le procedure numeriche per approssimare l’integrale definito:

Date da:

Sono note come formule di quadratura numerica.[a,b] è un intervallo chiuso e limitato. Gli n+1 punti distinti sono i nodi e gli sono i pesi della quadratura.Il problema è determinare ed in modo che un’ ampia classe di funzioni.

INTEGRAZIONE NUMERICA

b

adxxffI )()(

ix

ix

ia

ia

in

iin xfafQ

0

)f(Iapprossimi)f(Q

Se è un polinomio interpolante la f(x) negli

la formula:

si dice formula di quadratura interpolatoria.I nodi e i pesi sono scelti in modo da minimizzare l’errore:

n np (x) P

ix

dx xpxfafQb

a n

n

oiiin

(f)QI(f)(f)E nn

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 . 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.

k

p P

Se i nodi sono fissati, i pesi si trovano risolvendo il sistema lineare:

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.

n

0i

riixa

b

a

rdxx nr0

FORMULE DI QUADRATURA INTERPOLATORIE

Siano

punti di interpolazione e costruiamo

il polinomio di interpolazione

per ovvero tale che :

i j

n n

j n j

x x , i j , i,j 0,..,n,

n 1

p (x) P f (x) C a,b

f (x ) p (x ), j 0,..,n.

0

b

n na

n

n nj jj

Q (f ) p x dx

p (x) l ( x )f ( x )

0

0

n

n nj jj

b

nj nja

nnj '

j n j

n

n ji

Q ( f ) a f x

a l x dx

W ( x ) l j 0,...,n

x x W x

W (x) x x

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:

0

0

j

n

b a h , x x jh,

n x a, x b.

FORMULA DEL TRAPEZIO

La f ormula di NEWTON - COTES a due punti in cui :

è detta f ormula del trapezio.

Ricaviamola per il generico intervallo e poi per

I l polinomio tale che :

0 1

0 1

1 1

x a, x b

-h, h a, b

x h, x h.

p (x) P

1

è dato da:

i i

1

h

1T h

p(x ) f (x ), i 0,1

h-x h+x p (x) f (-h) f (h)

2h 2h

Q (f ) I (p ) p x dx hf h hf h

1

0 0 1 1

0 1

Per ricavarlo per usiamo il

metodo dei coeffi cienti indeterminati.

imponiamo che il grado di precisione sia 1 e sia:

1

i iTi o

i a

a,b

Q ( f ) af ( x ) a f ( x ) af ( x )

f ( x ) ,x a a a xdx b a

2 2

0 1

0 1

2

2

1 b

i 0

1 b

i ai 0

T

b a af ( x ) aa ba xdx

b a da cui: a a

b-apertanto: Q (f ) f ( a) f ( b)

2

f(x)

a b

f(a)

f(b)

errore

Geometricamente:

1

1

Per ricavare l'errore ricordiamo che se

l'errore dell' interpolazione è :

1

1

n 1

( n ) n

n ii 0

( n )b b

n na a

f C [a,b]

fe(x) f (x)-p ( x ) W ( x ) W(x) ( x x )

n !

f per cui : e f ( x ) p ( x ) dx W ( x )dx

n !

ponen

1

1 1

2

( n )n 1

b

n n a

b ' 'T a

do : M max f ( x )

W ( x )si ha : e M dx

( n ) !

x a x b e f dx

Applichiamo ora il teorema del valor medio

sugli integrali per il quale , se

non cambia segno in

ponendo:

si ha:

b b

a a

' 'x

g, h C a,b e g(x)

a,b

n a,b : g(x)h x h x g x dx

g(x) x-a x b , h(x) f

e

Si può verifi care che il grado di precisione è 1.

' ' ' 'b 3

T a

f (n) f (n)x a x b b a

2 12

REGOLA DI SIMPSON

La formula di Newton-Cotes a 3 punti è detta regola diSimpson.

2 3

Poniamo:

ed imponiamo che :

1 2

0

23

0 1 2

2 h

2 i i hi 0

2

h

0 1 2 h

h

0 2 -h

h2 20 2 h

x h, x 0, x h

Q (f ) af (x ) f ( x ) dx

Per: f (x) 1, x, x

a a a dx h

-a h a h xdx

a h a h x dx h

Si può facilmente verificare che il grado diprecisione è 3 e ciò è sfruttato per determinare l’errore.

Infatti, poiché cambia segnoin [a,b] non si può procedere come prima.

0 1 2x-x x-x x-x

b f 2

ba f 4 (a) f

6a-b

(f)Q :b a, per e

f(h) f(0) 4 (-h) f 3h (f)Q

h34 a ,

3h a a

2

2

120

Si definisce invece il polinomio hermitiano

con le seguenti condizioni:

a cui può applicarsi il teorema del valore medio.

33 P (x)p

e poichè il grado di precisione è 3 :

a+b a+b a+b3 3 32 2 2

' 'a+b a+b2 2

I V2

3 0 1 2

p (a) f (a), p =f , p (b) f (b)

p f

ff x p x x-x x-x x-x

4!

Da cui:

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

5

90

I V b-a

S 2

-f ( x )e

Esempi:

Trapezio :

Simpson:

Generale:

2 nodi prec. =13 n 1, e h ,

3 nodi prec. =35n 2, e h ,

n disparin+1 nodi , prec.

n+1 pari2n 1 e h ,

Per aumentare la precisione si hanno 2 alternative:

i)Aumentare il numero di nodi in modo che sia integrale di un polinomio interpolante di alto

grado:Quadrature Gaussiane

ii)Si divide [a,b] in sottointervalli, in essi si usanoformule di bassa precisione, si sommano i risultati:Regole di Quadratura Composte.

Esaminiamo prima le quadrature composte

(f)Qn

Regole di Quadratura Composte

Suddividiamo [a,b] in n intervallini:

1-n

0j

x

xn1j

j

dx f(x) (f)Q

e usiamo in la regola del trapezio.

Sia :

Sommando si ha:

j 1

j

j j 1

j

x j 3j j 1x

n 1 nj 3

n j 0 nj 1 j 0

[x ,x ]

b-a h , x a jh, j 0,...., n

nf (n )h

f (x)dx f (x ) f (x ) h2 12

f (η )hT (f ) h f (x ) f (x ) f (x ) h

2 12

Per semplificare l’espressione dell’errore usiamo il

lemma:

1

0Sia

tutte dello stesso segno

n

j j

j

j

n 1 n-1

j j jj 0 i 0

g(x) C a,b e a

a

x a,b , j 0,...n-1 η a,b :

a g(x ) g a

2

12

12

I dentifi cando con con si ha:

e indicando con la sommatoria dimezzata agli estremi si ha:

3 3

3

j

n 2h h b at 12 12j

nh ( b a)

n jj 0

h f (η ) g(x) e a -

12

e -ff n f h

Σ

T (f ) h ' 'f (x ) f

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

ba,Cf I(f)fQ lim

Nn Ka : 0K Se

n. da dipendenti ainterpolat

quadratura della nodi i e pesi i sono x,a dove

xfa(f)Q ,ba,Cf Sia

nn

n

0j

(n)j

(n)j

(n)j

n

0j

(n)j

(n)jn

Dim.:

Per Weierstrass

Poichè la quadratura è interpolatoria:

N N

N

n N N

ε 0 q (x) P (N f ) :

f -q ε

Q (q ) I (q ) n N

0 0

Scegliendo si ha:

n N n N n N n N n

N N

n n(n) (n) (n) (n)

n N n j N j j N jj j

n

n N

I (f ) - Q (f ) I (f ) I (q ) Q (q ) Q (f ) I (f ) I (q ) Q (q ) Q (f )

I (f ) I (q ) f q (b a)

Q (q ) Q (f ) a q (x ) f (x ) f q a Kε

I (f ) Q (f ) k _

b a ε ε

Si può provare che è vero il viceversa

Se gli sono tutti la convergenza è garantita.

I nf atti , poichè il polinomio è integrato esattamente

si ha :

=

Quindi se i pesi sono tutti positivi :

(n)j

0

nb (n)n ja

j 0

(n)j

j 0

a 0

p (x) 1

0 I 1 dx Q ( 1 ) a

a

n n b(n)

j aj 0

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.

12

Si ha :

Perchè integri esattamente

si deve imporre

Per da cui :

n

n n jj 0

b an n

1 1

a+b1 2

Q (f ) a f (x )

f (x) 1

a

n 1: a b a x (b a)

Q (f ) (b a) f

MIDPOINT RULE

Metodo Midpoint

Integrazione esatta di un’approsimazione lineare di

Taylor dell’integranda.

Approssimazione lineare di Taylor ad f(x) in

b

adxxffI

2ba

c

1

1

poichè : 0

'

b

MP a

b '

a

p x f c x c f c

Q f p x dx b a f c

x c f c dx

Ricaviamo l’errore

che è la metà dell’errore del metodo dei trapezi.

212

' ' 'f x f c x c f c x c f

321 1

2 24

b ' ' ' 'MP a

E f x c f b a 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:

 

 

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.

n

0iiin )f(xa(f)Q

Per n =0 e [a,b] =[-1,1] :

 

 

imponendo che E0(f) = 0 per f(x) =1, x si ha:

 

  

1

010 01

01

1 22 0

0 2 0

dx aa ,x

xdx Q ( f ) f ( )

1

1

0 0 0

0 0

I ( f ) f ( x )dx

Q ( f ) a f ( x )

I ( f ) Q ( f ) E ( )

che per [a,b] generico dà:

 

che è la regola del punto di mezzoregola del punto di mezzo,, che quindi è di tipo Gaussiano.

 Per n=1:  

 

0 2a bQ (f ) ( b a)f ( )

Notiamo che tale formula ha grado di precisione 3 e usa 2 punti mentre la regola di Simpsonregola di Simpson per avere la stessa precisione usa 3 punti.

Quindi, in generale, si deve risolvere il sistema non lineare:

 

nelle 2n+2 incognite a0, …, an, x0, …, xn

Però, nell’ambito delle formule di Quadratura formule di Quadratura InterporlatorieInterporlatorie 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.

dxxxan

0i

b

a

rrii

1n2,...,0r

A tale scopo si ha: TEOREMATEOREMA

Se è una Formula di Quadratura di tipo

INTERPOLATORIO, ovvero:

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+1Tn+1 insieme dei polinomi ortogonali su [a,b], allora il grado di precisione della formula è 2n+1.

n

0iii )x(faQ (f)n

b

n naQ (f ) p ( x )dx

Dimostrazione

Sia:

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

 

b

a

n

0inii )f(E)x(fadx)x(f

pertanto:  

essendo la formula di tipo interpolatorio, essa ha almeno precisione n.

e poiché q(x)pn+1(x):

avendo imposto f P2n+1 la formula ha precisione 2n+1 Mostriamo ora che le formule Gaussiane hanno i pesi positivi.

n

0inii

b

a

b

a 1n

b

a)f(E)x(radx)x(rdx)x(q)x(pdx)x(f

b

a

n

0iii )x(radx)x(r

b

a n1n 0)f(E0dx)x(q)x(p

Se:

è Gaussiana, ha precisione 2n+1 e come f(x) prendiamo il quadrato dei polinomi di Lagrange:

, 0kn

(lk(x))2P2n e poiché: lk(xi)=ik si ha:

 

0kn, c.v.d.

n

0iiin )x(fa)f(Q

0

i

ik

nx xx xk

i ,i k

l ( x )

b

a k2

k adx(x))(l0

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 con

dove lnj sono i polinomi di Lagrange di grado n. Se

l’integrale da calcolare è del tipo:

2b

j jaa ln ( x )dx nj0

b

adx)x()x(f)f(I

Se in [-1,1] con p<1, q<1 i

polinomi sono quelli di JACOBI.

Se invece ovvero i polinomi

sono quelli di CHEBICHEV. Con tali polinomi i coefficienti

sono uniformi e per n nodi sono dati da: cioè:

dove (x) è una funzione peso tale che:

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).

b

adxx 0)(

qp x1x1)x(

21

qp )x1(

1)x(

2

n

Se i polinomi sono quelli di

LEGENDRE.

1)( x

1

021

n

ii

b

a)x(f

ndx

x

)x(f

Metodi di EstrapolazioneServono per prendere informazioni da poche Approssimazioni e usarle sia per stimare l’errore che per avere un’approssimazione migliore.Supponiamo che si abbia:

P=2 Trapezi, MidpointP=4 Simpsonnon valida per GaussServono per stimare p , l’errore, e migliorare l’approssimazione

pnn cfIfI

Estrapolazione di Richardson

Sia Q(h) una formula con accuratezza p ovvero:

Dove è un infinitesimo di ordine superiore a p usando un passo qh si ha:

Moltiplicando la (1°) per qp e sottraendo la (2°) si ottiene:

Con tale procedimento è possibile ottenere da una formula di basso ordine di accuratezza una formula di accuratezza maggiore.

( ) ( ) ( )p pI f Q h Ch o h ( )po h

( ) ( ) ( ) ( )p pI f Q qh C qh o h (2°)

(1°)

( ) ( )( ) ( )

1

pp

p

q Q h Q qhI f o h

q

Che per tanto ha un ordine più elevato. Se la formula di

partenza ammette uno sviluppo dell’errore del tipo:

Si ha:

1 21 2( ) ( ) .... ...kpp p

kI f Q h c h c h c h

1

( ) ( )( )

1

k

k

pk k

k p

q Q h Q qhQ h

q

Con

l’errore

1kphO

Stima di pSupponiamo di avere n, 2n, 4n punti e applichiamo la :

Consideriamo:

p

p

p

)n(n

)n(n

nn

cfIfI

cfIfI

cfIfI

44

22

nn

nnn II

IIr

42

24

ppp

p

pp

pp

)n(n

nn

nn

nnn

nn

nn

cc

cc

IIII

IIIIr

pp

pp

22412

24

2

42

2

42

24

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.

lg2

r lgp

r

4n

4n

p2

Per stimare l’errore si ha:

Inoltre

12

2

22

2

22

pnn

p

np

np

nn

III

IIccII pp

12

2 22

p

nnp

n

IIR

122

222

pnn

nnn

IIIRE

Integrazione di Romberg

L’integrale è:

Allora per due valori h1e h2 si ha:

e poiché:

Tale metodo si ottiene applicando l’estrapolazione di Richardson al metodo dei trapezi. In tale metodo l’errore è:

si ha:

2hhE

2

2

21

2

1

h

h

hE

hE

22

21

2 h

hhEh1E

)h(E)h(T)f(I

)h(E)h(T)h(E)h(T)f(I 2211

Quindi:

Da cui si ricava

che:

Sostituendo

:

2

11 2 2 22

2

hI ( h ) E ( h ) I ( h ) E ( h )

h

1 22 2

1

2

1

I ( h ) I ( h )E ( h )

hh

1 22 2 2 2

1

2

1

I ( h ) I ( h )I ( f ) I ( h ) E ( h ) I ( h )

hh

Siano

:

12 hI31

hI34

I

21

2

hh

2 1

2 2 1

4 14 1 3 3

I ( h ) I ( h )I ( f ) I ( h ) I h I ( h )

21

2

hh

Quindi indicando con il metodo dei trapezi si ha:

Trapezi

e l’errore :

0nT

0 0

1 22

43

n nn

T TT

14

41

21

12

j

jn

jn

jj

n

TTT

0

02

1

0 1

2 Kk k

( )n

( ) (1)n 2n

( )4n

( ) ( ) ( k )2 nn 2 n

T

T T

T

T T T

142

2

k

)k(n

)k(n)k(

n

TTE