Formule di quadratura - trapezi, Cavalieri-Simpson ed estrapolazione di Richardson

Molti integrali che servono in ingegneria non hanno una primitiva elementare: $\int e^{-x^2}dx$, $\int\frac{\sin x}{x}dx$, la lunghezza di un'ellisse. Altre volte la funzione non è nemmeno nota in forma chiusa: si hanno solo valori misurati (la forza $F$ di una molla a diversi allungamenti, la pressione $p$ di un cilindro a diverse posizioni dello stantuffo), e il lavoro $W=\int F\,dx$ o $\int p\,dV$ si può stimare solo con una somma pesata di quei valori. Questa nota costruisce le formule più semplici, ne calcola l'errore e mostra come migliorarle senza calcolare nuovi coefficienti.

Prerequisiti: l'integrale come limite di somme in Integrale definito di RiemannL'integrale definito di f limitata su [a,b] è il limite delle somme di Cauchy-Riemann (somme di aree di rettangoli), se esiste finito e uguale per ogni scelta dei punti; equivale a sup delle somme inferiori = inf delle somme superiori. Se f ≥ 0 è l'area del trapezoide, in generale area sopra l'asse meno area sotto. Sono integrabili le funzioni continue, monotone, continue a tratti; la funzione di Dirichlet no.Integrale definito di Riemann →, il teorema della media in Proprietà dell'integrale definito e teorema della mediaL'integrale definito è lineare, monotono (f ≤ g implica ∫f ≤ ∫g), additivo rispetto all'intervallo (con ∫ da a ad a = 0 e ∫ da b ad a = −∫ da a a b) e verifica |∫f| ≤ ∫|f|. Teorema della media integrale: se f è continua su [a,b] esiste c con f(c) = (1/(b−a))∫f, cioè un rettangolo di altezza f(c) ha la stessa area del trapezoide.Proprietà dell'integrale definito e teorema della media →, l'errore di interpolazione in Interpolazione polinomiale - forma di Lagrange, differenze divise ed erroreDati n+1 punti con ascisse distinte esiste un unico polinomio p_n di grado ≤ n che passa per tutti. Forma di Lagrange p = Σ y_i L_i(x) (L_i(x_j) = δ_ij); forma di Newton p = Σ f[x_0..x_k] Π_{j<k}(x − x_j) con le differenze divise (si aggiunge un nodo senza rifare tutto). Errore f(x) − p_n(x) = f^(n+1)(ξ)/(n+1)! · ω(x), ω(x) = Π(x − x_i): dipende da f e dalla scelta dei nodi. Nodi equispaziati con n grande danno il fenomeno di Runge; i nodi di Chebyshev minimizzano max|ω| = 2^(−n) su [−1,1].Interpolazione polinomiale - forma di Lagrange, differenze divise ed errore → e il resto di Lagrange in Formula di Taylor con resto di LagrangeSe f ha n+1 derivate continue su [a,b], l'errore tra f e il suo polinomio di Taylor si scrive esattamente come f^(n+1)(c)(x-x0)^(n+1)/(n+1)! per un c opportuno: serve a stimare l'errore su tutto un intervallo e a scegliere il grado del polinomio.Formula di Taylor con resto di Lagrange →.

1. Idea generale: integrare il polinomio interpolante

Si scelgono $n+1$ nodi $x_0,\dots,x_n$ in $[a,b]$, si costruisce il polinomio $p_n$ di grado $\le n$ che interpola $f$ nei nodi e si pone $$I(f)=\int_a^bf(x)\,dx\ \approx\ \int_a^bp_n(x)\,dx=\sum_{i=0}^nw_i\,f(x_i),\qquad w_i=\int_a^b\ell_i(x)\,dx ,$$ dove $\ell_i$ è l'$i$-esimo polinomio di Lagrange (vale 1 nel nodo $x_i$ e 0 negli altri). I pesi $w_i$ dipendono solo dai nodi e dall'intervallo, non da $f$: per questo, una volta calcolati, la formula si applica a qualsiasi funzione. Si chiama $Q(f)=\sum w_if(x_i)$ la formula e $E(f)=I(f)-Q(f)$ l'errore (positivo se la formula sottostima).

