Quadratura gaussiana e formule composite adattive

Con trapezi e Simpson (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 →) i nodi sono fissati (equispaziati) e si ottimizzano solo i pesi. Ma anche i nodi sono incognite: una formula a $m$ nodi ha $2m$ parametri (nodi e pesi), non $m$. Sfruttandoli tutti si ottiene una formula con grado di esattezza $2m-1$ invece di $m-1$ (o $m$): è la quadratura gaussiana. Nella seconda parte si affronta un problema pratico: una funzione che varia poco quasi ovunque e molto in una zona piccola (un picco, un tratto ripido) non si integra bene con un passo uniforme. La soluzione è un passo che si adatta da solo alla funzione.

1. Quadratura di Gauss-Legendre

Si lavora su $[-1,1]$ e si cerca $$\int_{-1}^1f(t)\,dt\approx\sum_{i=1}^mw_if(t_i)$$ con $2m$ incognite $(t_i,w_i)$ da scegliere in modo che la formula sia esatta per $1,t,t^2,\dots,t^{2m-1}$ (sono $2m$ condizioni).

Caso $m=1$. Esattezza su $1,t$: $w_1=2$, $w_1t_1=0\Rightarrow t_1=0$. È la formula del punto medio, grado $1$.

Caso $m=2$. Incognite $w_1,w_2,t_1,t_2$; condizioni su $1,t,t^2,t^3$: $$w_1+w_2=2,\qquad w_1t_1+w_2t_2=0,\qquad w_1t_1^2+w_2t_2^2=\tfrac23,\qquad w_1t_1^3+w_2t_2^3=0 .$$ Per simmetria si cerca $t_2=-t_1$, $w_1=w_2$: la prima dà $w_1=w_2=1$, la seconda e la quarta sono soddisfatte, la terza dà $2t_1^2=\frac23$, cioè $t_1^2=\frac13$: $$t_{1,2}=\mp\frac1{\sqrt3}=\mp0{,}57735,\qquad w_1=w_2=1 .$$ Con due soli valori di $f$ si integra esattamente ogni cubico. Controllo: $t^4$ dà $Q=2\cdot\frac19=\frac29\ne\frac25$, quindi il grado è esattamente $3$.

Caso $m=3$. Con lo stesso metodo (o con il teorema qui sotto): $t=0,\pm\sqrt{3/5}=\pm0{,}774597$ con pesi $\frac89,\frac59,\frac59$, grado $5$.

I nodi sono le radici dei polinomi di Legendre $P_2(t)=\frac{3t^2-1}2$, $P_3(t)=\frac{5t^3-3t}2$, ... (in generale $P_m$ ha $m$ radici reali distinte in $(-1,1)$). Il perché è nel teorema seguente.

Teorema (D). Siano $t_1,\dots,t_m$ le radici di $P_m$ e $w_i=\int_{-1}^1\ell_i(t)\,dt$ i pesi interpolatori. La formula ha grado di esattezza $2m-1$.

Dimostrazione. Sia $p$ un polinomio di grado $\le2m-1$. Si divide per $P_m$ (grado $m$): $p=q\,P_m+r$ con $q,r$ di grado $\le m-1$. Allora:

• $\int_{-1}^1q\,P_m\,dt=0$, perché $P_m$ è ortogonale a tutti i polinomi di grado $<m$;
• $\int_{-1}^1r\,dt=\sum w_ir(t_i)$, perché la formula con pesi interpolatori è esatta sui polinomi di grado $\le m-1$ (il polinomio interpolante di $r$ in $m$ nodi è $r$ stesso);
• $p(t_i)=q(t_i)P_m(t_i)+r(t_i)=r(t_i)$, perché $P_m(t_i)=0$.

Quindi $\int p=\int r=\sum w_ir(t_i)=\sum w_ip(t_i)$. $\blacksquare$

Inoltre i pesi sono positivi (si vede applicando la formula al quadrato $\ell_j^2$ di grado $2m-2$: $0<\int\ell_j^2=\sum_iw_i\ell_j(t_i)^2=w_j$). I pesi positivi significano nessuna amplificazione degli errori sui dati (come per trapezi e Simpson), mentre le formule di Newton-Cotes con molti nodi equispaziati hanno pesi di segno alterno e diventano instabili.

Altro intervallo. Con il cambio di variabile $x=\frac{a+b}2+\frac{b-a}2t$, $dx=\frac{b-a}2dt$: $$\int_a^bf(x)\,dx\approx\frac{b-a}2\sum_{i=1}^mw_i\,f\Big(\frac{a+b}2+\frac{b-a}2t_i\Big).$$

