\usepackage{tikz}

13  Stime a posteriori

Autore/Autrice
Affiliazione

Enrico Bertolazzi

University of Trento, Department of Industrial Engineering

Una formula di quadratura ha una formula per l’errore del tipo

\mathcal{E}(f;a,b) = \dfrac{C_q}{2^{p+1}}(b-a)H^p f^{(p)}(z) \tag{13.1}

f^{(p)}(z) è una media campionata della derivata p-esima. Se [a,b] è un intervallo piccolo possiamo approssimare

f^{(p)}(z)\approx \dfrac{1}{b-a}\int_a^b f^{(p)}(x)\mathrm{d}x := F_p

da cui

\mathcal{E}(f;a,b) \approx \dfrac{C_q}{2^{p+1}}(b-a)H^p F_p

Usando un solo intervallo nella formula di quadratura abbiamo

\int_a^b f(x)\mathrm{d}x \approx I_1(f;a,b) + \dfrac{C_q}{2^{p+1}}(b-a)^{p+1} F_p \tag{13.2}

dove I_1(f;a,b) è il valore della formula di quadratura applicata alla funzione f(x) sull’intervallo [a,b] senza suddividerlo (H=b-a). Usando due intervalli nella formula di quadratura abbiamo

\int_a^b f(x)\mathrm{d}x \approx I_2(f;a,b) + \dfrac{C_q}{2^{p+1}}\dfrac{(b-a)^{p+1}}{2^p} F_p \tag{13.3}

dove I_2(f;a,b) è il valore della stessa formula di quadratura applicata alla funzione f(x) sull’intervallo [a,b], suddiviso in due parti (H=(b-a)/2). Sottraendo 13.2 da 13.3 si ottiene

I_2(f;a,b)-I_1(f;a,b) \approx \dfrac{C_q}{2^{p+1}}(b-a)^{p+1}F_p \left(1-\dfrac{1}{2^p}\right)

da cui

F_p \approx \dfrac{2^{2p+1}}{C_q(b-a)^{p+1}(2^p-1)}\left(I_2(f;a,b)-I_1(f;a,b)\right)

Sostituendo questa stima in 13.1, con H=(b-a)/2, si ottiene una stima dell’errore

\begin{aligned} \mathcal{E}(f;a,b)&= \dfrac{C_q}{2^{2p+1}}(b-a)^{p+1} f^{(p)}(z) \\ &\approx \dfrac{C_q}{2^{2p+1}}(b-a)^{p+1} F_p \\ &\approx \dfrac{C_q}{2^{2p+1}}(b-a)^{p+1} \dfrac{2^{2p+1}}{C_q(b-a)^{p+1}(2^p-1)}\left(I_2(f;a,b)-I_1(f;a,b)\right) \\ &\approx \dfrac{I_2(f;a,b)-I_1(f;a,b)}{2^p-1} \end{aligned}

Importante

Schema geometrico della formula di quadratura

Schema geometrico di quadratura: quadratura_08

Stima a posteriori

Data una formula di quadratura di ordine p, si calcoli l’integrale su [a,b] usando un solo intervallo, ottenendo I_1(f;a,b), e poi suddividendo l’intervallo in due, ottenendo I_2(f;a,b). L’errore si può stimare con la formula \mathcal{E}(f;a,b)\approx \dfrac{I_2(f;a,b)-I_1(f;a,b)}{2^p-1}

Si può usare questa stima per costruire formule di quadratura adattative.

13.1 Dalla stima a posteriori a un algoritmo adattativo

La stima ottenuta sopra è una stima asintotica, non un limite rigoroso in assenza di ulteriori ipotesi. Sia Q_1(I) la formula su I=[a,b], e sia Q_2(I) la somma delle due formule sui sottointervalli I_L=[a,m] e I_R=[m,b], dove m=(a+b)/2. Per una regola composita di ordine globale p (trapezi: p=2; Simpson: p=4), assumendo che il termine principale dell’errore sia regolare e non si annulli per cancellazione, vale

I-Q_1\simeq C h^{p+1},\qquad I-Q_2\simeq\frac{C h^{p+1}}{2^p},\qquad \boxed{\widehat E_I=\frac{|Q_2-Q_1|}{2^p-1}}.

In particolare, l’errore stimato della formula raffinata Q_2 è

\widehat E_{T,I}=\frac{|T_2-T_1|}{3},\qquad \widehat E_{S,I}=\frac{|S_2-S_1|}{15}.

Il denominatore 15 riguarda l’errore di S_2, non quello di S_1. L’eventuale correzione di Richardson Q_{\rm corr}=Q_2+(Q_2-Q_1)/(2^p-1) è un’approssimazione ulteriore, ma non è automaticamente certificata dalla stessa stima.

Figura 13.1: Un intervallo e la sua suddivisione in due metà.

Algoritmo ricorsivo con budget locale