Definizione (grado di esattezza). Una formula di quadratura ha grado di esattezza $r$ se $E(p)=0$ per ogni polinomio di grado $\le r$ ed $E(x^{r+1})\ne0$.

Basta controllare le potenze $1,x,x^2,\dots,x^r$, perché sia $I$ sia $Q$ sono lineari in $f$. Il grado di esattezza dice quanti termini del polinomio di Taylor di $f$ sono integrati senza errore, quindi quanto la formula è precisa per funzioni regolari.

2. Le tre formule semplici

Sia $[a,b]$ l'intervallo, $m=\frac{a+b}2$ il punto medio, $H=b-a$.

Formula del punto medio. $Q_{PM}=H\,f(m)$ (un nodo, $p_0$ costante). Formula del trapezio. $Q_T=\dfrac H2\big[f(a)+f(b)\big]$ (due nodi, $p_1$ retta per $(a,f(a))$ e $(b,f(b))$: l'integrale è l'area del trapezio). Formula di Cavalieri-Simpson. $Q_S=\dfrac H6\big[f(a)+4f(m)+f(b)\big]$ (tre nodi $a,m,b$, $p_2$ parabola).

Perché i pesi di Simpson sono $\frac H6,\frac{4H}6,\frac H6$. Si cercano $w_0,w_1,w_2$ tali che la formula sia esatta su $1,x,x^2$. Per semplicità $[a,b]=[-1,1]$, nodi $-1,0,1$, $H=2$:

• $f=1$: $w_0+w_1+w_2=\int_{-1}^11\,dx=2$;
• $f=x$: $-w_0+w_2=\int_{-1}^1x\,dx=0$, quindi $w_0=w_2$;
• $f=x^2$: $w_0+w_2=\int_{-1}^1x^2dx=\frac23$, quindi $w_0=w_2=\frac13$ e $w_1=2-\frac23=\frac43$.

I pesi $\frac13,\frac43,\frac13$ sono $\frac H6\cdot(1,4,1)$ con $H=2$. Con un cambio di variabile lineare si ottiene la formula su $[a,b]$.

Grado di esattezza. Per la simmetria rispetto al punto medio, tutte le formule viste integrano esattamente le funzioni dispari rispetto a $m$. Per il punto medio: esatta su $1$ e su $x-m$ (dispari), non su $(x-m)^2$ ($\int=\frac{H^3}{12}$ contro $Q=0$): grado 1. Trapezio: esatto su $1$ e $x$ (la retta coincide col suo interpolante), non su $x^2$: grado 1. Simpson: esatto su $1,x,x^2$ per costruzione, ma anche su $x^3$ perché $(x-m)^3$ è dispari rispetto a $m$ (sia $I$ sia $Q$ danno 0); non su $x^4$: grado 3, un grado in più di quello che ci si aspetta da una parabola.

Esempio. Su $[0,1]$ Simpson dà $Q_S(x^3)=\frac16\big[0+4\cdot\frac18+1\big]=\frac16\cdot\frac32=\frac14=\int_0^1x^3dx$ (esatto) e $Q_S(x^4)=\frac16\big[0+4\cdot\frac1{16}+1\big]=\frac{5}{24}=0{,}20833\ne\frac15$.

3. L'errore delle formule semplici

Teorema (errore del trapezio). Se $f\in C^2([a,b])$ esiste $\xi\in(a,b)$ con $$E_T=I-Q_T=-\frac{(b-a)^3}{12}\,f''(\xi).$$

Dimostrazione. L'errore di interpolazione della retta $p_1$ è $f(x)-p_1(x)=\frac{f''(\xi_x)}2(x-a)(x-b)$. Il polinomio $(x-a)(x-b)$ è $\le0$ su $[a,b]$ (non cambia segno), quindi per il teorema della media integrale esiste $\xi$ con $$E_T=\int_a^b\frac{f''(\xi_x)}2(x-a)(x-b)\,dx=\frac{f''(\xi)}2\int_a^b(x-a)(x-b)\,dx=\frac{f''(\xi)}2\cdot\Big(-\frac{(b-a)^3}6\Big)$$ (l'integrale si calcola con $t=x-a$: $\int_0^H t(t-H)dt=\frac{H^3}3-\frac{H^3}2=-\frac{H^3}6$). $\blacksquare$

