Esercizio - tabella di Romberg e stima dell'ordine di convergenza

Teoria: Formule di quadratura - trapezi, Cavalieri-Simpson ed estrapolazione di RichardsonUna formula di quadratura approssima $\int_a^b f,dx$ con $\sum w_if(x_i)$ ottenuta integrando il polinomio che interpola $f$ nei nodi. Trapezio $\frac{b-a}2[f(a)+f(b)]$ (grado di esattezza 1, errore $-\frac{(b-a)^3}{12}f''(\xi)$), punto medio $(b-a)f(m)$ (grado 1, errore $+\frac{(b-a)^3}{24}f''(\xi)$), Cavalieri-Simpson $\frac{b-a}6[f(a)+4f(m)+f(b)]$ (grado 3, errore $-\frac{(b-a)^5}{2880}f^{(4)}(\xi)$). Composite con passo $h=\frac{b-a}n$: trapezi errore $-\frac{(b-a)h^2}{12}f''(\xi)$ (ordine 2, rapporto tra errori 4 se $h$ si dimezza), Simpson errore $-\frac{(b-a)h^4}{180}f^{(4)}(\xi)$ (ordine 4, rapporto 16, $n$ pari). Estrapolazione di Richardson: se $Q(h)=I+Ch^p+\dots$ allora $I\approx\frac{2^pQ(h/2)-Q(h)}{2^p-1}$; per i trapezi $p=2$ e $\frac{4T_{h/2}-T_h}3=S_{h/2}$; ripetuta dà la tabella di Romberg. Stima dell'errore a posteriori: $I-Q(h/2)\approx\frac{Q(h/2)-Q(h)}{2^p-1}$.Formule di quadratura - trapezi, Cavalieri-Simpson ed estrapolazione di Richardson →.

Testo. Per $I=\displaystyle\int_0^\pi\sin x\,dx$ si sono calcolati i trapezi compositi $$T_{\pi/2}=1{,}570796\ (n=2),\qquad T_{\pi/4}=1{,}896119\ (n=4),\qquad T_{\pi/8}=1{,}974232\ (n=8).$$

1. Verificare uno di questi valori a mano.
2. Stimare l'ordine di convergenza $p$ dai tre valori, senza usare il valore esatto.
3. Costruire la tabella di Romberg (due colonne di estrapolazione) e dare l'errore stimato dell'ultima riga.
4. Confrontare con il valore esatto $I=2$.
5. (Teoria) Una formula ha ordine $1$: $Q(0{,}5)=1{,}324361$ e $Q(0{,}25)=1{,}512437$. Qual è l'estrapolazione di Richardson?

Soluzione

1. Verifica. Per $n=4$, $h=\frac\pi4$, nodi $0,\frac\pi4,\frac\pi2,\frac{3\pi}4,\pi$ con $\sin=0,\frac{\sqrt2}2,1,\frac{\sqrt2}2,0$: $$T_{\pi/4}=\frac\pi4\Big[\frac02+0{,}707107+1+0{,}707107+\frac02\Big]=\frac\pi4\cdot2{,}414214=1{,}896119 .$$ Coincide col valore dato ✓.

2. Ordine. Si calcolano le differenze successive $$T_{\pi/4}-T_{\pi/2}=0{,}325323,\qquad T_{\pi/8}-T_{\pi/4}=0{,}078113 .$$ Se $T_h=I+Ch^p$, il rapporto vale $\frac{T_{h}-T_{h/2}}{T_{h/2}-T_{h/4}}=2^p$ (le differenze sono $Ch^p(1-2^{-p})\cdot\{1,\,2^{-p}\}$): $$\frac{0{,}325323}{0{,}078113}=4{,}165\ \Longrightarrow\ p=\log_24{,}165=2{,}06\approx2 .$$ L'ordine è 2, come previsto per i trapezi; lo scostamento da 2 esatto è dovuto ai termini $O(h^4)$ non ancora trascurabili con $h$ grande.

3. Romberg. Prima colonna di estrapolazione (fattore $4$, $p=2$): $$R_{1,1}=\frac{4T_{\pi/4}-T_{\pi/2}}3=\frac{4\cdot1{,}896119-1{,}570796}3=\frac{6{,}013680}3=2{,}004560,$$ $$R_{2,1}=\frac{4T_{\pi/8}-T_{\pi/4}}3=\frac{4\cdot1{,}974232-1{,}896119}3=\frac{6{,}000809}3=2{,}000270 .$$ Seconda colonna (fattore $16$, $p=4$): $$R_{2,2}=\frac{16R_{2,1}-R_{1,1}}{15}=\frac{16\cdot2{,}0002692-2{,}0045598}{15}=\frac{32{,}0043072-2{,}0045598}{15}=1{,}9999831 .$$ (Qui servono almeno sette cifre decimali: la sottrazione di due numeri vicini perde cifre significative.)

Errore stimato: i valori $R_{1,1}$ e $R_{2,1}$ hanno ordine $4$, quindi $I-R_{2,1}\approx\frac{R_{2,1}-R_{1,1}}{2^4-1}=\frac{2{,}0002692-2{,}0045598}{15}=-2{,}86\cdot10^{-4}$. Sommare questa correzione a $R_{2,1}$ dà proprio $R_{2,2}$: l'ultima colonna è «il valore più la sua stima d'errore».

4. Confronto con $I=2$. Errori $I-Q$: $T_{\pi/2}$: $0{,}4292$; $T_{\pi/4}$: $0{,}1039$; $T_{\pi/8}$: $0{,}0258$; $R_{2,1}$: $-2{,}7\cdot10^{-4}$ (la stima $-2{,}86\cdot10^{-4}$ del punto 3 è corretta); $R_{2,2}$: $+1{,}7\cdot10^{-5}$. Con gli stessi 9 valori di $\sin$ l'errore è passato da $2{,}6\cdot10^{-2}$ a $1{,}7\cdot10^{-5}$. La stessa stima a posteriori applicata ai trapezi ($p=2$), $\frac{T_{\pi/8}-T_{\pi/4}}3=0{,}0260$, coincide con l'errore vero $0{,}0258$.

5. Ordine 1. Con $p=1$ il fattore è $2^1=2$: $$I\approx\frac{2\,Q(0{,}25)-Q(0{,}5)}{2-1}=2\cdot1{,}512437-1{,}324361=1{,}700513 .$$ (Sono i rettangoli sinistri per $\int_0^1e^xdx=1{,}718282$: l'errore passa da $0{,}206$ a $0{,}0178$. Usare il fattore $4$ di $p=2$ sarebbe sbagliato: dà $\frac{4\cdot1{,}512437-1{,}324361}3=1{,}575129$, con errore $0{,}143$, peggiore.)

Risultati: $p\approx2$; $R_{2,1}=2{,}00027$; $\mathbf{R_{2,2}=1{,}99998}$ (esatto $2$).

Errori tipici. Usare $2^p-1$ al denominatore con $p$ sbagliato; confondere l'indice della colonna con l'ordine ($4^j$, non $2^j$); stimare $p$ con le differenze nell'ordine inverso (rapporto $<1$).