Esercizio - Simpson adattiva a mano e confronto con il passo uniforme

Teoria: Quadratura gaussiana e formule composite adattiveLe formule di Newton-Cotes a $m$ nodi equispaziati hanno grado di esattezza $m-1$ ($m$ se $m$ è dispari). Scegliendo anche i nodi, la formula di Gauss-Legendre a $m$ nodi $\int_{-1}^1f\approx\sum w_if(t_i)$ ha grado $2m-1$: i nodi sono le radici del polinomio di Legendre $P_m$ ($m=2$: $\pm\frac1{\sqrt3}$, pesi 1; $m=3$: $0,\pm\sqrt{3/5}$, pesi $\frac89,\frac59,\frac59$), i pesi sono positivi. Su $[a,b]$: $x_i=\frac{a+b}2+\frac{b-a}2t_i$, $Q=\frac{b-a}2\sum w_if(x_i)$. Errore $\frac{(b-a)^{2m+1}(m!)^4}{(2m+1)[(2m)!]^3}f^{(2m)}(\xi)$ ($m=2$: $\frac{(b-a)^5}{4320}f^{(4)}$). Composita: ordine $2m$ in $h$. Adattiva: si confronta Simpson su $[a,b]$ ($S_1$) con Simpson sulle due metà ($S_2$); se $|S_2-S_1|\le15\varepsilon$ si accetta $S_2+\frac{S_2-S_1}{15}$, altrimenti si divide e si riparte con $\varepsilon/2$ su ciascuna metà: pochi nodi dove $f$ varia poco, molti dove varia molto.Quadratura gaussiana e formule composite adattive →, 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. Si calcola $I=\displaystyle\int_0^1\frac{dx}{x+0{,}2}=\ln6=1{,}791759$ con il metodo di Simpson adattivo e tolleranza globale $\varepsilon=5\cdot10^{-3}$ (criterio: si accetta $[a,b]$ se $|S_2-S_1|\le15\varepsilon_{[a,b]}$ e in tal caso si restituisce $S_2+\frac{S_2-S_1}{15}$; altrimenti si divide a metà con tolleranza $\frac{\varepsilon_{[a,b]}}2$ per parte).

1. Eseguire a mano il primo passo su $[0,1]$.
2. Proseguire sulle due metà, dire quali vengono accettate, dare il valore finale e il numero di valutazioni di $f$.
3. Confrontare l'errore vero con la tolleranza e con Simpson composito a passo uniforme.
4. Ripetere (senza calcoli a mano) con $\varepsilon=5\cdot10^{-4}$ e dire come cambia la suddivisione.
5. (Teoria) Perché si dimezza la tolleranza a ogni divisione? Perché la stima $\frac{S_2-S_1}{15}$ può fallire?

Soluzione

1. Primo passo. Con $f(x)=\frac1{x+0{,}2}$: $f(0)=5$, $f(0{,}5)=\frac1{0{,}7}=1{,}428571$, $f(1)=\frac1{1{,}2}=0{,}833333$. $$S_1=\frac16\big[5+4\cdot1{,}428571+0{,}833333\big]=\frac16\,[11{,}547619]=1{,}924603 .$$ Servono i punti medi delle due metà: $f(0{,}25)=\frac1{0{,}45}=2{,}222222$ e $f(0{,}75)=\frac1{0{,}95}=1{,}052632$. $$S(0,\tfrac12)=\frac{0{,}5}6\big[5+4\cdot2{,}222222+1{,}428571\big]=\frac{0{,}5}6\cdot15{,}317460=1{,}276455,$$ $$S(\tfrac12,1)=\frac{0{,}5}6\big[1{,}428571+4\cdot1{,}052632+0{,}833333\big]=\frac{0{,}5}6\,[1{,}428571+4{,}210526+0{,}833333]=\frac{0{,}5}6\cdot6{,}472430=0{,}539369 .$$ $$S_2=1{,}276455+0{,}539369=1{,}815824,\qquad S_2-S_1=-0{,}108779 .$$ Confronto con la soglia: $15\varepsilon=0{,}075<0{,}108779$, quindi si divide.

2. Le due metà, ciascuna con $\varepsilon'=2{,}5\cdot10^{-3}$ e soglia $15\varepsilon'=0{,}0375$.