Esempio. $\int_0^1e^x\,dx=1{,}718282$ con $m=2$. Nodi: $x=\frac12\pm\frac12\cdot0{,}57735=0{,}211325;\ 0{,}788675$. Valori: $e^{0{,}211325}=1{,}235314$, $e^{0{,}788675}=2{,}200479$. $$Q=\frac{1-0}2\,(1\cdot1{,}235314+1\cdot2{,}200479)=\frac{3{,}435793}2=1{,}717896,\qquad I-Q=3{,}85\cdot10^{-4}.$$ Per confronto Simpson con tre valori ha errore $5{,}79\cdot10^{-4}$: la formula gaussiana con meno valori è più precisa. Con $m=3$ ($x=0{,}112702;\,0{,}5;\,0{,}887298$, valori $1{,}119298;\,1{,}648721;\,2{,}428560$): $Q=\frac12\big(\frac59\cdot1{,}119298+\frac89\cdot1{,}648721+\frac59\cdot2{,}428560\big)=1{,}718281$, errore $8{,}2\cdot10^{-7}$.

Errore. Per $f\in C^{2m}$: $$E_m=\frac{(b-a)^{2m+1}(m!)^4}{(2m+1)\,[(2m)!]^3}\,f^{(2m)}(\xi).$$ Per $m=2$ la costante è $\frac{16}{5\cdot24^3}=\frac1{4320}$; l'errore dipende da $f^{(4)}$, come Simpson (con costante $\frac1{2880}$) ma con un valore in meno. Per $m=3$ è $\frac1{2016000}$ con $f^{(6)}$.

Pro e contro. Più precisa a parità di valutazioni, converge per ogni $f$ continua quando $m\to\infty$ (pesi positivi) e dà accuratezza altissima per funzioni analitiche. Però i nodi di $m$ e $m+1$ non si ricoprono: per controllare l'errore bisogna rifare tutte le valutazioni. Per questo i programmi (integral di MATLAB, quad di SciPy) usano le coppie di Gauss-Kronrod, che aggiungono nodi a una formula gaussiana per ottenere una seconda stima, e la differenza dà l'errore.

2. Formule composite gaussiane

Come per i trapezi si divide $[a,b]$ in $n$ sottointervalli di ampiezza $h$ e si applica Gauss $m$ nodi su ciascuno. L'errore locale è $O(h^{2m+1})$, la somma di $n$ errori dà $$I-Q_{n,m}\approx\frac{(b-a)\,h^{2m}\,(m!)^4}{(2m+1)\,[(2m)!]^3}f^{(2m)}(\xi):$$ ordine $2m$. Per $m=2$ l'ordine è 4 come Simpson composito, ma con costante $\frac1{4320}$ invece di $\frac1{180}$. A parità di valutazioni $N$ la differenza si riduce: Gauss usa $\frac N2$ intervalli di ampiezza $\frac{2(b-a)}N$ e Simpson $N$ sotto-intervalli di ampiezza $\frac{b-a}N$, quindi gli errori valgono $\frac{16}{4320}=\frac1{270}$ e $\frac1{180}$ volte $\frac{(b-a)^5}{N^4}f^{(4)}$: Gauss è circa $1{,}5$ volte più preciso.

Esempio (numero di intervalli). $\int_0^1e^x\,dx$ con $m=2$ e $\varepsilon=10^{-6}$: $\frac{(b-a)h^4}{4320}e\le10^{-6}\Rightarrow h^4\le\frac{4320\cdot10^{-6}}e=1{,}589\cdot10^{-3}\Rightarrow h\le0{,}1997\Rightarrow n\ge5{,}008$, quindi $n=6$ (12 valutazioni; Simpson ne richiede 13, i trapezi 477). In realtà $n=5$ dà già errore $6{,}4\cdot10^{-7}$: la maggiorazione è prudente.

Dimezzando $h$ l'errore si divide per $2^{2m}$ ($=16$ per $m=2$): rapporti osservati $15{,}6;\ 15{,}9;\ 16{,}0$ per $n=1\to2\to4\to8$.

3. Formule adattive

Una funzione come $f(x)=\frac1{x+0{,}1}$ su $[0,1]$ varia moltissimo vicino a $0$ ($f(0)=10$, $f(1)\approx0{,}9$) e quasi niente altrove. Un passo uniforme abbastanza piccolo per il tratto ripido è sprecato sul resto. L'idea: partire da un solo intervallo, stimare l'errore e dividere solo dove non basta.