Con lo stesso metodo (e un po' più di lavoro per Simpson, dove il polinomio dei nodi non ha segno costante) si ottengono: $$E_{PM}=+\frac{(b-a)^3}{24}f''(\xi),\qquad E_S=-\frac{(b-a)^5}{2880}f^{(4)}(\xi).$$

Che cosa dicono: se $f''>0$ ($f$ convessa) il trapezio sovrastima ($E_T<0$: la corda sta sopra il grafico) e il punto medio sottostima; l'errore del punto medio è la metà di quello del trapezio e di segno opposto. Per Simpson l'errore dipende dalla derivata quarta: si annulla per i cubici, coerente col grado 3.

Procedura (limite superiore dell'errore). Per maggiorare $|E|$ si cerca $M_2=\max_{[a,b]}|f''|$ (o $M_4$ per Simpson) e si scrive $|E_T|\le\frac{(b-a)^3}{12}M_2$, $|E_S|\le\frac{(b-a)^5}{2880}M_4$.

Esempio. $I=\int_0^1e^xdx=e-1=1{,}718282$. Con $f=e^x$, $f''=f^{(4)}=e^x\le e$ su $[0,1]$.

• Trapezio: $Q_T=\frac{1+e}2=1{,}859141$, $E_T=-0{,}140859$; la stima dà $|E_T|\le\frac e{12}=0{,}2265$ ✓.
• Simpson: $Q_S=\frac16\big[1+4e^{0{,}5}+e\big]=\frac16[1+6{,}594885+2{,}718282]=1{,}718861$, $E_S=-0{,}000579$; la stima dà $|E_S|\le\frac e{2880}=0{,}000944$ ✓.

Simpson è precisa di circa 240 volte il trapezio con soli tre valori della funzione. Per un solo intervallo grande, però, nessuna delle due è affidabile: si suddivide.

4. Formule composite

Si divide $[a,b]$ in $n$ sottointervalli di ampiezza $h=\frac{b-a}n$, nodi $x_k=a+kh$, e si applica la formula semplice a ciascuno. Sommando i trapezi: $$T_h=h\Big[\tfrac12f(x_0)+f(x_1)+\dots+f(x_{n-1})+\tfrac12f(x_n)\Big].$$ (I nodi interni compaiono due volte, una per intervallo adiacente, ciascuna con peso $\frac h2$: in totale $h$.) Per Simpson si prendono gli intervalli a coppie ($n$ pari): ogni coppia $[x_{2j},x_{2j+2}]$ ha punto medio $x_{2j+1}$ e contribuisce $\frac{2h}6[f(x_{2j})+4f(x_{2j+1})+f(x_{2j+2})]$: $$S_h=\frac h3\Big[f(x_0)+4f(x_1)+2f(x_2)+4f(x_3)+\dots+2f(x_{n-2})+4f(x_{n-1})+f(x_n)\Big].$$ I pesi alternano $4,2,4,2,\dots$ con $4$ sui nodi di indice dispari e $2$ su quelli pari interni, e $1$ agli estremi.

Teorema (errore delle formule composite). Se $f\in C^2$ (rispettivamente $C^4$) esistono $\xi,\eta\in(a,b)$ con $$I-T_h=-\frac{(b-a)\,h^2}{12}f''(\xi),\qquad I-S_h=-\frac{(b-a)\,h^4}{180}f^{(4)}(\eta).$$

Perché. Si sommano gli errori dei $n$ sottointervalli (ampiezza $h$): per il trapezio $-\frac{h^3}{12}\sum_{k}f''(\xi_k)=-\frac{h^3}{12}\,n\cdot\overline{f''}$ con $\overline{f''}$ media degli $f''(\xi_k)$, che per il teorema dei valori intermedi è un valore $f''(\xi)$; e $nh^3=(b-a)h^2$. Per Simpson ci sono $\frac n2$ coppie di ampiezza $2h$, ciascuna con errore $\frac{(2h)^5}{2880}=\frac{h^5}{90}$: la somma dà $\frac n2\cdot\frac{h^5}{90}f^{(4)}=\frac{(b-a)h^4}{180}f^{(4)}$.

Il trapezio è di ordine 2 (errore $\propto h^2$), Simpson di ordine 4 ($\propto h^4$). Dimezzando $h$ (cioè raddoppiando le valutazioni di $f$) l'errore del trapezio si divide per $2^2=4$, quello di Simpson per $2^4=16$.

Procedura (scelta del numero di intervalli). Per avere $|I-Q|\le\varepsilon$:

1. si calcola $M_2=\max|f''|$ (o $M_4$);
2. si impone $\frac{(b-a)h^2}{12}M_2\le\varepsilon$ (o $\frac{(b-a)h^4}{180}M_4\le\varepsilon$) e si ricava $h$;
3. si prende $n\ge\frac{b-a}h$ intero (per Simpson anche pari).

Esempio. $\int_0^1e^xdx$ con $\varepsilon=10^{-6}$: $M_2=M_4=e$. Trapezi: $\frac{h^2}{12}e\le10^{-6}\Rightarrow h\le\sqrt{12\cdot10^{-6}/e}=0{,}00210\Rightarrow n\ge476$. Simpson: $\frac{h^4}{180}e\le10^{-6}\Rightarrow h\le\big(180\cdot10^{-6}/e\big)^{1/4}=0{,}0902\Rightarrow n\ge11{,}08$, cioè $n=12$ (pari). Simpson richiede 13 valutazioni contro 477.

Tabella di convergenza. Per $f=e^x$ su $[0,1]$ ($I=1{,}7182818$), con l'errore $I-Q$ e il rapporto con l'errore precedente (il passo si dimezza a ogni riga):

I rapporti tendono a $4$ e $16$: è la verifica sperimentale dell'ordine. In generale, se il rapporto è $\approx2^p$ l'ordine è $p=\log_2(\text{rapporto})$.

5. Estrapolazione di Richardson

Per $f$ abbastanza regolare vale uno sviluppo dell'errore del trapezio in potenze pari di $h$ (formula di Eulero-Maclaurin): $$T_h=I+c_2h^2+c_4h^4+c_6h^6+\cdots,$$ con costanti $c_2,c_4,\dots$ che dipendono da $f$ ma non da $h$. Si scrive la stessa relazione con $\frac h2$: $$T_{h/2}=I+\frac{c_2h^2}4+\frac{c_4h^4}{16}+\cdots .$$ Moltiplicando la seconda per $4$ e sottraendo la prima, il termine in $h^2$ scompare: $$4T_{h/2}-T_h=3I-\frac{3c_4h^4}{4}+\cdots\quad\Longrightarrow\quad R_h=\frac{4T_{h/2}-T_h}3=I-\frac{c_4h^4}4+\cdots .$$

Estrapolazione di Richardson. Se una quantità calcolata con passo $h$ si comporta come $Q(h)=I+Ch^p+O(h^{p+1})$ con $p$ noto, la combinazione $$\boxed{\,I\approx\frac{2^pQ(h/2)-Q(h)}{2^p-1}\,}$$ elimina il termine dominante dell'errore e ha ordine più alto.

Fatto notevole. Per i trapezi ($p=2$), $\frac{4T_{h/2}-T_h}3$ coincide con Simpson $S_{h/2}$. Verifica con un solo intervallo: $T_h=\frac{f_0+f_2}2$ (passo $h$, estremi $f_0,f_2$), $T_{h/2}=\frac{f_0+2f_1+f_2}4$; allora $\frac{4T_{h/2}-T_h}3=\frac13\big[f_0+2f_1+f_2-\frac{f_0+f_2}2\big]=\frac{f_0+4f_1+f_2}6$ ✓. Quindi Simpson è «due trapezi estrapolati», e il suo ordine 4 si spiega con l'eliminazione del termine in $h^2$.

Tabella di Romberg. Si può ripetere: dopo aver eliminato $h^2$ resta $h^4$, che si elimina con fattore $2^4=16$, poi $h^6$ con $2^6=64$ e così via. Con $R_{k,0}=T_{h/2^k}$: $$R_{k,j}=\frac{4^jR_{k,j-1}-R_{k-1,j-1}}{4^j-1}\qquad(j=1,\dots,k).$$ Per $\int_0^1e^xdx$ con $h=1$ (valore esatto $1{,}718281828$):

Con soli 9 valutazioni di $f$ ($n=8$) l'ultima cella ha errore $3\cdot10^{-10}$ (il trapezio con lo stesso $n$ ha $2\cdot10^{-3}$): nessun nuovo coefficiente da calcolare, solo combinazioni di numeri già disponibili.

Stima dell'errore a posteriori. Dalla relazione $Q(h)=I+Ch^p$ si ricava $Q(h)-Q(h/2)=Ch^p\big(1-2^{-p}\big)=Ch^p\,\frac{2^p-1}{2^p}$. L'errore del valore con passo $\frac h2$ è $Q(h/2)-I=\frac{Ch^p}{2^p}=\frac{Q(h)-Q(h/2)}{2^p-1}$, cioè $$I-Q(h/2)\approx\frac{Q(h/2)-Q(h)}{2^p-1}.$$ Non serve conoscere $f''$: basta confrontare due risultati. Per $e^x$: $T_1=1{,}859141$, $T_2=1{,}753931$, $\frac{T_2-T_1}3=-0{,}03507$ contro l'errore vero $-0{,}03565$ ✓.

Ordine ignoto. Con tre valori $Q(h),Q(h/2),Q(h/4)$ si stima $p$ da $$\frac{Q(h)-Q(h/2)}{Q(h/2)-Q(h/4)}\approx2^p .$$

6. Implementazione

import numpy as np

def trapezi(f, a, b, n):
    x = np.linspace(a, b, n + 1)
    y = f(x)
    h = (b - a) / n
    return h * (y.sum() - (y[0] + y[-1]) / 2)

def simpson(f, a, b, n):                      # n deve essere pari
    x = np.linspace(a, b, n + 1)
    y = f(x)
    h = (b - a) / n
    return h / 3 * (y[0] + y[-1] + 4 * y[1:-1:2].sum() + 2 * y[2:-1:2].sum())

f, I = np.exp, np.e - 1                       # integrale di e^x tra 0 e 1
print("  n   trapezi      errore    rapp.   Simpson      errore    rapp.")
et_old = es_old = None
for n in [2, 4, 8, 16]:
    T, S = trapezi(f, 0, 1, n), simpson(f, 0, 1, n)
    et, es = I - T, I - S
    rt = f"{et_old / et:6.2f}" if et_old else "     -"
    rs = f"{es_old / es:6.2f}" if es_old else "     -"
    print(f"{n:3d}  {T:.8f}  {et:9.2e}  {rt}  {S:.8f}  {es:9.2e}  {rs}")
    et_old, es_old = et, es

# Richardson ripetuto (tabella di Romberg): R[k][j]
R = [[trapezi(f, 0, 1, 2**k)] for k in range(4)]
for k in range(1, 4):
    for j in range(1, k + 1):
        R[k].append((4**j * R[k][j - 1] - R[k - 1][j - 1]) / (4**j - 1))
for riga in R:
    print("  ".join(f"{v:.9f}" for v in riga))

Output (le due tabelle di sopra):

  n   trapezi      errore    rapp.   Simpson      errore    rapp.
  2  1.75393109  -3.56e-02       -  1.71886115  -5.79e-04       -
  4  1.72722190  -8.94e-03    3.99  1.71831884  -3.70e-05   15.65
  8  1.72051859  -2.24e-03    4.00  1.71828415  -2.33e-06   15.91
 16  1.71884113  -5.59e-04    4.00  1.71828197  -1.46e-07   15.98
1.859140914
1.753931092  1.718861152
1.727221905  1.718318842  1.718282688
1.720518592  1.718284155  1.718281842  1.718281829

Equivalenti MATLAB: trapz(x, y) per i trapezi su dati tabulati (anche con nodi non equispaziati), integral(@(x) exp(x), 0, 1) per un integrale con controllo automatico dell'errore. In Python np.trapz (np.trapezoid nelle versioni recenti di numpy) fa lo stesso di trapz.

7. Dati tabulati

Quando $f$ è nota solo in punti equispaziati, si usa direttamente $T_h$ o $S_h$ sui valori misurati. Il lavoro di una forza $F(x)$ lungo un percorso è $W=\int_{x_0}^{x_n}F\,dx$ (in joule se $F$ è in newton e $x$ in metri); il lavoro di un gas che espande è $\int p\,dV$. L'errore non si può maggiorare con $M_2$ o $M_4$ (non si conoscono le derivate): si usa il confronto tra due passi, cioè la stima a posteriori di §5. Se i punti sono in numero dispari ($n$ pari) si può usare Simpson, altrimenti Simpson sui primi $n-1$ intervalli e trapezio sull'ultimo (o Simpson $\frac38$, non trattata qui).

Errori da evitare

• Usare Simpson con $n$ dispari: la formula composta lavora su coppie di intervalli.
• Pesi $4,2,4,2$ sbagliati: $4$ sui nodi di indice dispari, $2$ su quelli pari interni, $1$ sugli estremi; la somma dei pesi (dentro la parentesi) deve essere $3n$.
• Dimenticare il fattore $\frac12$ agli estremi nei trapezi.
• Scambiare $h$ e $H$ nell'errore: la formula semplice usa $(b-a)^3$ (o $^5$), la composita usa $(b-a)h^2$ (o $h^4$).
• Credere che dimezzare $h$ dimezzi l'errore: dimezza di $4$ (trapezi) o di $16$ (Simpson); se il rapporto osservato è diverso, la funzione non è abbastanza regolare (derivata che esplode in un estremo) o il passo è troppo grande.
• Applicare Richardson con l'esponente sbagliato: $2^p$ con $p=2$ per i trapezi, $p=4$ per Simpson; usare $p=1$ dove l'errore è $O(h^2)$ peggiora il risultato.

Esercizi svolti: Esercizio - trapezio, punto medio e Simpson semplici su 1 su x, Esercizio - numero di sottointervalli per una tolleranza assegnata, Esercizio - lavoro di una forza da dati misurati con trapezi, Simpson e Richardson, Esercizio - tabella di Romberg e stima dell'ordine di convergenza.

Domande d'esame

1. Che cos'è il grado di esattezza? Qual è quello di trapezio, punto medio e Simpson, e perché? Traccia. Massimo $r$ per cui si integrano esattamente tutti i polinomi di grado $\le r$; 1, 1, 3. Per Simpson, $(x-m)^3$ è dispari rispetto a $m$ e quindi $I=Q=0$.
2. Dimostrare che $I-Q_T=-\frac{(b-a)^3}{12}f''(\xi)$. Traccia. Errore di interpolazione lineare, $(x-a)(x-b)$ di segno costante, teorema della media integrale, $\int_a^b(x-a)(x-b)dx=-\frac{(b-a)^3}6$.
3. Quanti intervalli servono per avere errore $\le10^{-6}$ con i trapezi e con Simpson? Traccia. Si maggiora $|f''|$ (o $|f^{(4)}|$) e si risolve $\frac{(b-a)h^2}{12}M_2\le10^{-6}$ (o $\frac{(b-a)h^4}{180}M_4\le10^{-6}$); $n=\lceil\frac{b-a}h\rceil$, pari per Simpson.
4. Mostrare che $\frac{4T_{h/2}-T_h}3$ è la formula di Simpson. Traccia. Su un intervallo: $\frac{f_0+2f_1+f_2-\frac12(f_0+f_2)}{3}=\frac{f_0+4f_1+f_2}{6}$; è Richardson con $p=2$.