Integrazione numerica
Appunti di calcolo numerico
Problema dell’integrazione numerica
Data una funzione f:\Bbb{R}\mapsto\Bbb{R} e l’intervallo chiuso e limitato [a,\,b] si vuole stimare numericamente l’integrale
\int_a^b\,f(x)\,\mathrm{d}x
Si parla di stima numerica perché in realtà non si sta chiedendo di calcolare il valore esatto dell’integrale che è definito come differenza dei valori assunti negli estremi di integrazione dalla primitiva della funzione integranda f. Siamo infatti interessati a determinare un valore, ovviamente approssimato, dell’integrale che coinvolga soltanto la conoscenza della funzione integranda f e non della sua primitiva.
Il problema è di particolare interessa perché, pur esistendo finito l’integrale, la funzione integranda potrebbe non ammettere una primitiva esprimibile in maniera analitica, oppure la primitiva potrebbe essere molto complessa da determinare e molto costosa da calcolare. Addirittura, la stessa funzione integranda f potrebbe non essere nota in forma analitica, ma solo tabulata in certi punti dell’intervallo [a,b], rendendo quindi di fatto impossibile il cercarne una primitiva.
Riprendendo una vecchia idea dell’integrale come di una sorta di somma generalizzata, cercheremo di stimare il valore dell’integrale di f per mezzo di espressioni del tipo:
\int_a^b\,f(x)\,\mathrm{d}x \approx \sum_{k=0}^{n} w_k \,f( x_k )
in cui la funzione da integrare è valutata in un certo insieme di punti x_i\in[a,\,b] e l’integrale con una somma pesata con coefficienti w_k opportuni.
Definizione 1 (Formula di quadratura) L’espressione della forma
\mathcal{I}(f;a,b) = \sum_{k=0}^{n} w_k \,f( x_k )
definita da \{x_i,w_i\}_{i=0,1,\ldots,n} si chiama formula di quadratura. I punti \{x_i\}_{i=0,1,\ldots,n} sono i nodi della formula di quadratura ed i coefficienti \{w_i\}_{i=0,1,\ldots,n} sono i pesi.
Ripetiamo ancora che la formula di quadratura fornisce giusto una stima e non il valore dell’integrale, il che implica un errore di approssimazione, che potrà in certe situazioni essere significativo. Ci porremo dunque la questione di quanto una data formula di quadratura sia precisa nella stima dell’integrale di una funzione generica e da quali fattori sia influenzata l’accuratezza del risultato numerico.
Per definire una formula di quadratura è necessario
scegliere i nodi \{x_i\}_{i=0,1,\ldots,n} della formula di quadratura;
scegliere i pesi \{w_i\}_{i=0,1,\ldots,n} della formula di quadratura.
Quindi è ovvio che il criterio con cui si scelgono nodi e pesi è fondamentale nel caratterizzare l’accuratezza di una formula di quadratura.
Inoltre, è ragionevole aspettarsi che la precisione del valore stimato dipenda anche dalla funzione che stiamo cercando di integrare. In altre parole, la stessa formula di quadratura potrebbe produrre un risultato molto accurato – o addirittura il valore esatto – dell’integrale per alcune famiglie di funzioni e fornire un cattivo risultato in altri casi.
Questi argomenti saranno oggetto di discussione nel presente capitolo.
Strategia interpolatoria
approssimiamo la funzione f con un polinomio interpolatore di grado n scegliendo con un criterio da stabilire n+1 nodi di interpolazione con ascisse distinte x_i \in [a,\,b] per i=0,1,\ldots,n; \
stimiamo l’integrale di f calcolando l’integrale del polinomio interpolante:
\int_a^b\,f(x)\,\mathrm{d}x \approx \int_a^b\,p_n(x)\,\mathrm{d}x.
Definizione 2 (Formule di quadratura) Le formule di quadratura costruite con la strategia interpolatoria prendono il nome di formule di quadratura interpolatorie
Classificazione (largamente incompleta)
Formule di Newton-Cotes
i nodi di integrazione sono scelti equidistanti nell’intervallo [a,\,b]; distinguiamotraformule aperte: gli estremi a e b non sono compresi;
formule chiuse: x_0=a e x_n=b;
Formule Gaussiane
i nodi di integrazione sono eventualmente gli estremi dellintervallo di integrazione più gli zeri di polinomi appartenenti a famiglie speciali. Citiamo, per esempio, i polinomi ortogonali che producono le formule di Gauss-Legendre, di Gauss-Lobatto, di Gauss-Radau, etc.
Nella nostra esposizione tratteremo solo il caso delle formule di Newton-Cotes.
Formule di Newton-Cotes
I pesi si determinano immediatamente calcolando l’integrale del polinomio interpolatore. Esaminiamo il caso delle formule chiuse:
\begin{aligned} h &= \frac{b-a}{n}, &\quad&\textrm{passo di integrazione},\\ x_i &=a+i\,h,\qquad i=0,1,\ldots,n,&\quad&\textrm{nodo $i$},\\ x &=a+t\,h,\qquad t\in [0,\,n], &\quad&\textrm{generico punto $x\in[a,\,b]$}. \end{aligned}
Sia p_n(x) il polinomio interpolatore di Lagrange
p_n(x) = f(x_0)\boldsymbol{L}_0(x) + f(x_1)\boldsymbol{L}_1(x) + \ldots + f(x_n)\boldsymbol{L}_n(x),
dove gli \boldsymbol{L}_k(x) sono i polinomi elementari di Lagrange, definiti dalla relazione \boldsymbol{L}_i(x_j)=\delta_{ij}. Ricordiamo la definizione data nel capitolo sull’interpolazione,
\begin{aligned} \boldsymbol{L}_i(x) &= \dfrac{(x-x_0)\ldots(x-x_{i-1} (x-x_{i+1})\ldots( x - x_n )} {(x_i-x_0)\ldots(x_i-x_{i-1})(x-x_{i+1})\ldots( x_i - x_n )} \\ &= \prod_{j=0\atop j\neq i}^{n} \left(\dfrac{ x - x_j }{ x_i - x_j }\right). \end{aligned}
Dato che possiamo scrivere per ogni indice i e j
\begin{aligned} x - x_j &= (a+t\,h) - (a+j\,h) = (t-j)\,h \\ x_i - x_j &= (a+i\,h) - (a+j\,h) = (i-j)\,h, \end{aligned}
si ha dalla definizione dei polinomi di Lagrange un’espressione particolare per nodi di interpolazione equidistanti
\begin{aligned} \int_a^b\,f(x)\,\mathrm{d}x &\approx \int_a^b\,p_n(x) \,\mathrm{d}x, \\ &=\int_a^b\,\sum_{j=0}^{n}f( x_j )\mathcal{L}_j( x )\,\mathrm{d}x, \\ &= \sum_{j=0}^{n}f( x_j )\int_a^b\,\mathcal{L}_j( x )\,\mathrm{d}x, \\ &= \sum_{j=0}^{n}f( x_j )\int_a^b\, \prod_{j=0\atop j\neq i}^{n} \left(\dfrac{ x - x_j }{ x_i - x_j }\right)\,\mathrm{d}x, \\ &=\sum_{j=0}^{n}f( x_j )\,h\int_0^n\,\prod_{j=0\atop j\neq i}^{n} \left(\dfrac{t-j}{i-j}\right)\,\mathrm{d}t. \end{aligned}
Si ottiene quindi la formula di quadratura
\mathcal{I}(f;a,b) = \sum_{k=0}^{n} w_k \,f( x_k ) = \sum_{j=0}^{n}f( x_j ) \,h\int_0^n\,\prod_{j=0\atop j\neq i}^{n} \left(\dfrac{t-j}{i-j}\right)\,\mathrm{d}t.
definita dai nodi e dai pesi
x_i = a+i\,h,\quad i=0,1,\ldots n, \quad w_i = h\int_0^n\,\prod_{j=0\atop j\neq i}^{n}\left(\dfrac{t-j}{i-j}\right)\,\mathrm{d}t.
Accuratezza
Definizione 3 (Errore di integrazione) Si chiama errore di integrazione la differenza tra il valore esatto dell’integrale ed il suo valore stimato con la formula di quadratura
\mathcal{E}(f;a,b) = \int_a^b\,f(x)\,\mathrm{d}x-\sum_{k=0}^{n} w_k \,f( x_k )
Ovviamente \mathcal{E}(f;a,b) deve dipendere dalla funzione f(x) e dall’intervallo di integrazione [a,\,b].
Definizione 4 (Ordine di accuratezza) L’ordine di accuratezza o di precisione della formula di quadratura è dato dal massimo grado dei polinomi che sono integrati esattamente dal metodo. Più correttamente diremo che una formula di quadratura
ha ordine almeno k se integra esattamente tutti i polinomi di grado fino a k;
ha ordine k se ha ordine almeno k ed esiste almeno un polinomio di grado k+1 che non è integrato esattamente.
Metodo dei Coefficienti Indeterminati
Il metodo dei Coefficienti Indeterminati permette di ricondurre il calcolo dei pesi delle formule di quadratura, una volta stabiliti i nodi, alla risoluzione di un problema algebrico.
Osserviamo che l’integrale di un polinomio p_n(x) della forma
p_n ( x ) = a_n x^n + a_{n-1} x^{n-1} + \ldots a_1 x + a_0
si scrive come
\begin{aligned} \int_a^bp_n ( x )\,\mathrm{d}x &= \int_a^b\left(a_n x^n + a_{n-1} x^{n-1} + \ldots + a_1 x + a_0\right)\,\mathrm{d}x \\ &= a_n \int_a^b x^n \,\mathrm{d}x + a_{n-1}\int_a^b x^{n-1}\,\mathrm{d}x+ \cdots+ a_1\int_a^b x \,\mathrm{d}x + a_0\int_a^b\,1\,\mathrm{d}x \end{aligned}
Quindi l’integrale è esatto indipendentemente dai coefficienti a_0,a_1,\ldots,a_n se sono esatti tutti gli integrali della forma
\int_a^b x^i\,\mathrm{d}x,\qquad\textrm{per }i=0,1,\ldots n.
L’osservazione precedente suggerisce il seguente procedimento per la determinazione dei pesi di una formula di quadratura, una volta scelti i nodi (non necessariamente equidistanti).
- Si applica la formula di quadratura a tutti i monomi x^i, per i=0,1,\ldots,n e si impone che il risultato sia uguale al risultato dell’integrazione esatta.
In questo modo si ricava un sistema nelle n+1 equazioni per i=0,1,\ldots,n nelle n+1 incognite w_0,w_1,\ldots w_n.
Praticamente, si deve imporre la condizione
\sum_{j=0}^{n} w_j x_j^i = \int_a^b x^i\,\mathrm{d}x = \frac{b^{i+1}-a^{i+1}}{i+1},\qquad\textrm{per }i=0,1,\ldots,n,
per cui si ottiene un sistema lineare nelle w_i
\begin{aligned} i=0 & \qquad w_0+w_1+\ldots+w_n &=& b-a \\ i=1 & \qquad w_0 x_0 + w_1 x_1 +\ldots+ w_n x_n &=& (b^2-a^2)/2 \\ i=2 & \qquad w_0 x_0^2 + w_1 x_1^2+\ldots+ w_n x_n^2 &=& (b^3-a^3)/3\\ \cdots \\ i=n & \qquad w_0 x_0^n + w_1 x_1^n +\ldots+ w_n x_n^n &=& (b^{n+1}-a^{n+1})/(n+1). \end{aligned}
Il sistema lineare si può riscrivere in forma matriciale compatta \boldsymbol{V}\boldsymbol{w}=\boldsymbol{q} dove
\boldsymbol{V}= \begin{pmatrix} 1 & 1 & \cdots & 1 \\ x_0 & x_1 & \cdots & x_n \\ x_0^2 & x_1^2 & \cdots & x_n^2 \\ \vdots & & & \vdots \\ x_0^n & x_1^n & \cdots & x_n^n \\ \end{pmatrix},\qquad \boldsymbol{q}= \begin{pmatrix} b-a \\ (b^2-a^2)/2 \\ (b^3-a^3)/3\\ \vdots \\ (b^{n+1}-a^{n+1})/(n+1) \end{pmatrix},
nelle incognite w_0,w_1,\ldots,w_n.
La matrice \boldsymbol{V} è una matrice di Vandermonde, il cui determinante
V( x_0 , x_1 ,\ldots, x_n ) = \prod_{i>j}^n ( x_i - x_j )
è sicuramente diverso da zero se i nodi sono scelti a due a due distinti, cioè vale x_i\neq x_j per i\neq j. La non singolarità della matrice \boldsymbol{V} garantisce l’esistenza di una ed una sola soluzione formata da una n+1-upla di pesi w_0,w_1,\ldots,w_n.
Stima dell’errore di integrazione
In questo paragrafo riportiamo per completezza l’enunciato di un teorema assai importante, che permette di stimare l’ordine di accuratezza per le formule di Newton-Cotes 23.
Il teorema si riassume usualmente con l’affermazione che le formule di Newton-Cotes di grado pari guadagnano un ordine di accuratezza. Per esempio, la formula di Simpson, pur essendo costruibile come un formula interpolatoria che integra esattamente tutti i polinomi di grado fino a due, ha in realtà ordine quattro, e non tre come si potrebbe erroneamente pensare.
Questa proprietà dipende essenzialmente dal fatto che l’integrale di \omega_{n+1}(x), che compare nell’espressione dell’errore di fronte al termine n+1 dello sviluppo di Taylor, si annulla sull’intervallo [a,b] per ragioni di simmetria. L’errore è dominato dal termine successivo dello sviluppo di Taylor che coinvolge la derivata n+2-esima.
L’enunciato del teorema mostra la forma generale dell’errore di integrazione per formule interpolatorie alla Newton-Cotes che interpolano esattamente polinomi di ordine fino ad n sia pari che dispari. Ricordiamo tuttavia che nella pratica non si usano formule di Newton-Cotes con n grandi, maggiori di 7-8, per ragioni di stabilità numerica.
La dimostrazione si può trovare per esempio in {cite}Bini:1992.
Teorema 1 (Stima errore quadratura) Sia
\mathcal{E}(f;a,b)= \int_a^b\,f(x)\,\mathrm{d}x-\mathcal{I}(f;a,b)
l’errore di integrazione che si ha applicando una formula di quadratura di Newton-Cotes con n+1 nodi
\mathcal{I}(f;a,b)=\sum_{k=0}^{n} w_k \,f( x_k )
per il calcolo dell’integrale della funzione f(x) su [a\,b]:
\int_a^b\,f(x)\,\mathrm{d}x.
Allora valgono i seguenti risultati:
sia f\in\mathcal{C}^{n+2}(a,b), con n pari; allora esiste un \bar{x}\in (a,b) tale che
\mathcal{E}(f;a,b) = \frac{1}{(n+2)!}f^{n+2}(\bar{x}) \int_a^b\, x \,\omega_{n+1}(x)\,\mathrm{d}x
sia f\in\mathcal{C}^{n+1}(a,b), con n dispari; allora esiste un \bar{x}\in (a,b) tale che
\mathcal{E}(f;a,b) = \frac{1}{(n+1)!}f^{n+1}(\bar{x}) \int_a^b\,\omega_{n+1}( x )\,\mathrm{d}x
dove \omega_{n+1}(x) = ( x - x_0 )( x - x_1 )\cdots( x - x_n ).
Esempio di stima dell’errore della fomula di quadratura
Nel caso dello studio di una singola formula di quadratura al posto del teorema generale Teorema 1 si possono fare dell stime ad hoc. Consideriamo come esempio la regola di quadratura che approssima l’integrale in un intervallo usando il rettangolo con altezza il valore della funzione all’estremo a sinistra dell’intervallo:
\int_{a}^b f(x) \mathrm{d}x \approx f(a)(b-a)
Lemma 1 (Errore regola dei rettangoli) L’errore di integrazione con la regola dei rettangoli a sinistra è ilseguente
\mathcal{E}(f,a,b)=\dfrac{(b-a)^2}{2}f'(\omega),
Dimostrazione. Se introduciamo la primitiva di f(x) indicandola con F(x), possiamo definire la funzione ausiliaria
G(t) := F(a+t)-F(a)- tf(a)-\dfrac{t^2}{h^2}E, \qquad h = b-a
La funzione G si annulla per t=0 e per t=h, infatti:
\begin{aligned} G(0) &= F(a)-F(a)-0-0 = 0\\ G(h) &= F(a+h)-F(a)-hf(a)-\dfrac{h^2}{h^2}\mathcal{E}(f,a,b) \\ &= \underbrace{\int_a^bf(x)\mathrm{d}x -hf(a)}_{\textrm{per definizione}=\mathcal{E}(f,a,b)}-\mathcal{E}(f,a,b) = 0. \end{aligned}
Allora per il teorema di Rolle (o del valor medio) esiste un punto \alpha compreso in 0<\alpha<h tale che G'(\alpha)=0. Possiamo ripetere il procedimento precedente applicandolo a G':
G'(t)=f(a+t)-f(a)-\dfrac{2t}{h^2}\mathcal{E}(f,a,b).
Anche questa volta G'(0)=0, infatti G'(0) = f(a)-f(a) = 0, e per il teorema di Rolle esite \beta tale che 0<\beta<\alpha<h per cui G''(\beta)=0. La derivata seconda di G è
G''(t)=f'(a+t)-\dfrac{2}{h^2}\mathcal{E}(f,a,b),
dunque
G''(\beta)=f'(a+\beta)-\dfrac{2}{h^2}\mathcal{E}(f,a,b)=0 \Rightarrow \mathcal{E}(f,a,b) = \dfrac{h^2f'(a+\beta)}{2}.
e ponendo \omega=a+\beta e ricordando che h=b-a si conclude la dimostrazione.
In maniera molto simile si possono calcolare le stime dell’errore per la regola dei trapezi (interpolazione lineare)
\int_{a}^b f(x) \mathrm{d}x \approx \dfrac{f(a)+f(b)}{2}(b-a)
e la regola di Simpson (interpolazione quadratica)
\int_{a}^b f(x) \mathrm{d}x \approx \dfrac{f(a)+4f((a+b)/2)+f(b)}{6}(b-a)
Integrazione composita
Nella applicazione delle formule di quadratura è conveniente calcolare l’integrale approssimato usando la formula di quadratura su dei sottointervalli. Definendo
H = \dfrac{b-a}{N},\qquad x_0 = a, \quad x_1=a+H, \quad \ldots\quad x_N = a+NH = b
possiamo spezzare l’integrale sui sottintrevalli [x_{k-1},x_k] come segue:
\begin{aligned} \int_a^b f(x)\,\mathrm{d}x &= \sum_{k=1}^N \int_{x_{k-1}}^{x_k} f(x)\,\mathrm{d}x \\ &= \sum_{k=1}^N \Big( \mathcal{I}(f;x_{k-1},x_k)+\mathcal{E}(f;x_{k-1},x_k)\Big) \\ &= \underbrace{\sum_{k=1}^N\mathcal{I}(f;x_{k-1},x_k)}_{\textrm{formula composita}} +\underbrace{\sum_{k=1}^N\mathcal{E}(f;x_{k-1},x_k)}_{\textrm{errore globale}} \end{aligned}
ancora meglio è definire definire nodi e pesi con stima dell’errore in un intervallo di riferimento standard facendo un cambio di variabile. L’intervallo di riferimento è l’intervallo [-1,1]. Facendo il cambio di variabile:
\begin{aligned} x &= \frac{x_{k-1}+x_k}{2}+z\frac{H}{2},\\ \mathrm{d}x &= \frac{H}{2}\mathrm{d}z,\\ g_k(z) &= f\left(\frac{x_{k-1}+x_k}{2}+z\frac{H}{2}\right) \end{aligned}
e quindi
\int_a^b f(x)\,\mathrm{d}x = \sum_{k=1}^N \int_{x_{k-1}}^{x_k} f(x)\,\mathrm{d}x = \frac{H}{2}\sum_{k=1}^N \int_{-1}^1 g_k(z)\,\mathrm{d}z
La formula di quadratura per l’intervallo [-1,1] (grazie al teorema Teorema 1 avrà la forma
\int_{-1}^1 g(z)\mathrm{d}z = \sum_{j=1}^m w_j g(z_j) + \underbrace{\color{blue}C_q\,g^{(q)}(\omega)}_{\textrm{errore}},\qquad \omega\in(-1,1)
dove C_q è una costante che dipende dalla formula di quadratura, e quindi
\begin{aligned} \int_a^b f(x)\,\mathrm{d}x &= \frac{H}{2}\sum_{k=1}^N \int_{-1}^1 g_k(z)\,\mathrm{d}z \\ &= \underbrace{\frac{H}{2}\sum_{k=1}^N \left(\sum_{j=1}^m w_j g_k(z_j)\right) }_{\textrm{formula composita}} +\underbrace{\frac{H}{2} \sum_{k=1}^N{\color{blue}C_q\,g_k^{(q)}(\omega_k)}}_{\textrm{errore globale}} \\ &= \sum_{k=1}^N \left(\sum_{j=1}^m \left(\frac{H}{2}w_j\right) f\left(\frac{x_{k-1}+x_k}{2}+z_j\frac{H}{2}\right)\right) +\frac{NH}{2}\left(\dfrac{1}{N}\sum_{k=1}^N {\color{blue}C_q\,g_k^{(q)}(\omega_k)}\right) \end{aligned}
osservando
\begin{aligned} g_k^{(q)}(\omega_k) &= \left(\dfrac{\mathrm{d}}{\mathrm{d}x}\right)^q f\left(\frac{x_{k-1}+x_k}{2}+z\frac{H}{2}\right) \\ &= \frac{H^q}{2^q} f^{(q)}\left(\frac{x_{k-1}+x_k}{2}+z\frac{H}{2}\right) \end{aligned}
e definendo i nodi z_{k,j} = \frac{x_{k-1}+x_k}{2}+z_j\frac{H}{2}
e i pesi
W_j=\frac{H}{2}w_j
otteniamo la formula di quadratura composita
\begin{aligned} \int_a^b f(x)\,\mathrm{d}x &= \sum_{k=1}^N \sum_{j=1}^m W_j f(z_{k,j}) +C_q\frac{b-a}{2}\frac{H^q}{2^q} \left( \dfrac{1}{N}\sum_{k=1}^N f^{(q)}\left(\frac{x_{k-1}+x_k}{2}+\omega_j\frac{H}{2}\right) \right) \\ &= \underbrace{\sum_{k=1}^N \sum_{j=1}^m W_j f(z_{k,j})}_{\textrm{formula composita}} + \underbrace{C_q(b-a)\frac{H^q}{2^{q+1}} f^{(q)}(\omega)}_{\textrm{errore globale}} ,\qquad \omega\in(a,b) \end{aligned}
Stima a priori degli intervalli
Quando si usa una formula di quadratura (composita) è importante dare una stima a priori del numero di intervalli in cui suddividere l’intervallo [a,b] di intergrazione in modo da mantenere l’errore sotto una tolleranza assegnata. Cioè dato \epsilon tolleranza assegnata vogliamo calcolare N in modo che
\left\lvert\int_a^b f(x)\,\mathrm{d}x- \underbrace{\sum_{k=1}^N \sum_{j=1}^m W_j f(z_{k,j})}_{\textrm{formula composita}}\right\rvert = \left\lvert\underbrace{C_q(b-a)\frac{H^q}{2^{q+1}} f^{(q)}(\omega)}_{\textrm{errore globale}}\right\rvert \leq \epsilon
data la costante
M_q = \max_{x\in[a,b]} \left\lvert f^{(q)}(x)\right\rvert
otteniamo
\left\lvert C_q(b-a)\frac{H^q}{2^{q+1}} f^{(q)}(\omega)\right\rvert \leq \left\lvert C_q\right\rvert(b-a)\frac{H^q}{2^{q+1}} \left\lvert f^{(q)}(\omega)\right\rvert \leq \left\lvert C_q\right\rvert(b-a)\frac{H^q}{2^{q+1}} M_q
scegliendo H in modo che soddisfi
\left\lvert C_q\right\rvert(b-a)\frac{H^q}{2^{q+1}} M_q \leq \epsilon
allora possiamo garantire che l’errore globale è sotto la tolleranza \epsilon. Dalla identità NH=(b-a) possiamo ricavare una diseguaglianza per N:
\left\lvert C_q\right\rvert(b-a)\frac{\left(\frac{b-a}{N}\right)^q}{2^{q+1}} M_q \leq \epsilon, \qquad\Rightarrow\qquad \left\lvert C_q\right\rvert\frac{(b-a)^{q+1}}{\epsilon\,2^{q+1}} M_q \leq N^q
da cui
per cui il più piccolo intero N che soddisfa Equazione 1 può essere usato nella formula di quadratura composita per calcolare l’integrale con la tolleranza assegnata.
Trapezi e Simpson composite
La formula dei trapezi composita si ottiene approssimando l’integrale con polinomi lineari a tratti e prende la forma
\begin{aligned} \int_a^b f(x)\,\mathrm{d}x &= \sum_{k=1}^N \frac{H}{2}(f(x_{k-1})+f(x_{k}))-\dfrac{b-a}{12}H^2 f''(\omega) \\ &= \underbrace{\frac{H}{2}(f(a)+f(b))+H\sum_{k=1}^{N-1}f(x_{k})}_{ \textrm{formula composita}}- \underbrace{\dfrac{b-a}{12}H^2 f''(\omega)}_{ \textrm{errore} } \end{aligned}
Mentre la formula composita di Simpson ottenuta approssimando l’integrale con polinomi quadratici a tratti (parabole) prende la forma
\int_a^b f(x)\,\mathrm{d}x = \sum_{k=1}^N \frac{H}{6}(f(x_{k-1})+4f(x_{k-1/2})+f(x_{k}))-\dfrac{b-a}{2880}H^4 f''''(\omega)
dove x_{k-1/2} = (x_{k-1}+x_k)/2. Nel caso di Simpson si possono definire i nodi della formula composita come segue:
n = 2N, \quad h = \dfrac{b-a}{n} = \dfrac{b-a}{2N} = \dfrac{H}{2}, \qquad \tilde{x}_k = a+kh, \quad k=0,1,\ldots,n
con queste definizioni le formula di Simpson diventa
\int_a^b f(x)\,\mathrm{d}x = \sum_{k=1,3,5,\ldots}^N \frac{h}{3}(f(\tilde{x}_{k-1})+4f(\tilde{x}_k)+f(\tilde{x}_{k+1}))-\dfrac{b-a}{180}h^4 f''''(\omega)
che di solito viene scritto come segue
\begin{aligned} \int_a^b f(x)\,\mathrm{d}x &= \frac{h}{3}\Big( f(a)+4f(\tilde{x}_1)+2f(\tilde{x}_2)+4f(\tilde{x}_3)+2f(\tilde{x}_4)+\cdots+\\ & \qquad\qquad+2f(\tilde{x}_{n-2})+4f(\tilde{x}_{n-1})+4f(b) \Big) -\dfrac{b-a}{180}h^4 f''''(\omega) \end{aligned}
dove si riconosce il pattern dei pesi [1,4,2,4,2,\ldots,2,4,1] della formula quadratura.
Note
la numerazione degli esponenti parte da zero, essendo x^0=1 il primo monomio di ogni polinomio di qualsiasi ordine.↩︎