Scarica Integrazione numerica e più Dispense in PDF di Modelli E Metodi Numerici solo su Docsity!
Capitolo 9
Integrazione numerica
Le formule di quadratura, ossia le formule per il calcolo di integrali definiti, prendono il loro nome da un problema antico della matematica: la quadratura del cerchio consisteva infatti nella ricerca di un quadrato equivalente ad un cerchio assegnato, e quindi dell’area della figura. Queste formule sono basate sull’idea di determinare il valore numerico di un integrale definito con una combinazione lineare di un numero finito di valori assunti dalla funzione inte- granda, o, equivalentemente, approssimare la funzione integranda mediante una funzioni elementare e integrare analiticamente quest’ultima.
9.1 Formule di quadratura
Il problema che prendiamo in esame `e il calcolo dell’integrale definito
I(f ) =
∫ (^) b
a
f (x) dx. (9.1)
Definiamo formula di quadratura un’espressione del tipo
In(f ) =
∑^ n
j=
αj f (xj ), (9.2)
ossia una combinazione lineare dei valori assunti dalla funzione integranda su n + 1 punti. I coefficienti αj vengono detti pesi e i punti xj nodi della formula. Fissato n, intendiamo approssimare I(f ) mediante In(f ). Diremo che la formula ha precisione algebrica pari a r, dove r `e un intero non negativo, se In(p) = I(p), per ogni p ∈ Πr,
9–
9–2 CAPITOLO 9. INTEGRAZIONE NUMERICA
cioe se la formula integra esattamente un polinomio di grado r. Una conseguenza immediata del fatto che una formula ha almeno preci- sione 0 (In(1) = I(1))e che debba risultare
∑^ n
j=
αj = b − a.
Molte formule di quadratura vengono date in particolari intervalli di ri- ferimento, come [− 1 , 1] o [0, 1]. Cio none restrittivo, in quanto ci si pu`o sempre ricondurre a tali intervalli mediante un cambio di variabile
Esempio 9.1 Applicando il cambio di variabile
x =
b − a 2 t^ +^
a + b 2 (9.3)
si ottiene ∫ (^) b
a
f (x) dx =
b − a 2
∫ (^1)
− 1
f
( b − a 2 t^ +^
a + b 2
) dt =
b − a 2
∫ (^1)
− 1
f^ ˜ (t) dt.
Un modo per costruire una formula di quadratura e il metodo dei coef- ficienti indeterminati. Esso consiste nel fissare arbitrariamente i nodi e determinare i pesi in modo che la formula abbia precisione algebrica n. Per la linearita dell’integrale `e sufficiente imporre questa condizione sulla base canonica In(xi) = I(xi), i = 0, 1 ,... , n,
ottenendo (^) n ∑
j=
xij αj =
bi+1^ − ai+ i + 1
, i = 0, 1 ,... , n,
ossia il sistema lineare XT^ α = b
la cui soluzione fornisce i pesi della formula. Dal momento che la matrice di Vandermonde X, di componenti Xji = xij (si veda la (8.5)), e non singolare quando i nodi sono distinti, resta dimostrato che una formula costruita su n + 1 nodi puo sempre avere precisione algebrica almeno n. Il metodo dei coefficienti indeterminati non e pero conveniente dal punto di vista numerico per quanto riguarda complessita e stabilita. Una meto- dologia alternativa conduce alle cosiddette formule di quadratura inter- polatorie. Una formula di questo tipo si ottiene sostituendo la funzione da
9–4 CAPITOLO 9. INTEGRAZIONE NUMERICA
In queste formule i pesi αj hanno la propriet`a di dipendere solo da n, oltre che dal passo h. Infatti, applicando il cambio di variabile x = x 0 + ht, nel caso delle formule chiuse si ricava
αj =
∫ (^) b
a
Lj (x) dx = h
∫ (^) n
0
Lj (x 0 + ht) dt = h
∫ (^) n
0
ϕj (t) dt,
dove
ϕj (t) = Lj (x 0 + ht) =
∏^ n
k k=0 6 =j
t − k j − k
La formula di quadratura (9.2) diventa allora
In(f ) = h
∑^ n
j=
wj f (xj ), con wj =
∫ (^) n
0
ϕj (t) dt.
Nel caso delle formule aperte, lo stesso cambio di variabile porta a
αj =
∫ (^) b
a
Lj (x) dx = h
∫ (^) n+
− 1
Lj (x 0 + ht) dt = h
∫ (^) n+
− 1
ϕj (t) dt.
Diamo ora alcuni esempi di formule di Newton-Cotes.
- Formula dei trapezi: n = 1, formula chiusa,
ϕ 0 (t) = 1 − t, w 0 =
0
(1 − t) dt =
ϕ 1 (t) = t, w 1 =
0
t dt =
I 1 (f ) =
h 2
[f (a) + f (b)] =
b − a 2
[f (a) + f (b)].
- Formula del punto medio o dei rettangoli: n = 0, formula aperta,
ϕ 0 (t) = 1, w 0 =
− 1
dt = 2,
I 0 (f ) = 2hf
a + b 2
= (b − a)f
a + b 2
9.2. FORMULE DI NEWTON-COTES 9–
- Formula di Simpson: n = 2, formula chiusa,
ϕ 0 (t) =
(t − 1)(t − 2), w 0 =
0
(t − 1)(t − 2) dt =
ϕ 1 (t) = t(2 − t), w 1 =
0
t(2 − t) dt =
ϕ 2 (t) =
t(t − 1), w 2 =
0
t(t − 1) dt =
I 2 (f ) =
h 3
[
f (a) + 4f
a + b 2
]
b − a 6
[
f (a) + 4f
a + b 2
]
Esercizio 9.1 Dal punto di vista geometrico, quali aree vengono valutate da queste formule?
Osservazione 9.1 Le formule di Newton-Cotes hanno sempre pesi simme- trici (il primo uguale all’ultimo, il secondo al penultimo, etc.).
Le formule come quelle appena trovate vengono dette elementari. Enaturale aspettarsi che tali formule non siano in grado di raggiungere preci- sioni elevate se non in presenza di una funzione integranda molto semplice. D’altra parte, aumentare il valore di n potrebbe non portare a buoni risul- tati. Infatti, sappiamo che l’interpolazione su un numero elevato di nodi equispaziati non garantisce la convergenza dell’errore neanche in presenza di una funzione infinitamente derivabile. Per aumentare la precisione si ricorre quindi ad una strategia che sfrutta la proprieta di additivit`a dell’integrale e che conduce alle formule composte. Consideriamo su [a, b] una discretizzazione di n + 1 punti equispaziati xi = xo + ih, i = 0, 1 ,... , n, h = b−na. Si ha ∫ (^) b
a
f (x) dx =
∑^ n
i=
∫ (^) xi
xi− 1
f (x) dx.
Una formula composta e basata sul calcolo di ciascun integrale con una for- mula elementare di ordine basso e la successiva somma delle approssimazioni ottenute. In questo modo l’intervallo di integrazione viene scomposto in un grande numero di intervalli piu piccoli su cui si auspica che la funzione in- tegranda sia sufficientemente semplice da essere bene approssimata da un polinomio interpolante di grado basso. Costruiamo ora le formule composte corrispondenti alle formule elemen- tari gi`a determinate.
9.2. FORMULE DI NEWTON-COTES 9–
Consideriamo ora una formula In con n pari e un polinomio pn+1 ∈ Πn+1. Mediante il cambio di variabile (9.3) si ottiene
I(pn+1) =
∫ (^) b
a
pn+1(x) dx =
− 1
qn+1(t) dt =
− 1
(γtn+1^ + qn(t)) dt,
dove il polinomio qn+1(t) ingloba il fattore differenziale b− 2 a e γ e il suo coef- ficiente di grado massimo. Poich´e tn+1^e una funzione dispari e la formula ha precisione almeno n,
I(pn+1) = γ
− 1
tn+1^ dt +
− 1
qn(t) dt =
− 1
qn(t) dt = In(qn),
il che mostra che la formula `e esatta se applicata ad un polinomio di grado n + 1. Di conseguenza, mentre la formula dei trapezi (n = 1) ha precisione algebrica 1, le formule del punto medio (n = 0) e di Simpson (n = 2) hanno precisione 1 e 3, rispettivamente.
9.2.2 Errore delle formule di Newton-Cotes
Il Teorema 8.5 dimostra che per ogni x esiste ξx tale che l’errore di interpo- lazione `e
f (x) − pn(x) =
f (n+1)(ξx) (n + 1)!
ωn(x).
Applicando questo risultato alle formule di quadratura, otteniamo
En(f ) = I(f ) − In(f ) =
∫ (^) b
a
(f (x) − pn(x)) dx =
∫ (^) b
a
f (n+1)(ξx) (n + 1)!
ωn(x) dx.
Se Mn+1 = maxx∈[a,b] |f (n+1)(x)|, allora
|En(f )| ≤
Mn+ (n + 1)!
∫ (^) b
a
|ωn(x)| dx.
Dal momento che l’errore di interpolazione su nodi equispaziati non tende necessariamente a zero quando n → ∞, la stessa affermazione vale anche per l’errore delle formule di Newton-Cotes elementari. Possiamo calcolare facilmente l’errore esatto per la formula dei trapezi, ricordando che vale il
9–8 CAPITOLO 9. INTEGRAZIONE NUMERICA
Teorema 9.1 (Secondo teorema della media) Se f (x) e una funzione continua e g(x)e integrabile e non cambia di segno in [a, b], allora esiste ξ ∈ [a, b] tale che ∫ (^) b
a
f (x)g(x) dx = f (ξ)
∫ (^) b
a
g(x) dx.
Infatti, poich´e la funzione ω 1 (x) = (x − a)(x − b) non cambia segno su [a, b], esiste η ∈ [a, b] tale che
E 1 (f ) =
∫ (^) b
a
f ′′(ξx)ω 1 (x) dx =
f ′′(η)
∫ (^) b
a
ω 1 (x) dx.
Inoltre, applicando il cambio di variabile x = a + ht, troviamo ∫ (^) b
a
ω 1 (x) dx = h
∫ (^) b
a
ht · h(t − 1) dt = −
h^3 6
da cui si ottiene l’espressione dell’errore
E 1 (f ) = −
f ′′(η) 12
(b − a)^3.
Un risultato analogo puo essere ottenuto anche per le altre formule, ma poich´e in questo caso la funzione ωn(x) cambia segno su [a, b],e necessario ricorrere ad una tecnica piu fine che fa uso del Teorema di Peano. E semplice estendere alle formule composte l’errore trovato per le formule elementari. Per la formula composta dei trapezi, dato che su ogni intervallino della decomposizione l’errore `e
E 1 (i )(f ) = −
f ′′(ηi) 12
h^3 , ηi ∈ [xi− 1 , xi],
poich`e abbiamo supposto f ∈ C^2 [a, b], esiste un η ∈ [a, b] tale che
E 1 (c )(f ) =
∑^ n
i=
f ′′(ηi) 12
h^3
h^2 12
(b − a)
n
∑^ n
i=
f ′′(ηi)
h^2 12
(b − a)f ′′(η).
Una conseguenza immediata, essendo h = b−n a, e che l’errore tende a zero per n → ∞. Nella Tabella 9.1 i risultati appena trovati per la formula dei trapezi sono confrontati con quelli validi per le formule dei rettangoli e di Simspon. Si noti come l’ordine di convergenza per le formule dei rettangoli e dei trapezi sia lo stesso (entrambe hanno precisione algebrica 1), ma che l’errore relativo alla primae la met`a rispetto alla seconda.
9–10 CAPITOLO 9. INTEGRAZIONE NUMERICA
Esempio 9.2 Costruire una formula su 2 nodi per approssimare l’integrale ∫ (^1)
0
f √ (x) x
dx,
con f funzione regolare. In questo caso utilizziamo w(x) = √^1 x come funzione peso. La presenza di una singolarit`a sul primo estremo dell’intervallo di integrazione ci obbliga ad utilizzare una formula aperta del tipo I 1 (f ) = α 0 f (x 0 ) + α 1 f (x 1 ). Fissiamo x 0 = 14 e x 1 = 34 e applichiamo il metodo dei coefficienti indeterminati
α 0 + α 1 =
∫ (^1)
0
w(x) · 1 dx = 2 1 4 α^0 +^
3 4 α^1 =
∫ (^1)
0
w(x) · x dx =
2 3 ottenendo α 0 = 53 e α 1 = 13. La formula ha precisione algebrica 1, cioe integra esattamente funzioni del tipo p (^1) √(x) x , per ogni polinomio di primo grado^ p^1. Spesso si ricorre a funzioni peso anche quando l’intervallo di integrazionee infinito, come mostra il seguente esempio.
Esempio 9.3 Applicando il cambio di variabile x = (^1) t , si ottiene ∫ (^) ∞
1
f (x) dx =
∫ (^1)
0
1 t^2 f
( 1 t^2
) dt.
Poich´e f (x) deve tendere a zero per x → ∞ con ordine superiore a 1, f
( (^1) t^2
) tende a zero per t → 0 con lo stesso ordine. Di conseguenza la nuova funzione integranda tende all’infinito con ordine inferiore a 1 ed e possibile trovare una funzione peso su [0, 1] che incorpori questa singolarita. Le funzioni peso possono essere impiegate anche per l’integrazione nu- merica di funzioni rapidamente oscillanti, senza bisogno di ricorrere ad un elevato numero di nodi. Infatti in molti casi tali funzioni possono essere fattorizzate nella forma w(x) · f (x), con w(x) continua e rapidamente oscillante e f (x) regolare. Ponendo, allora, w 0 = max[a,b] |w(x)|, si puo scrivere w(x)f (x) = (w(x) + w 0 )f (x) − w 0 f (x). La funzione w 0 f (x)e regolare e puo essere integrata senza problemi, mentre w(x) + w 0e continua e non negativa e pu`o essere assunta come funzione peso per ottenere una particolare formula di quadratura pesata.
9.4. FORMULE A PRECISIONE OTTIMALE 9–
9.4 Formule a precisione ottimale
Abbiamo visto che l’approccio per la costruzione di formule di quadratura interpolatorie consiste nel fissare n + 1 nodi in modo arbitrario e determinare di conseguenza i pesi. Questo consente di raggiungere la precisione algebrica n, o al piu n + 1 nel caso di nodi equispaziati e n pari. Si potrebbe, pero, pensare di lasciare liberi sia i pesi che i nodi e costruire, con un procedimento analogo al metodo dei coefficienti indeterminati, una formula di precisione ottimale (presumibilmente pari a 2n + 1, dato che le grandezze da determinare sono 2n + 2). Questo procedimento porterebbe ad un sistema di equazioni non lineari di difficile soluzione per n moderatamente elevato.
Esempio 9.4 Per costruire una formula ottimale per n = 1 della forma
I 1 (f ) = α 0 f (x 0 ) + α 1 f (x 1 ).
dovremmo imporre
α 0 xk 0 + α 1 xk 1 =
∫ (^) b
a
xk^ dx, k = 0, 1 , 2 , 3.
Il problema di determinare una formula a precisione ottimale e stato in realta elegantemente risolto con un approccio differente, basato sull’uso dei polinomi ortogonali, che conduce alle formule di quadratura Gaussiane. Per fare una trattazione pi`u generale, consideremo il calcolo dell’integrale
I(f ) =
∫ (^) b
a
w(x)f (x) dx (9.4)
per una qualsiasi funzione peso w(x). Nell’eventualita del calcolo di un integrale non pesato, bastera porre w(x) = 1.
9.4.1 Polinomi ortogonali
Scelta una funzione peso w(x), lo spazio L^2 w[a, b] delle funzioni per cui
∫ (^) b
a
w(x)(f (x))^2 dx < ∞
col prodotto scalare pesato
〈f, g〉 =
∫ (^) b
a
w(x)f (x)g(x) dx (9.5)
9.4. FORMULE A PRECISIONE OTTIMALE 9–
dato che 〈xpk, pk−j 〉 = 〈pk, xpk−j 〉 = 0 per ogni j ≥ 2, per l’ortogonalita dei polinomi gia costruiti. Ponendo αk = ck,k e β k^2 = ck,k− 1 ,
resta dimostrata la validit`a della formula (9.6). Tale formula risulta valida anche per k = 0, ponendo α 0 = (^) 〈〈px,p 0 ,p^00 〉〉 , β^20 = 0 e introducendo il polinomio fittizio p− 1 (x) = 0. Inoltre, ricavando xpk− 1 dalla formula ricorsiva di pk, si ottiene
〈xpk, pk− 1 〉 = 〈pk, xpk− 1 〉 = 〈pk, pk〉,
da cui risulta l’espressione alternativa data per β^2 k. Il fatto che questa quan- tita sia il rapporto tra due norme assicura la sua positivita e giustifica la sua rappresentazione come il quadrato di un numero positivo.
Utilizzando la definizione ricorsiva dei polinomi ortogonali e possibile dimostrarne alcune importanti proprieta.
Teorema 9.3 Il polinomio ortogonale pn+1 `e ortogonale allo spazio Πn, ossia
〈p, pn+1〉 = 0 per ogni p ∈ Πn.
Teorema 9.4 Le radici {xi}n i=1+1 di pn+1(x) sono semplici, reali e contenute nell’intervallo aperto (a, b).
Teorema 9.5 Se {xi}n i=1+1 sono gli zeri di pn+1(x) e {wi}ni=1 sono gli zeri di pn(x), allora
a < x 1 < w 1 < x 2 < w 2 < · · · < xn < wn < xn+1 < b.
Teorema 9.6 Fissata una qualsiasi discretizzazione {t 0 , t 1 ,... , tn}, la ma- trice (n + 1) × (n + 1)
A =
p 0 (t 0 ) · · · pn(t 0 ) .. .
p 0 (tn) · · · pn(tn)
`e non singolare se e solo se ti 6 = tj per i 6 = j.
9–14 CAPITOLO 9. INTEGRAZIONE NUMERICA
Dimostrazione. Supponiamo che i punti della discretizzazione siano di- stinti. Se A fosse singolare esisterebbe un vettore non nullo c ∈ Rn+1^ tale che Ac = 0. Allora il polinomio di grado n
q(x) =
∑^ n
j=
cj pj (x)
avrebbe n + 1 radici e, quindi, sarebbe identicamente nullo. Sia il massimo indice per cui c 6 = 0. Allora possiamo scrivere
p`(x) = −
c`
∑`−^1
j=
cj pj (x),
ma questo e assurdo perch´e i polinomi ai due lati dell’uguaglianza hanno grado e − 1, rispettivamente. La dimostrazione dell’implicazione inversae banale.
L’ultimo teorema ci dice che il determinante di Haar associato ad una base di polinomi ortogonali `e non nullo, e quindi che tali polinomi verificano la condizione di unisolvenza e costituiscono un sistema di Chebychev per qualsiasi discretizzazione a nodi distinti (§8.1).
9.4.2 Formule Gaussiane
Una formula di quadratura Gaussiana per l’integrale (9.4) e una formula di quadratura avente come nodi gli zeri del polinomio pn+1 ortogonale rispetto al prodotto scalare (9.5). I pesi sono scelti in modo da renderla interpolatoria, garantendo cosı una precisione algebrica almeno n
αj =
∫ (^) b
a
Lj (x)w(x) dx =
∫ (^) b
a
pn+1(x) cj (x − xj )
w(x) dx,
con cj =
∏n k=0,k 6 =j (xj^ −^ xk). I pesi possono essere trovati anche per altra via. L’ortogonalit`a implica che 〈pi, p 0 〉 = 〈p 0 , p 0 〉δi 0 , i = 0,... , n,
e il fatto che la precisione algebrica sia n, che
〈pi, p 0 〉 =
∫ (^) b
a
pi(x)w(x) dx = In(pi) =
∑^ n
j=
αj pi(xj ).
9–16 CAPITOLO 9. INTEGRAZIONE NUMERICA
considerando le propriet`a richieste ad una funzione peso e la precisione algebrica 2n + 1, si ottiene
∫ (^) b
a
w(x)pi(x) dx =
∑^ n
j=
αj pi(xj ) = αi
∏n
h h=0 6 =i
(xi − xh)^2 ,
da cui segue αi > 0, i = 0,... , n.
- Supponiamo ora di avere a disposizione una qualsiasi formula a preci- sione 2n + 1. Applicandola ai polinomi ortogonali pk(x), k = 0,... , n, troviamo ∑^ n
j=
αj pk(xj ) =
∫ (^) b
a
w(x)pk(x) dx = 〈pk, p 0 〉 = ‖p 0 ‖^2 δk 0 ,
ossia il sistema (9.8), il che mostra che i pesi sono necessariamente quelli di una formula Gaussiana. Applichiamo ora la formula ai polinomi pk(x)pn+1(x), k = 0,... , n, di grado non superiore a 2n + 1. Essendo la formula esatta per tali polinomi, si ha ∑^ n
j=
αj pk(xj )pn+1(xj ) =
∫ (^) b
a
w(x)pk(x)pn+1(x) dx = 〈pk, pn+1〉 = 0.
Otteniamo quindi il sistema omogeneo AT^ v = 0 dove la matrice A, data da (9.7), `e non singolare e pertanto
v = (α 0 pn+1(x 0 ),... , αnpn+1(xn))T^ = 0.
Essendo tutti i pesi positivi, questo implica che i nodi siano zeri di pn+1(x), e quindi che la formula sia Gaussiana.
- Supponiamo, per assurdo, che la formula sia esatta per ogni polinomio di grado 2n + 2, e applichiamola al polinomio non negativo di grado 2 n + 2 p(x) =
∏^ n
i=
(x − xi)^2 ,
avente i nodi di quadratura come zeri doppi. Allora
∫ (^) b
a
w(x)p(x) dx =
∑^ n
j=
αj p(xj ) = 0,
il che `e assurdo.
9.4. FORMULE A PRECISIONE OTTIMALE 9–
Nella determinazione di una formula Gaussiana, il problema piu comples- so appare essere il calcolo dei nodi, ossia degli zeri del polinomio ortogonale pn+1. Tuttavia,e possibile dimostrare che pn+1(λ) `e il polinomio caratteristico della matrice tridiagonale simmetrica
Jn =
α 0 β 1 β 1 α 1 β 2 β 2 α 2
βn βn αn,
le cui componenti non nulle sono i coefficienti della formula ricorsiva (9.6) che definisce pn+1. E evidente, dunque, che gli zeri del polinomio dato sonogli autovalori di Jn. Essendo tali zeri reali, distinti e compresi nell’intervallo (a, b) (Teorema (9.4)), l’algoritmo QR (§6.4)e particolarmente efficace per il loro calcolo. Per quanto riguarda l’errore delle formule Gaussiane, vale il seguente teorema.
Teorema 9.8 Se f ∈ C^2 n+1[a, b], allora esiste ξ ∈ (a, b) tale che ∫ (^) b
a
w(x)f (x) dx −
∑^ n
j=
αj f (xj ) =
f (2n+1)(ξ) (2n + 1)!
〈pn+1, pn+1〉.
9.4.3 Alcune famiglie di polinomi ortogonali
- Polinomi di Legendre: w(x) = 1, [a, b] = [− 1 , 1],
P 0 (x) = 1, P 1 (x) = x,
Pn+1(x) = xPn(x) −
n^2 4 n^2 − 1
Pn− 1 (x), n ≥ 1.
- Polinomi di Chebychev: w(x) =
1 − x^2
, [a, b] = [− 1 , 1], { T 0 (x) = 1, T 1 (x) = x, Tn+1(x) = 2xTn(x) − Tn− 1 (x), n ≥ 1.
Il peso di Chebychev e in grado di assorbire in un integrale con singo- larita agli estremi avente ordine di infinito non superiore a 12. Inoltre, poich´e Tn+1(x) = cos(n + 1)θ, con x = cos θ,
