Integrazione numerica di funzioni con singolarit`a

UNIVERSITÀ DEGLI STUDI DELLA CALABRIA
Facoltà di Scienze Matematiche, Fisiche e Naturali
Corso di Laurea in Matematica
Integrazione numerica di
funzioni con singolarità
RELATORE
Dr. Francesco Dell’Accio
CANDIDATO
Contartese Fortunato
Matr. 80432
Anno Accademico 2004-05
i
Indice
0 Introduzione
1 Quadratura numerica classica
1.1 Introduzione . . . . . . . . . . . . . . . . . .
1.2 Formule di quadratura . . . . . . . . . . . .
1.3 Formule di quadrature di tipo interpolatorio
1.4 Le formule di Newton-Cotes . . . . . . . . .
1.5 Formule di Newton-Cotes composte . . . . .
1.6 Studio dell’errore . . . . . . . . . . . . . . .
1
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
3
3
4
5
8
11
13
2 Quadratura numerica per funzione integranda limitata
2.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Integrali di funzioni con discontinuità di prima specie . . .
2.3 Integrazione adattiva . . . . . . . . . . . . . . . . . . . . .
2.4 Formula di Simpson adattiva . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
20
20
21
22
23
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
3 Quadratura numerica per funzione integranda non limitata
integrazione non limitato
3.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Idee generali . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Metodi per il calcolo di integrali impropri . . . . . . . . . . .
3.4 Procedimento con limite . . . . . . . . . . . . . . . . . . . . .
3.5 Troncamento dell’intervallo . . . . . . . . . . . . . . . . . . .
3.6 Cambio di variabile . . . . . . . . . . . . . . . . . . . . . . . .
3.7 Integrazione di tipo interpolatorio . . . . . . . . . . . . . . .
3.8 Cambio di variabile su intervalli non limitati . . . . . . . . . .
3.9 Procedimento al limite per intervalli infiniti . . . . . . . . . .
3.10 Accelerazione della convergenza . . . . . . . . . . . . . . . . .
3.11 Integrazione Gaussiana su intervalli illimitati . . . . . . . . .
Bibliografia
o intervallo di
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
32
32
33
34
37
38
39
40
41
42
43
44
48
1
Capitolo 0
Introduzione
Per integrali singolari si intende sia la classe degli integrali di funzioni con
punti di discontinuità (della funzione o delle sue derivate) sia quella degli integrali
su domini illimitati. Le singolarità possono quindi presentarsi sotto varie forme e
non è ragionevole pensare che esista un unico modo generale e conveniente per trattarle. Piuttosto, vi possono essere alcune idee da adottare caso per caso. In questa
tesi analizzeremo alcune di queste idee, limitandoci al solo caso di integrali di funzioni univariate. Bisogna comunque sottolineare che è l’analisi del caso particolare
o una sua opportuna trasformazione a suggerire in generale la soluzione più idonea.
L’elaborato si sviluppa in tre capitoli. Nel primo sono richiamati alcuni concetti basilari di teoria della quadratura numerica cui si fa riferimento nei capitoli successivi.
Nel secondo capitolo viene analizzato il caso della funzione integranda limitata, con
un numero finito di discontinuità nel dominio di integrazione, e introdotto il metodo
adattivo per il calcolo numerico degli integrali. Nel terzo conclusivo capitolo è con-
2
siderato il caso della funzione integranda illimitata, e sono introdotte diverse idee per
un’opportuna quadratura numerica. Infine è anche trattato il caso dell’integrazione
numerica su domini illimitati.
3
Capitolo 1
Quadratura numerica classica
1.1
Introduzione
In questo primo capitolo richiamiamo metodi numerici classici per il calcolo
di integrali definiti del tipo
Z
b
I(f ) =
f (x)dx,
(1.1)
a
dove f è una funzione reale integrabile sull’intervallo [a, b] e sufficientemente regolare,
che saranno utili nel seguito della tesi. In particolare tratteremo il caso delle formule
di quadratura di tipo interpolatorio, cioè formule del tipo
Z
b
f (x)dx =
a
n
X
Ai f (xi ) + Rn (f )
i=1
che risultano esatte nello spazio dei polinomi in x di grado ≤ n. Allo scopo di aumentarne la precisione, queste formule vengono combinate fra loro ripartendo l’intervallo
[a, b] in sottointervalli, sovente equispaziati, a formare formule di quadratura di tipo
composito, parimenti trattate in questo primo capitolo.
4
1.2
Formule di quadratura
Sia f una funzione reale integrabile sull’intervallo [a, b]. Una formula di
quadratura è un’espressione del tipo
Z
b
f (x)dx =
a
n
X
Ai f (xi ) + Rn (f )
(1.2)
i=1
dove gli Ai ∈ R sono detti coefficienti o pesi, e gli xi i = 1, .., n sono punti distinti
appartenenti all’intervallo [a, b] detti nodi. La quantità
Qn (f ) =
n
X
Ai f (xi )
i=1
è una approssimazione dell’integrale (1.1) con un errore dato da
Z
Rn (f ) =
b
f (x)dx −
a
n
X
Ai f (xi ).
(1.3)
i=1
Dipendentemente dai valori che assumono i pesi e i nodi, incognite nella formula (1.2),
si ottengono formule di quadratura più o meno esatte, come specificato nella seguente
definizione.
Definizione 1 Diremo che la formula di quadratura (1.2) ha grado di esattezza m se
l’errore di troncamento (1.3) è nullo per ogni polinomio p(x) di grado ≤ m.
E’ possibile calcolare i nodi xi e i pesi Ai in modo tale che la formula di
quadratura abbia grado di esattezza massimo, ossia pari a 2n − 1, ([1, p.272]). In
particolare, riferendosi all’intervallo [−1, 1], per i pesi Ai si ha la formula
Z
Ai =
1
−1
ω(x)
dx,
(x − xi )ω 0 (xi )
(1.4)
5
dove con ω(x) = (x − x0 )(x − x1 ) · · · (x − xn ) è stato denotato il polinomio nodale,
mentre i nodi risultano essere le radici del polinomio
Pn (x) =
1 dn 2
(x − 1)n ;
2n n! dxn
(1.5)
formule generali si ottengono a partire da quelle precedentemente date mediante una
semplice trasformazione lineare. Il polinomio nella (1.5) è detto polinomio ortogonale
di Legendre, di grado n. Nel caso in cui i nodi siano fissati a priori, possiamo scegliere
i pesi Ai in modo tale che la formula di quadratura abbia grado di esattezza massimo.
Questa tecnica conduce al concetto di formula di quadratura di tipo interpolatorio
che ci apprestiamo a trattare.
1.3
Formule di quadrature di tipo interpolatorio
Definizione 2 Una formula di quadratura su n + 1 nodi distinti, si dice di tipo inter-
polatorio, se ha grado di esattezza almeno n.
Proposizione 3 Una formula di quadratura ha grado di esattezza n se solo se per
l’errore di troncamento Rn associato risulta Rn (xk ) = 0 per k = 0, 1, ..., n.
Dimostrazione. Qui e in seguito denotiamo con P n lo spazio dei polinomi
in x di grado ≤ n. Sia p ∈ P n un generico polinomio:
p(x) = a0 + a1 x + ... + an xn .
Poichè Rn è un funzionale lineare risulta:
Rn (p(x)) = R(a0 + a1 x + ... + an xn )
= a0 Rn (1) + a1 Rn (x) + ... + an Rn
(1.6)
(xn ).
6
Poichè l’identità (1.6) vale per ogni scelta dei coefficienti a0 , a1 , ..., an , possiamo concludere con la tesi.
Il seguente teorema garantisce l’esistenza di formule di quadratura di tipo
interpolatorio.
Teorema 4 Dati n + 1 punti distinti xi i = 0, ..., n appartenenti all’intervallo [a, b],
esiste ed è unica la formula di quadratura di tipo interpolatorio avente come nodi i
punti dati.
Dimostrazione. Essendo i nodi n + 1, la formula di quadratura (1.2) è di tipo
interpolatorio se, e solo se, sono verificate le condizioni Rn (xk ) = 0 k = 0, ..., n. Inponendo che tali condizioni siano verificate otteniamo il seguente sistema di equazioni
lineari
n
X
Ai xki = mk k = 0, ..., n
(1.7)
i=0
con
mk =
bk+1 − ak+1
k+1
k = 0, ..., n.
Il sistema (1.7) può essere scritto nella forma matriciale
VA=M
(1.8)
7
dove A = [A0 , ..., An ]T , M = [m0 , ...., mn ]T e









V =







1
1
x0
x1
:
:
:
:
xn0 xn1
.. ..
1 


.. .. xn 



: : : 




: : : 


.. .. xnn
Essendo la matrice V non singolare in quanto matrice di Vandermonde sui
nodi distinti [1, p.194], il sistema (1.8) ammette un’unica soluzione che fornisce i pesi
della formula di quadratura (1.2).
Fissati gli n + 1 nodi xi i = 0, ..., n, tutti distinti e appartenenti all’intervallo
[a, b], supponiamo che la formula di quadratura (1.2) sia di tipo interpolatorio. Poichè
ogni polinomio p(x) di grado ≤ n può essere scritto nella formula di Lagrange
p(x) =
n
X
i=0
ω(x)
p(xi ),
(x − xi )ω 0 (xi )
risulta con facili calcoli
Rn (p(x)) =
Z bX
n
a i=0
n
X
ω(x)
p(x
)dx
−
Ai p(xi ) = 0.
i
(x − xi )ω 0 (xi )
i=0
La linearità dell’integrale ci consente di scrivere la precedente equazione nella forma
n
X
i=0
·Z
p(xi )
b
a
¸
ω(x)
dx − Ai = 0
(x − xi )ω 0 (xi )
da cui, per l’arbitrarieta del polinomio p(x)
Z
Ai =
a
b
ω(x)
dx
(x − xi )ω 0 (xi )
(1.9)
8
1.4
Le formule di Newton-Cotes
Le formule di quadratura di tipo interpolatorio con nodi ugualmente spaziati
in [a, b] sono dette formule di Newton-Cotes. Una proprietà interessante delle formule di Newton-Cotes è quella di avere i coefficienti di integrazione Ai i = 0, ..., n
dipendenti solo da n e da h =
b−a
. Questa proprietà si verifica facilmente calcolando
n
esplicitamente i coefficienti Ai a partire dalla formula (1.9). Le formule di quadratura
di Newton-Cotes si dicono chiuse se tra i nodi figurano gli estremi dell’intervallo di
integrazione, altrimenti si dicono aperte. Nel seguito limiteremo la nostra attenzione
al solo caso delle formule chiuse, di maggiore rilevanza in questo lavoro. Allo scopo
di determinare i coefficienti Ai i = 0, .., n ed il resto Rn (f ), nel caso delle formule
chiuse, riscriviamo l’espressione della generica formula di quadratura sugli n + 1 nodi
xi = a + ih, i = 0, ..., n nella forma
Z
b
f (x)dx = (b − a)
a
n
X
Bkn f (xk ) + Rn (f )
(1.10)
k=0
dove è stato posto
Bkn = (b − a)−1 Ak .
(1.11)
Come precedentemente mostrato i pesi Ak delle formule di quadratura di tipo interpolatorio si ricavano mediante la formula
Z
Ak =
a
b
ω(x)
dx, k = 0, ..., n
(x − xk )ω 0 (xi )
dove
ω(x) = (x − x0 )(x − x1 ) · · · (x − xn ).
(1.12)
9
Dalle uguaglianze precedenti risulta quindi
Z
Bkn
−1
= (b − a)
a
b
ω(x)
dx, k = 0, ..., n.
(x − xk )ω 0 (xk )
(1.13)
Introducendo la nuova variabile t, mediante la trasformazione x = x0 + th, si ha
x − xk = h(t − k) k = 0, ..., n,
ω(x) = hn+1 πn (t) con πn (t) = t(t − 1) · · · (t − n)
e
ω 0 (xk ) = ω 0 (x0 + kh) = (−1)n−k hn k!(n − k)!,
per cui, sostituendo nella (1.13) risulta
Bkn =
(−1)n−k
nk!(n − k)!
Z
n
t(t − 1) · · · (t − i + 1)(t − i − 1) · · · (t − n)dt.
(1.14)
0
Pertanto i coefficienti Bkn possono essere facilmente tabulati. Nella tabella che segue
n , come si può facilmente
sono riportati i valori Bkn per n = 1, ..., 10; poichè Bkn = Bn−k
verificare, ci limitiamo ad indicare solo i valori per k ≤
hni
2
+ 1.
10
n
B0n
B1n
B2n
B3n
1
1
2
2
1
6
4
6
3
1
8
3
8
4
7
90
32
90
12
90
5
19
288
75
288
50
288
6
41
840
216
840
27
840
272
840
7
751
17280
3577
17280
1323
17280
2989
17280
8
989
28350
5888
28350
9
2857
89600
15741
89600
10
16067
598752
106300
598752
−
928
28350
10496
28350
1080
89600
19344
89600
48525
598752
272400
598752
−
B4n
−
B5n
4540
28350
5778
89600
−
260550
598752
427368
598752
Dalla tabella precedente risulta che i coefficienti per n = 8 e n = 10 non sono tutti
positivi; in generale è possibile dimostrare [3, pag. 151] che i pesi sono tutti positivi
solo per n ≤7 e n =9. In conseguenza di ciò lo scarso utilizzo delle formule di NewtonCotes per n ≥ 8. Le formule di Newton Cotes più semplici si ottengono per n =1, 2 e
prendono il nome rispettivamente di formula del trapezio e di Simpson. L’espressione
della formula del trapezio è
Q1 (f ) =
b−a
[f (a) + f (b)] .
2
(1.15)
11
f (x )
x
b = x1
a = x0
Figura 1.1: Formula del Trapezio
L’errore di quadratura per questa formula, nel caso in cui f ∈ C 2 ([a, b]), è dato da
R1 (f ) = −
1
(b − a)3 f (2) (ξ), ξ ∈ (a, b) .
12
(1.16)
La formula di Simpson è
·
µ
¶
¸
b−a
a+b
Q2 (f ) =
f (a) + 4f
+ f (b)
6
2
(1.17)
con errore di quadratura, nel caso f ∈ C 4 ([a, b]), dato da
R2 (f ) = −
b−a
h5 (4)
f (ξ), h =
, ξ ∈ (a, b) .
90
2
(1.18)
Le espressioni dei resti delle formule del trapezio e dei Simpson saranno esplicitamente
ricavate nel paragrafo 1.6.
1.5
Formule di Newton-Cotes composte
La procedura generale consiste nel suddividere l’intervallo di integrazione
12
f (x)
a = x0
x1 =
a +b
2
b = x2
Figura 1.2: Formula di Simpson
[a, b] in N sottointervalli. Ponendo:
yi = a + iH i = 0, ..., N, H =
b−a
N
risulta, per la proprietà di additività dell’integrale,
Zb
N
−1
X
f (x)dx =
f (x)dx;
(1.19)
k=0 yi
a
ciascuno degli integrali
yZi+1
yZi+1
f (x)dx
(1.20)
yi
può essere approssimato con una delle formule di quadratura precedentemente discusse; se si usa la stessa per ogni i = 0, ..., N − 1 si ottiene una formula di quadratura
composta per il calcolo di
Rb
f (x)dx. Applicando alla (1.20) la formula di Newton-Cotes
a
chiusa su n + 1 nodi abbiamo
yZi+1
f (x)dx = H
yi
n
X
£ n
¤
Bk f (xk ) + Rni (f )
k=0
(1.21)
13
H
e sostituendo nella (1.19) risulta
n
·
¸
NP
−1
n
Rb
P
n
i
f (x)dx =
H
Bk f (xk ) + Rn (f )
con xk = yi + kh, k = 0, ..., n, h =
i=0
k=0
NP
−1 P
n
a
Bkn f (xk ) +
=H
=H
i=0 k=0
n
NP
−1
P
Bkn [
f (yi
i=0
k=0
NP
−1
i=0
Rni (f )
+ kh)] +
NP
−1
i=0
(1.22)
Rni (f ).
Inoltre, poichè yi = yi−1 +nh, i valori f (yi ), i = 1, ..., N −1 appaiono due volte nell’ultimo
membro di (1.22), per cui possiamo scrivere
Rb
a
(
NP
−1
B0n f (y0 ) + Bnn f (yN ) + (B0n + Bnn )
f (yi )
j=0
"
#)
n−1
N
P n P
+
Bk
f (yj−1 + kh)
+ RN (f )
f (x)dx = H
k=1
(1.23)
j=1
con
RN (f ) =
N
−1
X
Rni (f ).
(1.24)
i=0
La (1.22) prende il nome di formula di quadratura di Newton-Cotes composta. Nel
prossimo paragrafo, dopo aver esplicitimante ricavato espressioni per l’errore nelle
formule di quadratura di Newton-Cotes, mostreremo come da queste è possibile ottenere espressioni analitiche per l’errore (1.24) delle corrispondenti formule composte
(1.23).
1.6
Studio dell’errore
Scopo di questo paragrafo è lo studio dell’errore nei metodi di quadratura,
con l’obiettivo particolare di evidenziare l’ordine di convergenza. Ci poniamo quindi
nel seguente quadro generale: la funzione w(x) è una funzione peso positiva su (a, b)
(limitato o meno) e tale che xn w(x) ∈ L1 (a, b), ∀n ∈ N. Per funzioni f (x) tali che
14
w(x)f (x) ∈ L1 (a, b) consideriamo formule di quadratura del tipo:
Zb
w(x)f (x)dx ≈
n
X
Ai f (xi ), Ai ∈ R, xi ∈ (a, b).
(1.25)
i=0
a
In tale quadro rientrano sia le formule di quadratura elementari che quelle composte.
Per il seguito indicheremo con [α, β] un intervallo contenente a, b; l’errore della formula
di quadratura sarà indicato, al solito, con
Zb
R(f ) =
w(x)f (x)dx −
n
X
Ai f (xi ).
i=0
a
Come noto, R(f ) è un funzionale lineare e continuo su C 0 ([α, β]), rispetto alla norma
del massimo. Infine, indicando con u+ = max(u, 0), porremo, per t fissato e m intero
non negativo:
(x − t)m
+ =



 (x − t)m per x ≥ t,


 0
per x ≤ t.
Si ha il seguente
Teorema 5 Supponiamo che la formula di quadratura (1.25) sia esatta per polinomi
di grado ≤ N, N ≥ 0 e che f ∈ C N +1 ([α, β]). Allora
1
R(f ) =
N!
Zβ
KN (t)f (N +1) (t)dt,
(1.26)
α
dove la funzione nucleo
KN (t) = R((x − t)N
+ ),
corrispondente all’errore della formula di quadratura applicata alla funzione particolare x −→ (x − t)N
+ , è detta nucleo di Peano della formula (1.25).
15
Dimostrazione. Dall’ uguaglianza
Zx
f 0 (t)dt
f (x) = f (α) +
α
per integrazione per parti si ha
Zx
0
(t − α)f 00 (t)dt.
f (x) = f (α) + (x − α)f (x) −
α
Analogamente


Zx
(x − α)f 0 (x) = (x − α) f 0 (α) +
f 00 (t)dt
α
per cui, sostituendo nella precedente uguaglianza, otteniamo con facili calcoli
Zx
0
(x − t)f 00 (t)dt.
f (x) = f (α) + (x − α)f (α) +
α
Procedendo in questo modo, dopo N − 1 ulteriori integrazioni per parti, otteniamo
f (N ) (α)
1
f (x) = f (α) + f 0 (α) (x − α) + ... +
(x − α)N +
N!
N!
= TN [f, α](x) +
1
N!
Zβ
Zx
(x − t)N f (N +1) (t)dt
α
(N +1) (t)dt
(x − t)N
+f
α
Essendo TN [f, α](x) un polinomio di grado minore o uguale ad N , per la proprietà di
esattezza della formula di quadratura, risulta
 β

Z
1 
(N +1) (t)dt
R(f ) = R(TN [f, α]) +
R
(x − t)N
+f
N!
α


 b β
Zβ
Z Z
n
P
1  
(N +1) (t)dt
(N +1) (t)dt w(x)dx −
Ai (xi − t)N
(x − t)N
=
+f
+f
N!
i=0
a
α
α
da cui, dopo aver scambiato l’ordine di integrazione nell’integrale doppio ad ultimo
16
membro dell’uguaglianza precedente, otteniamo
R(f ) =
1
N!
=

Zβ Zb
α
1
N!

Zβ
(x − t)N
+ w(x)dx −
a
n
P
i=0
Ai (xi − t)N
+



f (N +1) (t)dt
KN (t)f (N +1) (t)dt.
α
Corollario 6 Se la funzione t → K N (t) non cambia segno su [α, β], allora esiste un
ξ ∈ (α, β) tale che:
R(f ) =
f (N +1) (ξ)
R(xN +1 )
(N + 1)!
Dimostrazione. Applicando alla (1.26) il primo teorema del valor medio per
gli integrali [4, p. 7] si ha:
R(f ) =f
(N +1)
1
(ξ)
N!
Zβ
KN (t)dt,
ξ ∈ (α, β) .
(1.27)
α
Prendendo, in particolare, f (x) = xN +1 si ha:
Zβ
Zβ
R(xN +1 ) = (N + 1)
KN (t)dt =
KN (t)dt =⇒
α
α
R(xN +1 )
.
(N + 1)
(1.28)
Da cui, sostituendo la (1.28) nella (1.27) otteniamo la tesi.
Applichiamo il corollario precedente per calcolare espressioni dell’errore di alcune formule di quadratura. Oltre alle già citate formule del trapezio e di Simpson, prendiamo
in esame la formula del rettangolo
Z
b
f (x)dx ≈ f (a)
a
e quella del punto medio
Z
µ
b
f (x)dx ≈ f
a
a+b
2
¶
.
17
Metodo
Errore
Metodo del rettangolo
f (ξ)
Metodo del punto medio
f (ξ)
Metodo del trapezio
−f (ξ)
Metodo di Simpson
−f
0
(b − a)2
2
00
00
(4)
(b − a)3
24
(b − a)3
12
(ξ)
[(b − a)2 /2]5
90
Nota un’espressione dell’errore per una data formula elementare, è del tutto
immediato lo studio dell’errore della corrispondente formula composta. A tale scopo
è d’ausilio il seguente
Lemma 7 Sia u ∈ C 0 ([a, b]) e siano dati in [a, b] punti xj j = 0, ..., s + 1 e in corrispon-
denza costanti δj , tutte dello stesso segno. Esiste allora almeno un punto ξ ∈ [a, b]
tale che
s
X
δj u(xj ) = u(ξ)
j=0
s
X
δj .
(1.29)
j=0
Dimostrazione. Siano um = min u(x) = u(xm ) e uM = max u(x) = u(xM )
x∈[a,b]
x∈[a,b]
con xm , xM ∈ [a, b]. Dato che i coefficienti δj per ipotesi devono essere tutti dello stesso
segno, procederemo nella dimostrazione supponendoli positivi. Per definizione di um ,
uM si ha:
s
s
s
X
X
X
um δj ≤
δj u(xj ) ≤ uM
δj .
j=0
Poniamo σs =
j=0
(1.30)
j=0
s
s
X
X
δj u(xj ) e consideriamo la funzione continua U (x) = u(x) δj . In
j=0
j=0
virtù della (1.30) risulta U (xm ) ≤ σs ≤ U (xM ). Per il teorema del valor medio esiste
almeno un punto ξ compreso tra a e b tale che U (ξ) =σ . Una dimostrazione del tutto
18
simile può eseguirsi nel caso in cui i coefficienti δj siano negativi.
Suddividiamo l’intervallo (a, b) negli intervalli [yi , yi+1 ], con yi = a + iH ,
H = (b − a)/N ; come constatato nel paragrafo precedente per l’errore nella formula
composita risulta
RN (f ) =
N
−1
X
Rni (f ),
i=0
dove con
Rni (f )
abbiamo indicato l’errore della formula di quadratura elementare
quando applicata al generico intervallo [yi , yi+1 ] i = 0, ..., N − 1. Consideriamo come
esemplificazione il caso della formula dei trapezi. Sappiamo che se la funzione f (x) è
sufficientemente regolare l’errore che si commette utilizzando tale formula sul generico
intervallo [xi , xi+1 ] è dato da
R1i (f ) =
h3 (2)
f (ξi ), ξ ∈ ]yi , yi+1 [ .
12
Si ha, pertanto, denotando con T (H) la formula del trapezio composita (o dei trapezi)
Z
T (H) −
b
f (x)dx =
a
N −1
N −1
H 3 X (2)
H2
1 X (2)
f (ξi ) =
(b − a)
f (ξi ).
12
12
N
i=0
i=0
Se f ∈ C 2 ([a, b]) per il lemma precedente, assumendo δi =
1
i=0, ..., N esiste almeno
N
un punto ξ nell’intervallo [a, b] tale che
f
(2)
N −1
1 X (2)
(ξ) =
f (ξi )
N
i=0
In definitiva, si ha quindi
Zb
f (x)dx − T (H) = −
a
b − a 2 (2)
H f (ξ), H = (b − a)/N.
12
Si procede in modo del tutto analogo per il calcolo degli errori nelle formule composite
ottenute a partire dalle formule di quadratura discusse nel paragrafo precedente; in
particolare nella tabella che segue riportiamo alcune di queste espressioni:
19
Metodo
Errore
Metodo dei rettangoli
f 0 (ξ)
Metodo dei punti medi
f 00 (ξ)
Metodo dei trapezi
−f 00 (ξ)
Metodo di Simpson composito
−f (4) (ξ)
(b − a)
h
2
(b − a) 2
h
24
(b − a) 2
h
12
(b − a) 4
h
180
20
Capitolo 2
Quadratura numerica per funzione
integranda limitata
2.1
Introduzione
In questo capitolo analizzeremo un metodo di integrazione numerica auto-
matica: un insieme, cioè, di algoritmi numerici che restituisce un’approssimazione
dell’integrale I(f ) =
Rb
a
f (x)dx, entro i limiti di una tolleranza ε precisata dall’utente
e tale da essere capace di modificare in modo automatico il passo di integrazione,
dipendentemente dal comportamento locale della funzione. Svilupperemo tale idea
utilizzando in particolare la formula di Simpson per la quale riporteremo anche una
implementazione dimostrativa.
21
2.2
Integrali di funzioni con discontinuità di prima specie
Se la funzione integranda è limitata, ma non sufficientemente regolare, in-
tendento con ciò che la funzione stessa o alcune sue prime derivate presentano salti
in punti la cui collocazione è nota a priori, è possibile decomporre il dominio di integrazione in maniera che tali punti risultino estremi degli intervalli di integrazione
parziale. Più precisamente, sia c un punto noto all’interno di [a, b] e sia f una funzione
continua e sufficientemente regolare in [a, c) e (c, b], con salto f (c+ )−f (c− ) finito; si ha
allora
Zb
I(f ) =
Zc
f (x)dx =
a
Zb
f (x)dx +
a
f (x)dx
c
Utilizzando separatamente su [a, c− ) e (c+ , b] una qualsiasi delle formule di integrazione numerica precedentemente considerate, si può approssimare correttamente
I(f ). Le stesse considerazioni valgono se f ammette un numero finito di discontinuità
all’interno di [a, b] note. Più difficile, ovviamente, si presenta la situazione quando i
punti di discontinuità non sono noti a priori. Se la funzione è data in forma analitica, per localizzare almeno approssimativamente tali punti può essere utile un esame
della funzione mediante le tecniche dell’analisi. Un prezioso supporto, se usato correttamente, è lo studio grafico. Quando la funzione non è data in forma esplicita e
comunque risulta troppo complicata da studiare, si può affidare al metodo di integrazione il compito di scoprire eventuali punti di discontinuità. Evidentemente tali
metodi, noti col nome di metodi adattivi, dovranno essere capaci di modificare in
modo automatico il passo di integrazione, dipendentemente dal comportamento lo-
22
cale della funzione. E infatti, quando si applica un metodo adattivo, la presenza di
punti di discontinuità della funzione è segnalata da un eccessiva riduzione del passo
di integrazione. Naturalmente il successo di questa idea presuppone che il metodo
adattivo sia efficiente.
2.3
Integrazione adattiva
L’obiettivo di un integratore adattivo è primariamente quello di fornire una
Zb
f (x)dx entro i limiti di una prefissata tolleranza ε. La
approssimazione di I(f ) =
a
scelta della suddivisione in parti uguali dell’intervallo di integrazione in una formula
di quadratura composta, può non essere appropriata quando si integra una funzione
con comportamento altamente variato nell’intervallo di integrazione, cioè con grandi
variazioni su alcune parti e piccole variazioni su altre. Il significato di variazione è
precisato dal comportamento delle derivate della funzione. Fissiamo l’attenzione sul
generico sottointervallo [α, β] ⊂ [a, b]. L’analisi delle stime dell’errore per le formule
di Newton-Cotes, alcune delle quali sono state riportate nel capitolo precedente, evidenzia che bisognerebbe valutare le derivate di f sino ad un certo ordine per potere
scegliere un passo h di integrazione tale da garantire un’accuratezza prefissata, diciamo ²
(β − α)
. Tale procedimento, che richiederebbe uno sforzo computazionale ec(b − a)
cessivo, risultando cosı̀ impraticabile nelle modalità appena descritte, viene realizzato
in un integratore automatico adattivo.
23
2.4
Formula di Simpson adattiva
Per fissare le idee ci riferiamo alla formula di Simpson (1.17), sebbene il
metodo che ci apprestiamo a considerare possa essere applicato anche ad altre formule
di quadratura. Poniamo
Zβ
If (α, β) =
α
·
µ
¶
¸
h
α+β
β−α
f (x)dx, Sf (α, β) =
f (α) + 4f
+ f (β) , h =
.
3
2
2
L’espressione dell’errore nella formula di Simpson (1.18), nell’ipotesi che f sia sufficientemente regolare in [α, β] , è
If (α, β) − Sf (α, β) = −
h5 (4)
f (ξ),
90
(2.1)
essendo ξ un punto interno all’intervallo [α, β]. Allo scopo di stimare l’errore If (α, β)−
Sf (α, β) senza calcolare esplicitamente f (4) (x) utilizziamo nuovamente la formula di
Simpson (1.17) su ciascuno dei due sottointervalli [α, (α + β)/2] e [(α + β)/2, β], ciascuno di ampiezza h/2:
If (α, β) − Sf,2 (α, β) = −
(h/2)5 (4)
(f (ξ) + f (4) (η)),
90
dove ξ ∈ (α, (α + β)/2), η ∈ ((α + β)/2, β) e
Sf,2 (α, β) = Sf (α, (α + β)/2) + Sf ((α + β)/2, β).
Nell’ulteriore ipotesi che la funzione f (4) (x) non sia troppo variabile su [α, β] (e ciò
non è sempre vero, in generale) possiamo assumere f (4) (ξ) ' f (4) (η) e concludere che
If (α, β) − Sf,2 (α, β) = −
1 h5 (4)
f (ξ),
16 90
(2.2)
24
con una riduzione di un fattore 16 rispetto all’errore (2.1), corrispondente alla scelta
di un passo doppio. Confrontando la (2.2) e la (2.1), si ricava la stima
h5 (4)
16
f (ξ) ' δf (α, β)
90
15
(2.3)
dove abbiamo posto
δf (α, β) = Sf (α, β) − Sf,2 (α, β).
Confrontando le uguaglianze (2.2) e (2.3) otteniamo, quindi, la stima
|If (α, β) − Sf,2 (α, β)| '
|δf (α, β)|
.
15
Abbiamo dunque ottenuto una stima dell’errore commesso se l’integrale esatto, nella
parte relativa al sottointervallo [α, β], viene approssimato con la quadratura di Simpson composita su quel sottointervallo; tale stima deve essere verificata se la funzione
f è sufficientemente regolare e non troppo variabile nel sottointervallo [α, β]. Eviden-
temente la procedura di calcolo dell’integrale approssimato, che comunque richiederà
di decomporre l’intervallo di integrazione iniziale [a, b] in sottointervalli [α, β] più o
meno spaziati, richiederà uno sforzo computazionale maggiore soltanto nella parte in
cui la funzione subisce forti variazioni, nella maggior parte dei casi abbastanza localizzate. In pratica è più conveniente assumere una stima dell’errore più conservativa,
ad esempio
|If (α, β) − Sf,2 (α, β)| '
|δf (α, β)|
;
10
inoltre, in virtù della proprietà di additività dell’integrale, per garantire un’accuratezza
complessiva su [a, b] pari alla tolleranza ε prefissata, basterà imporre che su ogni sin-
25
a
α Α β
S
N
b
N
β =b
(I)
a
a
S
S
α
α+β
α A 2
N
( II )
b
golo sottointervallo [α, β] ⊂ [a, b] la stima dell’errore δf (α, β) verifichi
|δf (α, β)|
β−α
≤ε
.
10
b−a
(2.4)
L’algoritmo automatico per l’integrazione adattiva può essere organizzato come segue.
Indichiamo con:
1. A l’intervallo di integrazione attivo, ovvero dove si sta calcolando l’integrale;
2. S l’intervallo di integrazione già esaminato, per il quale il test sull’errore
è stato superato con successo;
3. N l’intervallo di integrazione che deve essere ancora esaminato.
All’inizio del processo di integrazione abbiamo N = [a, b], A =N e S =∅, mentre la
situazione al generico passo dell’algoritmo è visualizzata in figura 2.1. Indichiamo con
Js (f ) l’approssimazione della parte di integrale
Rα
a
f (x)dx già calcolata; evidentemente
Js (f ) = 0 all’inizio del processo; se l’algoritmo di calcolo giunge a buon fine, Js (f )
fornisce l’approssimazione cercata di I(f ). Indichiamo inoltre con J(α,β) (f ) l’integrale
approssimato di f sull’intervallo attivo [α, β] . Quest’ultimo è evidenziato in grassetto
26
nella figura 2.1. Ad ogni passo del metodo di integrazione adattiva si procede secondo
lo schema rappresentato in figura 2.1 e cioè:
1. Se il test sull’errore locale (2.4) è superato:
(i) si incrementa Js (f ) della quantità J(α,β) (f ), ovvero Js (f ) = Js (f )+J(α,β) (f );
(ii) si pone S = S ∪ A, A = N (corrispodente al percorso (I) nella figura 2.1),
e β = b.
2. Se il test sull’errore locale (2.4) non è superato:
·
α+β
(j) si dimezza A, ponendo il nuovo intervallo attivo pari a A = α,
2
¸
(corrispondente al percorso (II) nella figura 2.1);
·
(jj) si pone N =
¸
α+β
α+β
, β ∪ N, e β =
;
2
2
(jjj) si procede ad una nuova stima dell’errore.
Onde evitare che il passo di integrazione diventi troppo piccolo, conviene introdurre
un controllo sull’ampiezza di A e segnalare, in caso di eccessiva riduzione, la presenza
di un eventuale punto di singolarità della funzione integranda.
Esercizio 8 Applichiamo l’algoritmo adattivo di Simpson al calcolo del seguente in-
tegrale
Z4
tan−1 (10x)dx = 4 tan−1 (40) + 3 tan−1 (−30) − (1/20) log(16/9) ' 1.54201193
I(f ) =
−3
L’esecuzione del programma sottostante, assumendo tol = 10−4 e hmin = 10−3 , fornisce una approssimazione dell’integrale con un errore assoluto di 2.104·10−5 . L’algoritmo
esegue 77 valutazioni funzionali corrispondenti ad una suddivisione dell’intervallo [a, b]
in 38 sottointervalli di ampiezza non uniforme. La corrispondente formula composita
27
a passo uniforme avrebbe richiesto 128 suddivisioni con un errore assoluto pari a
2.413 · 10−5 .
L’algoritmo adattivo sopra descritto è implementato nel programma sottostante. Tra i
parametri di ingresso, hmin è il valore minimo ammissibile per il passo di integrazione.
In uscita vengono restituiti il valore approssimato dell’integrale integ , il numero totale
di valutazioni funzionali eseguite nf v e l’insieme dei punti di valutazione xf v .
Programma simpadpt : Integratore adattivo con la formula di Simpson
function [integ,xfv,nfv]=simpadpt(a,b,tol,fun,hmin)
integ=0;
level=0;
i=1;
alfa(i)=a;
beta(i)=b;
step=( beta(i)-alfa(i))/4;
nfv=0;
for k=1: 5
x=a+(k-1)*step;
f(i,k)=eval(fun);
nfv=nfv+1;
end
while(i > 0)
S=0;
28
S2=0;
h=(beta(i)-alfa(i))/2;
S=(h/3)*(f(i,1)+4*f(i,3)+f(i,5));
h=h/2;
S2=(h/3)*((f(i,1)+4*f(i,2)+f(i,3));
S2=S2+(h/3)*(f(i,3)+4*f(i,4)+f(i,5));
tolrv=tol*(beta(i)-alfa(i))/(b-a);
errrv=abs(S-S2)/10;
if(errrv > tolrv)
i=i+1;
alfa(i)=alfa(i-1);
beta(i)=(alfa(i-1)+beta(i-1))/2;
f(i,1)=f(i-1,1);
f(i,3)=f(i-1,2);
f(i,5)=f(i-1,3);
len=abs(beta(i)-alfa(i));
if(len >= hmin)
if(len <= 11*hmin)
str=sprintf(’Funzione singolare in x=%12.7f’, alfa(i));
disp(str);
end
step=len/4;
29
x=alfa(i)+step;
f(i,2)=eval(fun);
nfv=nfv+1;
x=beta(i)-step;
f(i,4)=eval(fun);
nfv=nfv+1;
else
xfv=xfv’;
disp(’h e”minore di hmin’)
str=sprintf(’L”integrale sinora calcolato vale%12.7e’, integ);
disp(str);
return
end
else
integ=integ+S2;
level=level+1;
if(level==1)
for k=1: 5
xfv(k)=alfa(i)+(k-1)*h;
end
ist=5;
else
30
for k=1: 4
xfv(ist+k)=alfa(i)+k*h;
end
ist=ist+4;
end
if(beta(i)==b)
xfv=xfv’;
str=sprintf(’L’integrale vale%12.7e’,integ);
disp(str);
return
end
i=i-1;
alfa(i)=beta(i+1);
f(i,1)=f(i+1,5);
f(i,3)=f(i,4);
step=abs(beta(i)-alfa(i))/4;
x=alfa(i)+step;
f(i,2)=eval(fun);
nfv=nfv+1;
x=beta(i)-step;
f(i,4)=eval(fun);
nfv=nfv+1;
31
end
end
32
Capitolo 3
Quadratura numerica per funzione
integranda non limitata o intervallo
di integrazione non limitato
3.1
Introduzione
Col termine ”integrale improprio” si indica generalmente un integrale in cui
la funzione integranda non è limitata nell’intervallo di integrazione; tuttavia, occasionalmente è utile includere in questa categoria il caso di integrali di funzioni
che posseggono una singolarità apparente o eliminabile nell’intervallo di integrazione,
ed anche il caso degli integrali su intevalli infiniti. Gli integrali impropri appaiono
frequentemente in processi computazionali e devono essere trattati mediante metodi
speciali. I primi paragrafi del presente capitolo sono dedicati al caso in cui la funzione
33
integranda non è limitata nell’intervallo di integrazione finito. Assumeremo a volte,
senza perdere in generalità, che l’integrale da valutare sia dato nella forma
R1
f (x)dx,
0
dove la funzione f (x) è continua in 0 < x ≤ 1 ma non in 0 ≤ x ≤ 1; ad esempio f (x) può
essere non limitata nelle vicinanze di x = 0. Il caso degli integrali su domini illimitati,
parimenti trattato nei paragrafi finali del presente capitolo, può essere ricondotto al
caso di integrale di funzione non limitata su intervallo finito mediante una opportuna trasformazione di variabile. Ulteriori metodi di quadratura in questo caso sono
organizzati mediante procedimento limite, o anche mediante l’utilizzo di formule interpolatorie di tipo Gaussiano che assumono come nodi di quadratura rispettivamente
gli zeri dei polinomi ortogonali di Hermite e di Laguerre.
3.2
Idee generali
Per fissare le idee, supponiamo che la funzione integranda f (x) diventi infinita
per x −→ a+ , con a estremo sinistro dell’intervallo di integrazione; l’analisi è del tutto
simile qualora f (x) diventi infinita per x −→ b− . Supponiamo, ad esempio, che la
funzione integranda si possa scrivere nella seguente forma:
f (x) =
φ(x)
, 0<θ<1
(x − a)θ
dove φ(x) è una funzione regolare sull’intervallo chiuso e limitato [a, b] , in modulo non
superiore a M . Sotto tali ipotesi possiamo dare la seguente stima dell’integrale:
Zb
|I(f )| ≤ M lim
t→a+
t
1
(b − a)1−θ
dx
=
M
.
1−θ
(x − a)θ
34
La restrizione su θ assicura, come noto, l’esistenza dell’integrale come limite di integrali di Cauchy-Riemann su sottointervalli che costituiscono una esaustione dell’intervallo
[a, b]. Per il calcolo numerico dell’integrale I(f ) si possono considerare le idee espresse
nel paragrafo 3.2.
3.3
Metodi per il calcolo di integrali impropri
Un primo metodo consiste nel fissare un ² tale che 0 < ² < (b − a) in modo
tale da scrivere l’integrale improprio come somma di due integrali I(f ) = I1 + I2 , dove
I1 =
a+²
R
a
Rb φ(x)
φ(x)
dx
I
=
dx
2
θ
(x − a)θ
a+² (x − a)
L’integrale I2 non presenta più problemi, essendo la funzione integranda regolare
sull’intervallo di integrazione; quindi, per la sua approssimazione si può procedere
con metodi di quadratura classici. Allo scopo di calcolare I1 sostituiamo φ con il suo
sviluppo in serie di Taylor attorno ad x = a arrestato all’ordine p :
φ(x) = φp (x) +
(x − a)p+1 (p+1)
φ
(ξ(x)), p ≥ 0
(p + 1)!
dove abbiamo posto
φp (x) =
p
X
k=0
φ(k) (a)
(x − a)k
.
k!
Si ha allora
I1 (f ) = ε1−θ
p
X
k=0
εk φ(k) (a)
1
+
k!(k + 1 − θ) (p + 1)!
Za+²
(x − a)p+1−θ φ(p+1) (ξ(x))dx
a
(3.1)
35
per cui, sostituendo I1 (f ) con la sommatoria a secondo membro dell’equazione precedente, l’errore corrispondente R1 (f ) è maggiorato da
|R1 (f )| ≤
¯
¯
¯
¯
max ¯φ(p+1) (x)¯ ,
a≤x≤a+²
p ≥ 0.
Supponiamo di fissare ε < 1; se le derivate di φ sono limitate o non aumentano troppo
all’aumentare di p, si ha che R1 (f ) decresce in modulo al crescere di p. L’integrale
relativo alla parte I2 (f ) si può approssimare utilizzando una formula di Newton Cotes
composita a m sottointervalli e n + 1 nodi di quadratura per sottointervallo, con n
pari. Usando le note formule per l’errore nel caso delle formule chiuse con n pari [2,
p. 281]
Rn,m (f ) =
(b − a) Mn n+2 (n+2)
H
f
(ξ),
(n + 2)! nn+3
H=
b−a
,
m
Z
Mn =
0
n
tπn+1 (t)dt < 0
e volendo equidistribuire l’errore δ tra I1 e I2 , si ha
|R2 (f )| ≤ M(n+2) (ε)
dove
(n+2)
M
(b − a − ε) |Mn |
(n + 2)!nn+3
µ
b−a−ε
m
¶n+2
=
δ
2
(3.2)
¯ n+2 ·
¸¯
¯d
¯
φ(x)
¯.
(ε) = max ¯¯ n+2
θ
a+²≤x≤b dx
(x − a) ¯
Si osserva che il valore della costante M(n+2) (ε) cresce rapidamente al tendere a zero
di ε; conseguentemente la (3.2) può comportare un numero di suddivisioni mε = m (ε)
talmente elevato da rendere il metodo in esame di scarsa utilità pratica.
Esempio 9 Si debba calcolare l’integrale generalizzato (noto come integrale di Fresnel)
π
Z2
I(f ) =
0
cos(x)
√ dx.
x
36
Sviluppando in serie di Taylor la funzione integranda nell’intorno dell’origine e applicando il teorema di integrazione per serie, si ottiene
I(f ) =
∞
X
(−1)k
k=0
1
(π/2)2k+1/2 .
(2k)! (2k + 1/2)
Troncando la serie ai primi 10 termini, si ha un valore dell’integrale pari a 1.9549.
Utilizzando la formula di Simpson composita, la stima a priori (3.2) fornisce, ponendo
n = 2, |Mn | =
4
e facendo tendere ε a zero
15
·
¸1/4
´5
0.018 ³ π
mε '
− ε ε−9/2
.
δ
2
Per δ = 10−4 , prendendo ε = 10−2 , si ricava che sono necessarie 1140 suddivisioni
(uniformi) mentre per ε = 10−4 e ε = 10−6 ne servono rispettivamente 2 · 105 e 3.6 · 107 .
Per confronto, eseguendo il programma scritto nel capitolo precedente (integrazione
adattiva con formula di Simpson) con a = ε = 10−10 , hmin = 10−12 e tol = 10−4 , si
ottiene la stima 1.955 dell’integrale con 1057 valutazioni funzionali, corrispondenti a
h πi
528 suddivisioni non uniformi dell’intervallo 0,
.
2
Un altro metodo si ricava decomponendo l’integrale I(f ) nella somma
Zb
I(f ) =
a
φ(x) − φp (x)
dx +
(x − a)θ
Zb
a
φp (x)
dx = I1 + I2
(x − a)θ
con l’ausilio dello sviluppo di Taylor (3.1). Il calcolo esatto dell’integrale I2 , dà la
seguente espressione
I2 = (b − a)1−θ
p
X
(b − a)k φ(k) (a)
k=0
L’integrale I1 si scrive, per p ≥ 0
k!(k + 1 − θ)
.
37
Zb
I1 =
(p+1) (ξ(x))
p+1−θ φ
(x − a)
(p + 1)!
a
Zb
g(x)dx
dx =
a
dove I1 ha una funzione integranda con una singolarità su derivate di ordine più
elevato. La scelta del valore p più opportuno dipenderà quindi dalla formula di
quadratura utilizzata.
3.4
Procedimento con limite
La definizione di base è
Z1
Z1
f (x)dx = lim
t→0+
0
f (x)dx.
t
Analizziamo un primo modo di procedere. Sia 1 > r1 > r2 > ... > 0 una sequenza di
punti che converge a zero, ad esempio rn =
Z1
0
Zr2
Zr1
Z1
f (x)dx =
1
. Scriviamo
2n
r3
r2
r1
f (x)dx + ....
f (x)dx +
f (x)dx +
Ogni integrale
al secondo
membro è un integrale proprio e la valutazione termina
¯r
¯
n+1
¯Z
¯
¯
¯
¯
quando ¯
f (x)dx¯¯≤ ². Questo è un criterio pratico di stop, non corretto teorica¯
¯
rn
mente. Infatti, se applicato all’integrale divergente
rn =
1
indica convergenza, essendo
n
R1 dx
relativamente alla successione
0 x
1
Zn
dx
n+1
= log
x
n
1
(n+1)
e
¯
¯
¯
¯
n
+
1
¯ = log 1 = 0.
lim ¯¯log
n→∞
n ¯
38
Esempio 10
Z1
I=
0
n
1
2
4
8
16
32
40
Z1
dx
1
2
x +x
1
3
, In =
rn
x +x
1
3
, rn = 2−n
Numero di valutazioni della funzione
In
0.28492598
0.47448022
0.68323927
0.81280497
0.84029678
0.84111612
0.84111663
0.84111692
valore esatto
dx
1
2
9
18
44
80
176
344
432
Rrn
In questo caso ogni integrale
f (x)dx è stato calcolato mediante una par-
rn−1
ticolare modifica adattiva delo schema di integrazione di Romberg.
3.5
Troncamento dell’intervallo
Rr
f (x)dx senza
0
¯r
¯
¯R
¯
difficoltà eccessive. Se ¯¯ f (x)dx¯¯ ≤ ε, allora possiamo limitarci a calcolare l’integrale
In qualche caso è possibile ottenere una stima dell’integrale
proprio
R1
0
f (x)dx.
r
Esempio 11 ([5, p. 163]) Stimiamo
Zr
0
g(x)
1
2
1
x + x3
dx
1
1
dove g(x) ∈ C [0, 1] e inoltre |g(x)| ≤ 1. Poichè x 2 ≤ x 3 in [0, 1], abbiamo
¯
¯
¯ g(x) ¯
1
¯ 1
¯
1 ¯ ≤
1 .
¯ 2
x + x3
2x 2
quindi
¯ r
¯
¯Z
¯
Zr
¯
¯ 1
1
g(x)
dx
¯
¯
2
1
1 dx¯ ≤
1 = r .
¯
x2
¯ x2 + x3 ¯ 2
0
0
Ciò suggerisce che dobbiamo prendere r ≤ 10−6 per ottenere una accuratezza di 10−3 .
39
3.6
Cambio di variabile
A volte è possibile trovare un cambio di variabile che elimina la singolarità.
Per esempio, se f (x) ∈ C [0, 1], il cambio di variabile tn = x trasforma l’integrale
Z1
1
x− n f (x)dx,
n≥2
0
Z1
in n f (tn )tn−2 dt che è un integrale proprio. Se l’integrale improprio
0
Z
1
f (x) log xdx,
0
dove f (x) ∈ C [0, 1], f (0) 6= 0 è trattato mediante la sostituzione ovvia t = − log x,
Z∞
otteniamo − te−t f (e−t )dt un integrale su un intervallo di integrazione infinito. In
0
questo caso però il cambiamento di variabile produce soltanto una trasformazione
della difficoltà iniziale in un altro tipo.
Se la funzione integranda è limitata, ma ha un basso ordine di continuità,
può essere preferibile effettuare un cambio di variabile. Per esempio, consideriamo il
caso dell’integrale
Za
I=
p
f (x)x q dx
0
dove p e q sono interi positivi. Ponendo x = tq otteniamo
1
Za q
I = q tp+q−1 f (tq )dt.
0
Altre trasformazioni utili sono:
Z1
−1
f (x)dx
(1 − x2 )1/2
Zπ
=
f (cos t) dt,
0
40
e
Z1
f (x)dx
(x (1 − x))1/2
0
3.7
Zπ/2
¡
¢
=2
f sin2 t dt.
0
Integrazione di tipo interpolatorio
Sia w(x) una funzione con una singolarità nelle vicinanze di x = 0, ma tale
Z1
w(x)xk dx esiste per k = 0, 1, 2, .., n. Per una data sequenza di ascisse
che l’integrale
0
0 < x0 < x1 < x2 < ... < xn ≤ 1, possiamo determinare pesi wi tali che
Z1
w(x)f (x)dx =
n
X
wk f (xk )
k=0
0
se f (x) ∈ Pn . Ciò induce una approssimazione della formula di integrazione
Z
1
w(x)f (x)dx ≈
0
n
X
wk f (xk ).
k=0
1
Esempio 12 Sia w(x) = x− 2 , x0 = 13 , x1 = 23 , x2 = 1. Allora i pesi wk sono determinati
dal seguente sistema lineare

R1 1



w1 + w2 + w3 = 0 x− 2 dx = 2




R 1 −1
2
2
1
2 xdx =
w
+
w
+
w
=
1
2
3
3
3
3
0 x





R

 1 w1 + 4 w2 + w3 = 1 x− 12 x2 dx = 2
9
9
5
0
Otteniamo quindi la formula di quadratura
Z
1
− 12
x
0
14
f (x)dx ≈ f
5
µ ¶
µ ¶
1
8
2
4
− f
+ f (1)
3
5
3
5
che è più conveniente usare nella forma
Z
r
·
− 21
x
0
f (x)dx ≈ r
1
2
14
f
5
µ ¶
µ ¶
¸
1
8
2
4
− f
+ f (1)
3
5
3
5
41
Esempio 13 ([5, p. 141]) Riportiamo una formula dovuta a A. Young
Z
1
−1
3.8
·
µ ¶
µ ¶
¸
1
1
π
f (−1) + 2f −
+ 2f
+ f (1) .
dx ≈
6
2
2
(1 − x2 )1/2
f (x)
Cambio di variabile su intervalli non limitati
La sostituzione x = e−y trasforma l’intervallo 0 ≤ y ≤ ∞ in 0 ≤ x ≤ 1. Quindi
possiamo scrivere la seguente formula
Z∞
Z1
f (y)dy =
0
0
f (− log x)
dx =
x
Z1
0
g(x)
dx
x
(3.3)
che riconduce un integrale su un intervallo illimitato ad un integrale su un intervallo
limitato. Se
g(x)
è limitata nelle vicinanze di x = 0, allora il secondo integrale sarà
x
proprio. Ma se ciò non accade vuol dire che abbiamo cambiato solo il tipo di difficoltà.
Una forma alternativa di questa trasformazione è
R∞
0
e−x f (x)dx =
¢
R1 ¡
1
0 f log x dx. Es-
istono ragioni di natura numerica per le quali la trasformazione (3.3) è valida, se f (y)
è di tipo esponenziale, e |f (y)| ≤ e−ky , 0 ≤ y < ∞. La trasformazione (3.3) è un caso
particolare di una procedura generale. Sia t(x) una qualsiasi funzione appartenente a
C 1 [0, ∞) e monotona in questo intervallo, tale che t(x) soddisfa le condizioni t(0) = 1,
t(∞) = 0, oppure t(0) = 0, t(∞) = 1. Allora abbiamo
Z
Z
∞
1
f (x)dx =
0
0
¯ ¯
¯ dx ¯
f (x(t)) ¯¯ ¯¯ dt.
dt
La formula (3.3) usa la sostituzione t(x) = e−x . Altre possibili sostituzioni sono
t(x) =
x
e t(x) = tanh(x). Come nella ( 3.3), le risultanti integrazioni usualmente
1+x
sono improprie. Altre trasformazioni che sono usualmente utilizzate sono
Z
Z
b
f (x)dx = (b − a)
a
µ
∞
f
0
a + bt
1+t
¶
dt
,
(1 + t)2
42
π
Z2
Z∞
f (x)dx =
0
f (x)
0
sin(x)
dx,
x
a patto che
f (x + π) = f (x) e f (x) = f (−x),
e
π
Z2
Z∞
f (x) cos(x)dx =
0
f (x)
0
sin(x)
dx,
x
a patto che
f (x + π) = −f (x) e f (x) = f (−x),
3.9
Procedimento al limite per intervalli infiniti
La definizione di base
Z
Z
∞
r
f (x)dx = lim
f (x)dx
r→∞ 0
0
suggerisce un primo modo di procedere. Sia 0 < r0 < r1 < ... una sequenza di numeri
tale che lim rn = ∞. Scriviamo
n→∞
Z
Z
∞
f (x)dx =
0
Z
r0
r1
f (x)dx +
0
f (x)dx + ...
(3.4)
r0
Ogni integrale a secondo membro della (3.4) è proprio, e le valutazioni terminano
¯
¯
¯
¯ rn+1
¯
¯ R
f (x)dx¯ ≤ ε. Come constatato in precedenza, questo particolare metodo
quando ¯
¯
¯ rn
di stop non è corretto teoricamente; per esempio l’integrale
Z∞
1
dx
x
43
è divergente, ma il test di stop applicato al caso rn = n indica convergenza. Di solito
il sottointervallo di integrazione finita è frequentemente raddoppiato ad ogni passo,
cioè si pone rn = 2n ; l’idea che supporta questa scelta è che se si usa una sequenza
aritmetica (rn = cn), il contributo di ogni singola addizione ad ogni passo può essere
eccessivamente piccolo, a volte minore di ε già dopo pochi passi, interrompendo cosı̀
il procedimento di calcolo troppo presto.
Esempio 14
Zrn
In =
0
e−x
dx, rn = 2n
1 + x4
n
In
Numero di valutazioni della funzione
0
0.57202582
35
1
0.62745952
52
2
0.63043990
100
3
0.63407761
178
4
0.63407766
322
valore esatto 0.63407783
3.10
Accelerazione della convergenza
Il metodo descritto sopra può essere accelerato se riusciamo ad ottenere uno
sviluppo asintotico di
R∞
r
f (x)dx per r → ∞. Per esempio, usando il primo termine
della (3.6) si ha
Z
r
∞
ce−r
e−x
dx
v
(1 + x4 )
(1 + r4 )
44
da cui applicando l’estrapolazione di Richardson nella forma
In0 =
In Φ(rn+1 ) − In+1 Φ(rn )
,
Φ(rn+1 ) − Φ(rn )
Φ(r) =
e−r
,
(1 + r4 )
rn = 2n
si ottengono i seguenti valori
n
0
1
2
3
In0
0.62996722
0.63046682
0.63047765
0.63047766
Si noti che I10 è più grande di I2 e I20 è quasi uguale a I4 .
3.11
Integrazione Gaussiana su intervalli illimitati
Per l’integrazione sulla semiretta o sull’intera retta, si possono usare formule
interpolatorie di tipo Gaussiano che assumono come nodi di quadratura rispettivamente gli zeri dei polinomi ortogonali di Laguerre e di Hermite.
I polinomi di Laguerre.
Sono polinomi algebrici, ortogonali sull’intervallo [0, +∞) rispetto alla funzione peso w(x)= e−x , cioè se indichiamo con Ln (x) l’n-simo polinomio allora
Z∞
Ln (x)Lm (x)e−x dx = 0,
n 6= m.
0
Essi sono definiti nel seguente modo
Ln (x) = ex
dn −x n
(e x ),
dxn
n ≥ 0,
e soddisfano la seguente relazione ricorsiva a tre termini



 Ln+1 (x) = (2n + 1 − x)Ln (x) − n2 Ln−1 (x), n ≥ 0,


 L−1 = 0,
L0 = 1.
45
I primi polinomi hanno le seguenti espressioni:
k Lk (x)
0 1
1 −x + 1
2 x2 − 4x + 2
3 −x3 + 9x2 − 18x + 6
4 x4 − 16x3 + 72x2 − 96x + 24
5 −x5 + 25x4 − 200x3 + 600x2 − 600x + 120
6 x6 − 36x5 + 450x4 − 2400x3 + 5400x2 − 4320x + 720
7 −x7 + 49x6 − 882x5 + 7350x4 − 29400x3 + 52920x2 − 35280x + 5040
Per ogni funzione f , definiamo ϕ(x) = f (x)ex . Allora
Z
I(f ) =
Z
∞
f (x)dx =
0
∞
e−x ϕ(x)dx,
0
e basterà applicare a quest’ultimo integrale le seguenti formule, note col nome di
formule di quadratura di Gauss-Laguerre, ottenendo, per n ≥ 1 e f ∈ C 2n ([0, +∞))
I(f ) =
n
X
αk ϕ(xk ) +
k=1
(n!)2 (2n)
ϕ
(ξ),
(2n)!
0<ξ<∞
(3.5)
dove i nodi xk , per k = 1, ..., n, sono gli zeri di Ln (x) e i pesi sono dati da
αk =
(n!)2 xk
.
[Ln+1 (xk )]2
Dalla (3.5) si evince che le formule di Gauss-Laguerre integrano esattamente funzioni
f del tipo ϕ(x)e−x , dove ϕ(x) ∈ P2n−1 . In senso generalizzato possiamo affermare che
esse hanno grado di esattezza ottimale, pari a 2n − 1.
46
Esempio 15 Usando la formula di quadratura di Gauss-Laguerre con n =12 per cal-
colare l’integrale I(f ) =
R∞
0
cos2 (x)e−x dx otteniamo il valore 0.5997 con un errore as-
soluto rispetto all’integrale esatto pari a 2.96 · 10−4 . La formula composita del trapezio
richiederebbe 277 nodi per ottenere la stessa accuratezza.
Mediante il cambio di variabile x→x + a nella precedente formula, si ottiene
la formula di Laguerre modificata:
Z
∞
e−x ϕ(x)dx = e−a
a
n
P
αk ϕ(xk + a) +
k=1
(n!)2 (2n)
ϕ
(ξ),
(2n)!
− ∞ < a < ξ < ∞.
(3.6)
I polinomi di Hermite:
Sono polinomi ortogonali sull’intera retta reale rispetto alla funzione peso
2
w(x) = e−x , cioè se indichiamo con Hn (x) l’n-simo polinomio allora
Z∞
2
Hn (x)Hm (x)e−x dx = 0,
n 6= m.
−∞
Essi sono definiti da
2
Hn (x) = (−1)n e−x
dn −x2
(e ),
dxn
n≥0
e si possono generare ricorsivamente nel modo seguente



 Hn+1 (x) = 2xHn (x) − 2nHn−1 (x), n ≥ 0,


 H−1 = 0,
H0 = 1.
I primi polinomi hanno le seguenti espressioni:
47
k Hk (x)
0 1
1 2x
2 4x2 − 2
3 8x3 − 12x
4 16x4 − 48x2 + 12
5 32x5 − 160x3 + 120x
6 64x6 − 480x4 + 720x2 − 120
7 128x7 − 1344x5 + 3360x3 − 1680x
2
Analogamente al caso precedente, posto ϕ(x) = f (x)ex , si ha
Z∞
I(f ) =
Z∞
2
e−x ϕ(x)dx.
f (x)dx =
−∞
−∞
Applichiamo a quest’ultimo integrale le seguenti formule, note col nome di formule di
quadratura di Gauss-Hermite, ottenendo, per n ≥ 1 e f ∈ C 2n (R)
Z∞
I(f ) =
e
−x2
ϕ(x)dx =
n
X
k=1
−∞
√
n! π (2n)
αk ϕ(xk ) + n
ϕ
(ξ)
2 (2n)!
(3.7)
dove i nodi xk per k = 1, ..., n, sono gli zeri di Hn (x) e i pesi sono dati da
√
2n+1 n! π
αk =
.
[Hn+1 (xk )]2
Dalla (3.7) possiamo concludere che le formule di Gauss-Hermite integrano esatta2
mente funzioni f del tipo ϕ(x)e−x , dove ϕ(x) ∈ P2n−1 ; esse hanno pertanto grado di
esattezza ottimale pari a 2n − 1.
48
Bibliografia
[1] F. Costabile, Appunti di Calcolo Numerico con software didattico, Liguori, Napoli, 1992.
[2] A. Quarteroni, R. Sacco, F. Saleri, Matematica Numerica, Springer-Verlag Italia, Milano, 2000.
[3] A. Ralston, A first course in Numerical Analysis, New York, 1965.
[4] P.J. Davis, Interpolation & Approximation, Dover Publications, New York, 1975.
[5] P.J. Davis, P. Rabinowitz, Methods of Numerical Integration, Academic Press, New
York, 1975.
49
Al termine, desidero ringraziare il Dr. Francesco Dell’Accio che con molta serietà e
pazienza mi ha seguito nella stesura di questa tesi.