Data una tolleranza assoluta globale \mathrm{tol}, all’intervallo iniziale si assegna il budget \tau_{[a,b]}=\mathrm{tol}. Su ogni intervallo I si calcolano Q_1(I), Q_2(I) e \widehat E_I. Se \widehat E_I\le\tau_I si accetta Q_2(I); altrimenti si divide I in due metà e si assegna a ciascuna \tau_I/2. L’algoritmo termina quando tutte le foglie sono accettate (o si raggiunge un limite di sicurezza).

adapt(f, a, b, tau):
    Q1 = quadratura(f, a, b)
    m  = (a+b)/2
    Q2 = quadratura(f, a, m) + quadratura(f, m, b)
    est = abs(Q2-Q1)/(2^p-1)
    if est <= tau:
        return Q2, est
    else:
        L = adapt(f, a, m, tau/2)
        R = adapt(f, m, b, tau/2)
        return L.valore + R.valore, L.stima + R.stima

Se \mathcal L è l’insieme degli intervalli terminali, la costruzione dei budget assicura \sum_{I\in\mathcal L}\tau_I=\mathrm{tol}. Per la disuguaglianza triangolare, se ciascuna stima locale fosse un maggiorante rigoroso, allora

\left|\int_a^b f(x)\,\mathrm dx-\sum_{I\in\mathcal L}Q_2(I)\right| \le\sum_{I\in\mathcal L}|E_I|\le\sum_{I\in\mathcal L}\tau_I=\mathrm{tol}.

Con lo stimatore per differenza di quadrature si ottiene invece un controllo euristico/asintotico: errori con segno opposto possono compensarsi, e in presenza di oscillazioni non risolte Q_1 e Q_2 possono essere vicini pur essendo entrambi inaccurati. Sono utili limiti sulla profondità, sul numero di valutazioni e test contro la stagnazione numerica.

13.2 Perché cercare di equidistribuire l’errore

Sia h_i l’ampiezza dell’intervallo I_i e supponiamo che, localmente, l’errore abbia modulo |E_i|\simeq c(x_i)h_i^{p+1}, con c(x)>0 lentamente variabile. Un reticolo efficiente usa intervalli piccoli dove c(x) è grande e grandi dove c(x) è piccolo.

Teorema 13.1 (Equidistribuzione asintotica per una partizione ottimale) Nel modello continuo, tra le partizioni che minimizzano il numero di intervalli per un budget globale fissato \int_a^b c(x)h(x)^p\,\mathrm dx=\mathrm{tol}, la condizione ottimale è

\boxed{c(x)h(x)^{p+1}=\text{costante}.}

Pertanto gli errori locali (in modulo, nel modello asintotico) sono approssimativamente uguali.

Dimostrazione. La densità di intervalli è 1/h(x), per cui il costo è N\simeq\int_a^b 1/h(x)\,\mathrm dx. L’errore globale stimato è \int_a^b c(x)h(x)^p\,\mathrm dx perché ciascun intervallo contribuisce c(x)h^{p+1} e ce ne sono circa 1/h per unità di lunghezza. Minimizzando il funzionale

\int_a^b\left[\frac1{h(x)}+\lambda c(x)h(x)^p\right]\,\mathrm dx

rispetto a h(x), la condizione di stazionarietà dà

-\frac1{h(x)^2}+\lambda p c(x)h(x)^{p-1}=0 \quad\Longrightarrow\quad c(x)h(x)^{p+1}=\frac1{\lambda p}.

L’equidistribuzione è dunque una proprietà dell’ottimo nel modello continuo, non una condizione che ogni algoritmo a bisezione soddisfa esattamente: l’algoritmo con budget dimezzato garantisce un criterio di accettazione locale, ma può produrre errori terminali molto diversi. Per Simpson p=4, dunque h(x)\propto c(x)^{-1/5}; per i trapezi p=2, h(x)\propto c(x)^{-1/3}. In prossimità di zeri del coefficiente principale o di singolarità servono precauzioni ulteriori.

13.3 Esperimenti numerici: stima ed errore effettivo

Per ciascun esperimento si usa l’algoritmo sopra descritto. Nelle tabelle, Q_{2,i} è la somma delle due quadrature sui semintervalli della foglia I_i, \widehat E_i è lo stimatore per differenza e |E_i|=|\int_{I_i}f-Q_{2,i}| è calcolato confrontando con un integrale di riferimento ad alta accuratezza. Le colonne \tau_i e \widehat E_i/|E_i| permettono di verificare la decisione di accettazione e la qualità dello stimatore. I valori numerici sono arrotondati; l’algoritmo usa valori non arrotondati.

Esempio 1: picco localizzato

Intervallo [0,5], tolleranza \mathrm{tol}=2.0e-04, metodo Simpson. La partizione adattativa contiene 16 intervalli terminali, con 15 suddivisioni. I colori alternati distinguono quattro regioni iniziali dell’intervallo; le linee verticali identificano i sottointervalli terminali.