• $[0,\frac12]$: $S_1=1{,}276455$ (già calcolato). Servono $f(0{,}125)=\frac1{0{,}325}=3{,}076923$ e $f(0{,}375)=\frac1{0{,}575}=1{,}739130$. $S(0,\frac14)=\frac{0{,}25}6\big[5+4\cdot3{,}076923+2{,}222222\big]=0{,}813746$; $S(\frac14,\frac12)=\frac{0{,}25}6\big[2{,}222222+4\cdot1{,}739130+1{,}428571\big]=0{,}441971$. Quindi $S_2=1{,}255717$ e $S_2-S_1=-0{,}020738$. Poiché $0{,}0207<0{,}0375$, accettato, con valore $1{,}255717-\frac{0{,}020738}{15}=1{,}254335$.
• $[\frac12,1]$: $S_1=0{,}539369$; con $f(0{,}625)=\frac1{0{,}825}=1{,}212121$ e $f(0{,}875)=\frac1{1{,}075}=0{,}930233$ si trova $S(\frac12,\frac34)=0{,}305404$ e $S(\frac34,1)=0{,}233621$, $S_2=0{,}539025$, $S_2-S_1=-0{,}000344$; $0{,}000344<0{,}0375$, accettato, valore $0{,}539025-0{,}000023=0{,}539002$.

Il valore finale è $Q=1{,}254335+0{,}539002=1{,}793337$. Valutazioni: i nodi sono $0,\frac18,\frac14,\frac38,\frac12,\frac58,\frac34,\frac78,1$, cioè 9 (dopo il primo passo, $f(0)$ e $f(1)$ già noti, ogni mezzo intervallo esaminato richiede 2 nuovi punti).

3. Errore e confronto. $I-Q=1{,}791759-1{,}793337=-1{,}58\cdot10^{-3}$, in modulo minore di $\varepsilon=5\cdot10^{-3}$ ✓ (di circa 3 volte: la stima è prudente). Simpson composito uniforme con $n=8$ (9 valutazioni, $h=0{,}125$) ha errore $I-S=-2{,}98\cdot10^{-3}$, con $n=16$ (17 valutazioni) $-2{,}6\cdot10^{-4}$. A parità di 9 valutazioni la versione adattiva è quindi quasi due volte più precisa (grazie alla correzione $\frac{S_2-S_1}{15}$, che è un passo di Richardson). Qui i due sotto-intervalli accettati sono entrambi larghi $\frac12$ perché la tolleranza è larga; con tolleranze più strette la suddivisione diventa non uniforme (punto 4) e il vantaggio cresce.

4. $\varepsilon=5\cdot10^{-4}$ (eseguito con il programma): gli intervalli accettati diventano $[0,\frac18]$, $[\frac18,\frac14]$, $[\frac14,\frac12]$, $[\frac12,1]$, con $4\cdot4+1=17$ valutazioni ed errore $-8{,}5\cdot10^{-6}$. La divisione è concentrata vicino a $0$, dove $f=\frac1{x+0{,}2}$ ha derivate quarte più grandi ($f^{(4)}=\frac{24}{(x+0{,}2)^5}$ è $24\cdot5^5=7{,}5\cdot10^4$ in $x=0$ e $\approx9{,}6$ in $x=1$); a destra l'intervallo resta largo $\frac12$. Simpson uniforme con lo stesso costo ($n=16$, 17 valutazioni) ha errore $-2{,}6\cdot10^{-4}$: rispetta la tolleranza, ma è circa $30$ volte meno preciso; per scendere sotto $10^{-5}$ servono $n=38$ (39 valutazioni: con $n=32$ l'errore è ancora $-1{,}9\cdot10^{-5}$).

5. Teoria. Dimezzare $\varepsilon$: le due metà contribuiscono con errori che si sommano; se ciascuna ha un errore $\le\frac\varepsilon2$ il totale è $\le\varepsilon$. Senza dimezzare, dopo $k$ divisioni l'errore ammesso potrebbe crescere a $2^k\varepsilon$. Perché la stima può fallire: $\frac{S_2-S_1}{15}$ vale se $f^{(4)}$ è circa costante sull'intervallo. Vicino a una singolarità non è così: per $\int_0^1\sqrt x\,dx$ con $\varepsilon=10^{-3}$ l'intervallo $[0,\frac12]$ ha $|S_2-S_1|=0{,}00652<15\cdot5\cdot10^{-4}=0{,}0075$ e viene accettato, ma l'errore vero su quel tratto è $3{,}1\cdot10^{-3}$, oltre il limite $5\cdot10^{-4}$. Per questo si limita sempre la profondità della ricorsione.

Risultati: $\mathbf{Q=1{,}79334}$ (errore $-1{,}6\cdot10^{-3}$), 2 intervalli accettati, 9 valutazioni; con $\varepsilon=5\cdot10^{-4}$: 17 valutazioni, errore $-8{,}5\cdot10^{-6}$.

Errori tipici. Confrontare $|S_2-S_1|$ con $\varepsilon$ invece che con $15\varepsilon$; non dimezzare la tolleranza; dimenticare di riusare i valori già calcolati; non sommare la correzione $\frac{S_2-S_1}{15}$ al risultato accettato.