Stima dell'errore con due Simpson. Su $[a,b]$ con punto medio $c$ si calcolano $$S_1=\frac{b-a}6\big[f(a)+4f(c)+f(b)\big],\qquad S_2=S(a,c)+S(c,b)$$ (Simpson su tutto l'intervallo e sulle due metà, con passo dimezzato). Poiché Simpson ha ordine $p=4$, per la stima di Richardson (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 →): $$I-S_2\approx\frac{S_2-S_1}{2^4-1}=\frac{S_2-S_1}{15}.$$

Procedura (Simpson adattivo). Dati $[a,b]$ e la tolleranza $\varepsilon$ per questo intervallo:

1. calcola $S_1$ e $S_2$;
2. se $|S_2-S_1|\le15\,\varepsilon$ accetta e restituisci $S_2+\dfrac{S_2-S_1}{15}$ (valore estrapolato, ordine 6);
3. altrimenti dividi in $[a,c]$ e $[c,b]$ e applica la procedura a ciascuna con tolleranza $\frac\varepsilon2$; il risultato è la somma.

La tolleranza dimezzata a ogni divisione garantisce che, se tutte le stime sono corrette, la somma degli errori resti $\le\varepsilon$. Il costo: sull'intervallo iniziale servono $5$ valutazioni ($a$, $b$, $c$ e i punti medi delle due metà). Ogni intervallo figlio, per esempio $[a,c]$, eredita dal padre i valori in $a$, nel suo punto medio e in $c$, quindi richiede solo $2$ nuove valutazioni (i punti medi delle sue due metà). Nell'esempio sotto si esaminano $9$ intervalli: $5+2\cdot8=21$ valutazioni.

Esempio (a mano). $I=\int_0^1\frac{dx}{x+0{,}1}=\ln11=2{,}397895$, $\varepsilon=10^{-3}$. Ogni riga è un intervallo esaminato:

Primo passo a mano: $S_1=\frac16\big[10+4\cdot\frac1{0{,}6}+\frac1{1{,}1}\big]=\frac16\,[10+6{,}666667+0{,}909091]=2{,}929293$; $S_2=S(0,\tfrac12)+S(\tfrac12,1)$ con $S(0,\tfrac12)=\frac{0{,}5}6\big[10+\frac4{0{,}35}+\frac1{0{,}6}\big]=1{,}924603$ e $S(\tfrac12,1)=\frac{0{,}5}6\big[\frac1{0{,}6}+\frac4{0{,}85}+\frac1{1{,}1}\big]=0{,}606803$, somma $2{,}531406$. Il risultato finale somma i cinque intervalli accettati col valore estrapolato: $$Q=2{,}397915,\qquad I-Q=-1{,}9\cdot10^{-5},$$ con $5$ intervalli, $4\cdot5+1=21$ valutazioni, e quindi molto meglio dei $10^{-3}$ richiesti. I punti sono addensati dove serve: gli intervalli accettati sono $[0,0{,}0625]$, $[0{,}0625,0{,}125]$, $[0{,}125,0{,}25]$, $[0{,}25,0{,}5]$, $[0{,}5,1]$ (larghezze $0{,}0625,0{,}0625,0{,}125,0{,}25,0{,}5$). Simpson uniforme per lo stesso errore avrebbe bisogno di $n=64$ (65 valutazioni).

Attenzione. La stima $\frac{S_2-S_1}{15}$ si basa sull'ipotesi che $f^{(4)}$ sia circa costante sull'intervallo: vicino a una singolarità o a un picco molto stretto può sottostimare l'errore, e l'algoritmo accetta troppo presto. Per $\int_0^1\sqrt x\,dx$ con $\varepsilon=10^{-3}$ l'intervallo $[0,\frac12]$ ha $|S_2-S_1|=0{,}0065<15\cdot5\cdot10^{-4}=0{,}0075$ e viene accettato, ma l'errore vero è $3{,}1\cdot10^{-3}$, tre volte la tolleranza: $f^{(4)}\propto x^{-7/2}$ non è limitata in $0$. Per questo i programmi aggiungono un limite al numero di divisioni e segnalano gli intervalli che non convergono.

4. Implementazione

import numpy as np

def gauss_legendre(f, a, b, m):
    t, w = np.polynomial.legendre.leggauss(m)     # nodi e pesi su [-1, 1]
    x = (a + b) / 2 + (b - a) / 2 * t              # cambio di intervallo
    return (b - a) / 2 * np.sum(w * f(x))

def gauss_composita(f, a, b, n, m):
    xs = np.linspace(a, b, n + 1)
    return sum(gauss_legendre(f, xs[i], xs[i + 1], m) for i in range(n))