Figura 13.2: Partizione adattativa e nodi per l’esempio 1; colori alternati per le regioni iniziali.
I_i=[a_i,b_i] h_i Q_{2,i} \widehat E_i |E_i| \tau_i
[0.00000,1.25000] 1.25000 0.0082919552 1.67e-05 3.26e-05 5.00e-05
[1.25000,1.56250] 0.31250 0.0092179833 1.62e-06 2.02e-06 1.25e-05
[1.56250,1.71875] 0.15625 0.011691209 6.92e-07 7.77e-07 6.25e-06
[1.71875,1.87500] 0.15625 0.033317659 5.93e-06 5.11e-06 6.25e-06
[1.87500,1.91406] 0.03906 0.018614328 4.57e-08 4.87e-08 1.56e-06
[1.91406,1.95312] 0.03906 0.027157332 1.82e-07 1.74e-07 1.56e-06
[1.95312,1.99219] 0.03906 0.0360371 9.73e-08 1.07e-07 1.56e-06
[1.99219,2.03125] 0.03906 0.038085486 3.68e-07 3.36e-07 1.56e-06
[2.03125,2.07031] 0.03906 0.030993378 1.78e-07 1.55e-07 1.56e-06
[2.07031,2.10938] 0.03906 0.021732322 9.48e-08 9.71e-08 1.56e-06
[2.10938,2.18750] 0.07812 0.025069459 1.02e-07 1.34e-09 3.13e-06
[2.18750,2.34375] 0.15625 0.020689176 2.71e-06 3.02e-06 6.25e-06
[2.34375,2.50000] 0.15625 0.0085703575 3.04e-07 3.35e-07 6.25e-06
[2.50000,3.12500] 0.62500 0.010890084 1.00e-05 1.61e-05 2.50e-05
[3.12500,3.75000] 0.62500 0.003157862 2.95e-07 3.52e-07 2.50e-05
[3.75000,5.00000] 1.25000 0.0023766176 4.95e-07 6.39e-07 5.00e-05

Bilancio globale: \sum_i\widehat E_i=3.9875e-05; errore effettivo |I-Q|=6.0983e-05; integrale di riferimento I\approx0.305831326199. La stima a posteriori è asintotica, non una maggiorazione rigorosa.

Esempio 2: funzione regolare

Intervallo [0,1], tolleranza \mathrm{tol}=1.0e-06, metodo Simpson. La partizione adattativa contiene 4 intervalli terminali, con 3 suddivisioni. I colori alternati distinguono quattro regioni iniziali dell’intervallo; le linee verticali identificano i sottointervalli terminali.

Figura 13.3: Partizione adattativa e nodi per l’esempio 2; colori alternati per le regioni iniziali.
I_i=[a_i,b_i] h_i Q_{2,i} \widehat E_i |E_i| \tau_i
[0.00000,0.25000] 0.25000 0.28402544 2.40e-08 2.41e-08 2.50e-07
[0.25000,0.50000] 0.25000 0.36469588 3.09e-08 3.09e-08 2.50e-07
[0.50000,0.75000] 0.25000 0.46827879 3.96e-08 3.97e-08 2.50e-07
[0.75000,1.00000] 0.25000 0.60128186 5.09e-08 5.09e-08 2.50e-07

Bilancio globale: \sum_i\widehat E_i=1.4538e-07; errore effettivo |I-Q|=1.4559e-07; integrale di riferimento I\approx1.71828182846. La stima a posteriori è asintotica, non una maggiorazione rigorosa.

Esempio 3: oscillazioni localizzate

Intervallo [0,3], tolleranza \mathrm{tol}=2.0e-02, metodo trapezi. La partizione adattativa contiene 10 intervalli terminali, con 9 suddivisioni. I colori alternati distinguono quattro regioni iniziali dell’intervallo; le linee verticali identificano i sottointervalli terminali.

Figura 13.4: Partizione adattativa e nodi per l’esempio 3; colori alternati per le regioni iniziali.
I_i=[a_i,b_i] h_i Q_{2,i} \widehat E_i |E_i| \tau_i
[0.00000,0.75000] 0.75000 0.75000818 2.73e-06 7.94e-06 5.00e-03
[0.75000,1.50000] 0.75000 0.80949697 4.20e-03 5.58e-02 5.00e-03
[1.50000,1.54688] 0.04688 0.048064812 1.00e-04 8.71e-05 3.13e-04
[1.54688,1.57031] 0.02344 0.019815627 6.14e-05 6.06e-05 1.56e-04
[1.57031,1.59375] 0.02344 0.020272781 3.68e-05 3.65e-05 1.56e-04
[1.59375,1.64062] 0.04688 0.046705993 1.50e-04 1.36e-04 3.13e-04
[1.64062,1.68750] 0.04688 0.049716074 1.54e-04 1.51e-04 3.13e-04
[1.68750,1.87500] 0.18750 0.18838794 7.32e-04 1.69e-03 1.25e-03
[1.87500,2.25000] 0.37500 0.3750912 3.01e-05 1.37e-04 2.50e-03
[2.25000,3.00000] 0.75000 0.75 1.29e-09 2.21e-09 5.00e-03

Bilancio globale: \sum_i\widehat E_i=5.4660e-03; errore effettivo |I-Q|=5.7560e-02; integrale di riferimento I\approx2.99999995924. La stima a posteriori è asintotica, non una maggiorazione rigorosa.