11 Integrazione composita
Nella applicazione delle formule di quadratura è conveniente calcolare l’integrale approssimato usando la formula di quadratura su dei sottointervalli. Definendo
H = \dfrac{b-a}{N},\qquad x_0 = a, \quad x_1=a+H, \quad \ldots\quad x_N = a+NH = b
possiamo spezzare l’integrale sui sottointervalli [x_{k-1},x_k] come segue:
\begin{aligned} \int_a^b f(x)\,\mathrm{d}x &= \sum_{k=1}^N \int_{x_{k-1}}^{x_k} f(x)\,\mathrm{d}x \\ &= \sum_{k=1}^N \Big( \mathcal{I}(f;x_{k-1},x_k)+\mathcal{E}(f;x_{k-1},x_k)\Big) \\ &= \underbrace{\sum_{k=1}^N\mathcal{I}(f;x_{k-1},x_k)}_{\text{formula composita}} +\underbrace{\sum_{k=1}^N\mathcal{E}(f;x_{k-1},x_k)}_{\text{errore globale}} \end{aligned}
ancora meglio è definire nodi e pesi con stima dell’errore in un intervallo di riferimento standard facendo un cambio di variabile. L’intervallo di riferimento è l’intervallo [-1,1]. Facendo il cambio di variabile:
\begin{aligned} x &= \dfrac{x_{k-1}+x_k}{2}+z\dfrac{H}{2},\\ \mathrm{d}x &= \dfrac{H}{2}\mathrm{d}z,\\ g_k(z) &= f\left(\dfrac{x_{k-1}+x_k}{2}+z\dfrac{H}{2}\right) \end{aligned}
e quindi
\int_a^b f(x)\,\mathrm{d}x = \sum_{k=1}^N \int_{x_{k-1}}^{x_k} f(x)\,\mathrm{d}x = \dfrac{H}{2}\sum_{k=1}^N \int_{-1}^1 g_k(z)\,\mathrm{d}z
La formula di quadratura per l’intervallo [-1,1] (grazie al teorema Teorema 10.1 avrà la forma
\int_{-1}^1 g(z)\mathrm{d}z = \sum_{j=1}^m w_j g(z_j) + \underbrace{\color{blue}C_q\,g^{(q)}(\omega)}_{\text{errore}},\qquad \omega\in(-1,1)
dove C_q è una costante che dipende dalla formula di quadratura, e quindi
\begin{aligned} \int_a^b f(x)\,\mathrm{d}x &= \dfrac{H}{2}\sum_{k=1}^N \int_{-1}^1 g_k(z)\,\mathrm{d}z \\ &= \underbrace{\dfrac{H}{2}\sum_{k=1}^N \left(\sum_{j=1}^m w_j g_k(z_j)\right) }_{\text{formula composita}} +\underbrace{\dfrac{H}{2} \sum_{k=1}^N{\color{blue}C_q\,g_k^{(q)}(\omega_k)}}_{\text{errore globale}} \\ &= \sum_{k=1}^N \left(\sum_{j=1}^m \left(\dfrac{H}{2}w_j\right) f\left(\dfrac{x_{k-1}+x_k}{2}+z_j\dfrac{H}{2}\right)\right) +\dfrac{NH}{2}\left(\dfrac{1}{N}\sum_{k=1}^N {\color{blue}C_q\,g_k^{(q)}(\omega_k)}\right) \end{aligned}
osservando
\begin{aligned} g_k^{(q)}(\omega_k) &= \left(\dfrac{\mathrm{d}}{\mathrm{d}x}\right)^q f\left(\dfrac{x_{k-1}+x_k}{2}+z\dfrac{H}{2}\right) \\ &= \dfrac{H^q}{2^q} f^{(q)}\left(\dfrac{x_{k-1}+x_k}{2}+z\dfrac{H}{2}\right) \end{aligned}
e definendo i nodi z_{k,j} = \dfrac{x_{k-1}+x_k}{2}+z_j\dfrac{H}{2}
e i pesi
W_j=\dfrac{H}{2}w_j
otteniamo la formula di quadratura composita
\begin{aligned} \int_a^b f(x)\,\mathrm{d}x &= \sum_{k=1}^N \sum_{j=1}^m W_j f(z_{k,j}) +C_q\dfrac{b-a}{2}\dfrac{H^q}{2^q} \left( \dfrac{1}{N}\sum_{k=1}^N f^{(q)}\left(\dfrac{x_{k-1}+x_k}{2}+\omega_j\dfrac{H}{2}\right) \right) \\ &= \underbrace{\sum_{k=1}^N \sum_{j=1}^m W_j f(z_{k,j})}_{\text{formula composita}} + \underbrace{C_q(b-a)\dfrac{H^q}{2^{q+1}} f^{(q)}(\omega)}_{\text{errore globale}} ,\qquad \omega\in(a,b) \end{aligned}
11.1 Stima a priori degli intervalli
Quando si usa una formula di quadratura (composita) è importante dare una stima a priori del numero di intervalli in cui suddividere l’intervallo [a,b] di integrazione in modo da mantenere l’errore sotto una tolleranza assegnata. Cioè dato \epsilon tolleranza assegnata vogliamo calcolare N in modo che
\left\lvert\int_a^b f(x)\,\mathrm{d}x- \underbrace{\sum_{k=1}^N \sum_{j=1}^m W_j f(z_{k,j})}_{\text{formula composita}}\right\rvert = \left\lvert\underbrace{C_q(b-a)\dfrac{H^q}{2^{q+1}} f^{(q)}(\omega)}_{\text{errore globale}}\right\rvert \leq \epsilon
data la costante
M_q = \max_{x\in[a,b]} \left\lvert f^{(q)}(x)\right\rvert
otteniamo
\left\lvert C_q(b-a)\dfrac{H^q}{2^{q+1}} f^{(q)}(\omega)\right\rvert \leq \left\lvert C_q\right\rvert(b-a)\dfrac{H^q}{2^{q+1}} \left\lvert f^{(q)}(\omega)\right\rvert \leq \left\lvert C_q\right\rvert(b-a)\dfrac{H^q}{2^{q+1}} M_q
scegliendo H in modo che soddisfi
\left\lvert C_q\right\rvert(b-a)\dfrac{H^q}{2^{q+1}} M_q \leq \epsilon
allora possiamo garantire che l’errore globale è sotto la tolleranza \epsilon. Dalla identità NH=(b-a) possiamo ricavare una diseguaglianza per N:
\left\lvert C_q\right\rvert(b-a)\dfrac{\left(\dfrac{b-a}{N}\right)^q}{2^{q+1}} M_q \leq \epsilon, \qquad\Rightarrow\qquad \left\lvert C_q\right\rvert\dfrac{(b-a)^{q+1}}{\epsilon\,2^{q+1}} M_q \leq N^q
da cui
per cui il più piccolo intero N che soddisfa 11.1 può essere usato nella formula di quadratura composita per calcolare l’integrale con la tolleranza assegnata.
11.2 Trapezi e Simpson composite
La formula dei trapezi composita si ottiene approssimando l’integrale con polinomi lineari a tratti e prende la forma
\begin{aligned} \int_a^b f(x)\,\mathrm{d}x &= \sum_{k=1}^N \dfrac{H}{2}(f(x_{k-1})+f(x_{k}))-\dfrac{b-a}{12}H^2 f''(\omega) \\ &= \underbrace{\dfrac{H}{2}(f(a)+f(b))+H\sum_{k=1}^{N-1}f(x_{k})}_{ \text{formula composita}}- \underbrace{\dfrac{b-a}{12}H^2 f''(\omega)}_{ \text{errore} } \end{aligned}
Mentre la formula composita di Simpson1 ottenuta approssimando l’integrale con polinomi quadratici a tratti (parabole) prende la forma
\int_a^b f(x)\,\mathrm{d}x = \sum_{k=1}^N \dfrac{H}{6}(f(x_{k-1})+4f(x_{k-1/2})+f(x_{k}))-\dfrac{b-a}{2880}H^4 f''''(\omega)
dove x_{k-1/2} = (x_{k-1}+x_k)/2. Nel caso di Simpson si possono definire i nodi della formula composita come segue:
n = 2N, \quad h = \dfrac{b-a}{n} = \dfrac{b-a}{2N} = \dfrac{H}{2}, \qquad \tilde{x}_k = a+kh, \quad k=0,1,\ldots,n
Con queste definizioni, la formula di Simpson diventa
\int_a^b f(x)\,\mathrm{d}x = \sum_{k=1,3,5,\ldots}^{n-1} \dfrac{h}{3}(f(\tilde{x}_{k-1})+4f(\tilde{x}_k)+f(\tilde{x}_{k+1}))-\dfrac{b-a}{180}h^4 f''''(\omega)
che di solito si scrive nella forma seguente
\begin{aligned} \int_a^b f(x)\,\mathrm{d}x &= \dfrac{h}{3}\Big( f(a)+4f(\tilde{x}_1)+2f(\tilde{x}_2)+4f(\tilde{x}_3)+2f(\tilde{x}_4)+\cdots+\\ & \qquad\qquad+2f(\tilde{x}_{n-2})+4f(\tilde{x}_{n-1})+f(b) \Big) -\dfrac{b-a}{180}h^4 f''''(\omega) \end{aligned}
dove si riconosce il pattern dei pesi [1,4,2,4,2,\ldots,2,4,1] della formula di quadratura.
In sintesi
- Invece di alzare il grado si suddivide l’intervallo e si applica una formula di basso grado su ogni sottointervallo: è la strategia composita.
- L’errore passa da \mathcal{O}(h^{p}) locale a \mathcal{O}(h^{p-1}) globale sull’intervallo, perché i contributi si sommano su \sim (b-a)/h pezzi.
- Trapezi composita è \mathcal{O}(h^2), Simpson composita \mathcal{O}(h^4): da qui la stima a priori del numero di sottointervalli necessari a una tolleranza data.
- La stima a priori richiede una limitazione delle derivate di f; quando non è disponibile si passa alle stime a posteriori e all’adattività.
11.3 Scelta a priori del numero di sottointervalli
Le formule di errore delle regole composite permettono di dimensionare la discretizzazione prima di calcolare l’integrale. Fissata una tolleranza assoluta \mathrm{tol}>0, cerchiamo il più piccolo intero N\geq1 per cui la maggiorazione teorica dell’errore non supera \mathrm{tol}. Questo è il minimo numero di sottointervalli certificato dalla stima, non necessariamente il minimo che produce effettivamente un errore inferiore alla tolleranza.
Ponendo L=b-a, H=L/N e
M_2\geq\max_{[a,b]}|f''(x)|,\qquad M_4\geq\max_{[a,b]}|f^{(4)}(x)|,
le formule composite danno
|E_T|\leq\frac{L^3M_2}{12N^2},\qquad |E_S|\leq\frac{L^5M_4}{2880N^4}.
Pertanto, per M_2>0 e M_4>0, arrotondando sempre per eccesso,
\boxed{N_T=\left\lceil\sqrt{\frac{L^3M_2}{12\,\mathrm{tol}}}\right\rceil}, \qquad \boxed{N_S=\left\lceil\left(\frac{L^5M_4}{2880\,\mathrm{tol}}\right)^{1/4}\right\rceil}.
Il simbolo \lceil x\rceil (parentesi quadre con il tratto inferiore omesso) indica l’arrotondamento all’intero superiore, detto ceiling o ceil: è il più piccolo intero maggiore o uguale a x. Ad esempio, \lceil 15.051\rceil=16, mentre \lceil 16\rceil=16.
Se la maggiorazione della derivata è nulla, la rispettiva regola è esatta e basta un solo intervallo.
Qui N_S conta gli intervalli di Simpson di ampiezza H: ciascun intervallo è diviso dal punto medio in due sottointervalli di ampiezza h=H/2. Se si conta invece il numero n di sottointervalli, allora n=2N_S deve essere pari. Le maggiorazioni richiedono rispettivamente f\in C^2([a,b]) e f\in C^4([a,b]).
Quattro esercizi su questo dimensionamento — trapezi ed esponenziale, Simpson e seno, trapezi e funzione razionale, Simpson e polinomio — sono svolti nel volume degli esercizi. Qui ne resta uno solo, che confronta fra loro tutte le regole viste.
Esercizio — Confronto tra trapezi, Simpson e Gauss-Legendre
Per approssimare
I=\int_0^1 e^x\,\mathrm{d}x=e-1,\qquad \mathrm{tol}=10^{-8},
determinare il minimo numero di intervalli uniformi certificato dalle stime a priori usando trapezi, Simpson e Gauss-Legendre con m=2,3,4,5 nodi per intervallo. Confrontare anche il numero di valutazioni di f.
Svolgimento
Dividiamo [0,1] in N intervalli di ampiezza H=1/N. Per Simpson ciascun intervallo contiene due sottointervalli di ampiezza h=H/2. Tutte le derivate di f(x)=e^x coincidono con e^x, quindi possiamo usare la maggiorazione M_k=e per ogni ordine k.
Le stime per trapezi e Simpson sono
|E_T|\leq\frac{e}{12N^2},\qquad |E_S|\leq\frac{e}{2880N^4}.
Per Gauss-Legendre con m nodi, il cambio di variabile da [-1,1] a un intervallo di ampiezza H dà la stima locale
|E_{G_m,\mathrm{loc}}|\leq C_m H^{2m+1}M_{2m}, \qquad C_m=\frac{(m!)^4}{(2m+1)((2m)!)^3}.
Sommando i contributi dei N intervalli, con NH=1, otteniamo
|E_{G_m}|\leq\frac{C_m e}{N^{2m}},\qquad N_{G_m}=\left\lceil\left(\frac{C_m e}{\mathrm{tol}}\right)^{1/(2m)}\right\rceil.
Applicando le formule e arrotondando per eccesso si trova:
| Regola | Soglia reale per N | Intervalli N | Valutazioni di f | Maggiorazione dell’errore |
|---|---|---|---|---|
| Trapezi | 4759.4483 | 4760 | 4761 | 9.998\cdot10^{-9} |
| Simpson | 17.5277 | 18 | 37 | 8.991\cdot10^{-9} |
| Gauss-Legendre, 2 nodi | 15.8381 | 16 | 32 | 9.601\cdot10^{-9} |
| Gauss-Legendre, 3 nodi | 2.2645 | 3 | 9 | 1.850\cdot10^{-9} |
| Gauss-Legendre, 4 nodi | 0.7908 | 1 | 4 | 1.529\cdot10^{-9} |
| Gauss-Legendre, 5 nodi | 0.4009 | 1 | 5 | 1.072\cdot10^{-12} |
Le valutazioni sono N+1 per i trapezi e 2N+1 per Simpson, riutilizzando i valori agli estremi comuni; per Gauss-Legendre sono mN, perché i nodi sono interni a ciascun intervallo.
In questo esempio, aumentando il numero di nodi di Gauss diminuisce il numero di intervalli richiesto. Tuttavia, il numero di intervalli da solo non misura il costo: Gauss a 4 e a 5 nodi richiedono entrambi un solo intervallo, ma la regola a 4 nodi soddisfa già la tolleranza con una valutazione in meno. Le quantità in tabella sono maggiorazioni teoriche, non errori effettivi; il confronto dipende dalla funzione e dalla tolleranza scelte.