def simpson_adattiva(f, a, b, eps, intervalli=None):
    """Restituisce (valore, lista degli intervalli accettati)."""
    c = (a + b) / 2
    S = lambda p, q: (q - p) / 6 * (f(p) + 4 * f((p + q) / 2) + f(q))
    S1, S2 = S(a, b), S(a, c) + S(c, b)
    if abs(S2 - S1) <= 15 * eps:
        return S2 + (S2 - S1) / 15, [(a, b)]
    vl, il = simpson_adattiva(f, a, c, eps / 2)
    vr, ir = simpson_adattiva(f, c, b, eps / 2)
    return vl + vr, il + ir

print("Gauss su e^x in [0,1], errore I - Q:")
I = np.e - 1
for m in [1, 2, 3]:
    print(f"  m = {m}: {gauss_legendre(np.exp, 0, 1, m):.9f}   errore {I - gauss_legendre(np.exp, 0, 1, m):.2e}")

f = lambda x: 1 / (x + 0.1)
v, iv = simpson_adattiva(f, 0, 1, 1e-3)
print(f"Simpson adattiva: {v:.6f}  (esatto {np.log(11):.6f}),", len(iv), "intervalli,", 4 * len(iv) + 1, "valutazioni")
print([(round(a, 4), round(b, 4)) for a, b in iv])

Output:

Gauss su e^x in [0,1], errore I - Q:
  m = 1: 1.648721271   errore 6.96e-02
  m = 2: 1.717896378   errore 3.85e-04
  m = 3: 1.718281004   errore 8.24e-07
Simpson adattiva: 2.397915  (esatto 2.397895), 5 intervalli, 21 valutazioni
[(0, 0.0625), (0.0625, 0.125), (0.125, 0.25), (0.25, 0.5), (0.5, 1)]

Equivalenti MATLAB: integral(@(x) 1./(x+0.1), 0, 1, 'AbsTol', 1e-3) (Gauss-Kronrod adattivo), quad (Simpson adattivo, versioni vecchie), quadgk (Gauss-Kronrod). MATLAB non ha una funzione per nodi e pesi di Gauss-Legendre: si tabulano o si calcolano come autovalori di una matrice tridiagonale (eig). In Python numpy.polynomial.legendre.leggauss(m) dà nodi e pesi, scipy.integrate.quad è l'equivalente di integral.

Errori da evitare

• Dimenticare il fattore $\frac{b-a}2$ nel cambio di intervallo (i pesi sono relativi a $[-1,1]$, lunghezza 2).
• Usare i nodi $t_i$ direttamente come $x_i$ senza la trasformazione affine.
• Credere che il grado di esattezza sia $m$: è $2m-1$ ($m=2$: grado 3, non 2).
• Riusare nodi gaussiani di $m$ per $m+1$: non sono annidati.
• Accettare ciecamente la stima adattiva vicino a punti dove $f$ non è liscia.
• Non dimezzare $\varepsilon$ a ogni divisione (o, in alternativa, non dividere $15\varepsilon$ di conseguenza).

Esercizi svolti: Esercizio - Gauss a due e tre punti su 1 su x e confronto con Simpson, Esercizio - periodo del pendolo con grande ampiezza e quadratura gaussiana, Esercizio - Simpson adattiva a mano e confronto con il passo uniforme, Esercizio - nodi e pesi di Gauss a due punti, composita e numero di intervalli.

Domande d'esame

1. Qual è il grado di esattezza della formula di Gauss a $m$ nodi? Perché non si può fare meglio? Traccia. $2m-1$; con $m$ nodi e $m$ pesi ci sono $2m$ parametri e $2m$ condizioni (monomi $1,\dots,t^{2m-1}$). Non esiste una formula a $m$ nodi esatta anche su $\prod(t-t_i)^2$ di grado $2m$, il cui integrale è $>0$ mentre $Q=0$.
2. Dimostrare che con i nodi radici di $P_m$ la formula è esatta fino al grado $2m-1$. Traccia. Divisione $p=qP_m+r$, ortogonalità di $P_m$ ai polinomi di grado $<m$, esattezza interpolatoria su $r$, $P_m(t_i)=0$.
3. Ricavare nodi e pesi di Gauss a due punti. Traccia. Quattro condizioni su $1,t,t^2,t^3$; per simmetria $t_2=-t_1$, $w_1=w_2=1$, $t_1^2=\frac13$.
4. Come funziona un metodo di quadratura adattivo e perché dimezzare la tolleranza? Traccia. $S_1$ contro $S_2$, stima $\frac{S_2-S_1}{15}$, dividere se $>15\varepsilon$; $\varepsilon/2$ a ciascuna metà perché le due metà si sommano.