Formule di Newton-Coves

==Formula dei trapezi==
Considero una funzione positiva crescente, considero la retta passante per <math>(a,f(a))</math> e <math>(b,f(b))</math> e calcolo l'area del trapezio sotteso, ponendo <math>h=b-a</math>.
Il trapezio ha area
<math display="block">h/2*(f(a)+f(b))*(b-a) = I_n(f)</math>

Osservo che il polinomio interpolante <math>f</math> nei nodi <math>x_1=a</math> e <math>x_2=b</math> è:
<math display="block">f(a)*l_1(x)+f(b)*l_2(x) = f(a)*\frac{x-a}{b-a}+f(b)*\frac{x-b}{a-b}</math>
quindi
<math display="block">\int_a^b f(x) \, dx</math>
si approssima con
<math display="block">I_n = \int_a^b f(a)*\frac{x-a}{b-a}+f(b)*\frac{x-b}{a-b} \, dx=</math><math display="block">I_n = \int_a^b f(a)*\frac{x-a}{b-a}-f(b)*\frac{x-b}{b-a} \, dx=</math><math display="block">I_n = \frac{1}{b-a}*\int_a^b f(a)*(x-a)-f(b)*(x-b) \, dx=</math><math display="block">I_n = \frac{1}{b-a}*\int_a^b x(f(a)-f(b))-a f(a)+b f(b) \, dx=</math><math display="block">I_n = \frac{1}{b-a}*[x^2/2 (f(a)-f(b))-a f(a) x+b f(b) x]_a^b=</math><math display="block">I_n = \frac{1}{b-a}*[b^2/2 f(a)-a^2/2 f(a)-b^2/2 f(b)+a^2/2 f(b)-ab f(a)+a^2 f(a)+b^2 f(b)-ab f(b)]=</math><math display="block">I_n = \frac{1}{b-a}*[b^2/2 f(a)+a^2/2 f(a)+b^2/2 f(b)+a^2/2 f(b)-ab f(a)-ab f(b)]=</math><math display="block">I_n = \frac{1}{2(b-a)}*[b^2 f(a)+a^2 f(a)+b^2 f(b)+a^2 f(b)-2ab f(a)-2ab f(b)]=</math><math display="block">I_n = \frac{1}{2(b-a)}*[(b-a)^2 f(a)+(b-a)^2 f(b)]</math><math display="block">I_n = 1/2 [(b-a) f(a)+(b-a) f(b)]  = (b-a)/2*(f(b)+f(a)) \, \quad \hbox{formula dei trapezi}</math>

Supponiamo che <math>f \in C^2([a,b])</math>, allora
<math display="block">f(x)-p_1(x) = \frac{f^{(2)}(\psi)}{2!}*w_{2}(x), \quad \psi = \psi(x)</math>
(deriva dall'espressione dell'errore analitico per l'interpolazione di polinomi)


Quindi si ha
<math display="block">E_{an} = I(f)-I_n(f) = \int_a^b f(x)-p_1(x) \, dx</math><math display="block">= 1/2*\int_a^b f^{(2)}(\psi(x))*W_2(x) \, dx</math>

Nel caso della formula dei trapezi, <math>w_2(x) = (x-a)*(x-b)</math> ha segno negativo sull'intervallo <math>[a,b]</math>. Posso quindi applicare il teorema del valor medio.

{{InizioTeorema|titolo=|number=6.1|anchor=Teorema6_1}}
Supponiamo di avere <math>G \colon [a,b] \to \mathbb R</math> continua, <math>\phi \colon [a,b] \to \mathbb R</math> integrabile, che non cambia segno in <math>[a,b]</math>. Allora esiste <math>\psi \in[a,b]</math> tale che
<math display="block">\int_a^b g(x) \phi(x) \, dx = g(\psi)*\int_a^b \phi(x) \, dx</math>
{{FineTeorema}}


Allora, riconsiderando la formula dell'errore, siccome la funzione peso non cambia segno si applica il teorema del valor medio, esisterà <math>\psi</math> indipendente da x tae che
<math display="block">E_{an} = 1/2*f^{(2)}(\psi)*\int_a^b w_2(x) \, dx</math><math display="block">E_{an} = 1/2*f^{(2)}(\psi)*\int_a^b (x-a)(x-b) \, dx</math><math display="block">E_{an} = 1/2*f^{(2)}(\psi)*\int_a^b x^2-(a+b)x+ab \, dx</math><math display="block">E_{an} = 1/2*f^{(2)}(\psi)*[x^3/3-(a+b)x^2/2+ab x]_a^b =</math><math display="block">E_{an} = 1/2*f^{(2)}(\psi)*[b^3/3-a^3/3+(a+b) a^2/2-(a+b) b^2/2+ab^2-a^2 b] =</math><math display="block">E_{an} = 1/2*f^{(2)}(\psi)*[b^3/3-a^3/3+a^3/2+ba^2/2-ab^2/2-b^3/2+ab^2-a^2 b] =</math><math display="block">E_{an} = 1/2*f^{(2)}(\psi)*[a^3/6-ba^2/2+ab^2/2-b^3/6] =</math><math display="block">E_{an} = \frac{1}{12}*f^{(2)}(\psi)*[a^3-3ba^2+3ab^2-b^3] =</math><math display="block">E_{an} = \frac{1}{12}*f^{(2)}(\psi)*(b-a)^3</math><math display="block">E_{an} = \frac{H^3}{12}*f^{(2)}(\psi)</math>

==Formula dei trapezi composita==
Se l'intervallo è grande, l'errore potrebbe essere elevato, quindi si pensa di passare alla formula composita.


Suddivido quindi <math>[a,b]</math> in <math>m</math> sottointervalli, considero la partizione
<math display="block">a = x_0 < x_1 < \dots < x_m=b</math>
Pongo <math>H=(b-a)/m</math> (ipotesi semplificativa: intervalli di uguale ampiezza). Si calcola l'area del trapezio su ogni intervallino, e sommando i vari contributi ottengo:
<math display="block">I_{k,m}(f) = H/2*\sum_{k=0}^{m-1} [f(x_k)+f(x_{k+1})]</math>
e siccome i nodi interni compaiono due volte posso riscriverlo come
<math display="block">I_{k,m}(f) = H*(f(a)/2+\sum_{i=1}^{m-1} f(x_i)+f(b)/2)</math>

<math>H</math> è la distanza tra due nodi successivi, e più la partizione è fine, più l'approssimazione è buona.


Sommando i vari contributi dell'errore analitico, si ha:
<math display="block">E_{1,m}(f) = -\sum_{k=0}^{m-1} \frac{f^{(2)}(\psi_k)}{12}*H^3</math>
(come dimostrato prima <math>H^3/6</math> è l'integrale del polinomio nodale)
e moltiplicando e dividendo per m:
<math display="block">= -m \frac{H^3}{12}*\sum_{k=0}^{m-1} f^^{(2)}(\psi(k))/m</math>
e per il teorema del valor medio discreto, assumendo che <math>f^{(2)}(x)</math> sia continuo,
<math display="block">= -H^3/12*m*f^{(2)}(\psi), \, \psi \in [a,b]</math><math display="block">= -H^2/12*H*m*f^{(2)}(\psi)</math>
e siccome <math>Hm =h</math> ottengo il risultato cercato.
<math display="block">E_{k,m}(f) = -h/12*H^2*f^{(2)}(\psi)</math>

==Formula di Simpson==
Considero la funzione <math>f \colon [a,b] \to \mathbb R</math>, sia <math>x_1</math> il punto medio dell'intervallo. Considerando il polinomio di grado 2 interpolante nei tre nodi
equispaziati, si ha
<math display="block">I_2(f) = H/3*(f(a)+4 f((a+b)/2)+f(b))</math>
con <math>h=(b-a)/2</math> (passo di equispaziatura dei nodi).
<math display="block">I_2(f) = \int_a^b p_2(x) \, dx</math><math display="block">= \int_a^b \sum_{i=0}^2 f(x_i)*L_i(x) \, dx</math><math display="block">= \sum_i f(x_i)*\int_a^b L_i(x) \, dx</math>
Considerando i polinomi di Lagrange di grado 2 ottengo i pesi <math>w_1=w_3=H/3</math> e <math>w_2=4H/3</math>.

==Modalità di scrittura dei pesi==
{{InizioOsservazione|titolo=|number=6.1|anchor=Osservazione6_1}}
Considero un esempio di formula di quadratura con nodi equispaziati. I pesi della formula di quadratura si possono sempre scrivere come <math>H \alpha_i</math>
dove gli <math>\alpha_i</math> non dipendono dall'intervallo <math>[a,b]</math>.
{{FineOsservazione}}


Posso scrivere i nodi come
<math display="block">x_i = a+b h, \, h=(b-a)/n</math>
Il polinomio generico di grado <math>n</math> si scrive come
<math display="block">\sum_{i=0}^n y_i*L_i(x)</math>
dove
<math display="block">L_i(x) = \prod_{j \neq i, j=0}^n \frac{x-x_j}{x_i-x_j}</math>
Considero un cambio di variabili.
<math display="block">x = \phi(t) = x_0+th</math>
tale che <math>\phi(0) = a, \phi(n) = b</math>.


L'espressione del polinomio di Lagrange diventa
<math display="block">\prod_{j \neq i} (x-x_j)/(x_i-x_j) = \prod_{j \neq i} \frac{x_0+th-(x_0+hj)}{x_0+ih-(x_0+jh)} = \prod_{j \neq i} (t-j)/(i-j)</math>
Tenendo conto che <math>dx = h dt</math>, il peso si può quindi riscrivere come
<math display="block">w_i = \int_a^b L_i(x) \, dx = \int_0^n L_i(t) h dt</math>
e definisco
<math display="block">\alpha_i = \int_0^n L_i(t) \, dt</math>
in modo che <math>w_i = \alpha(t) h</math> con <math>h=(b-a)/n</math> nel caso di nodi equispaziati.



Inoltre <math>\alpha_u = \alpha_{n-u}</math>, cioè i pesi sono simmetrici. Questo deriva dal fatto che anche i nodi equispaziati sull'intervallo <math>[a,b]</math> sono simmetrici.

==Metodo dei coefficienti indeterminati==
Fisso <math>x_0, x_1, \dots, x_k</math> nodi distinti, e impongo che la formula di quadratura abbia grado di precisione <math>k</math>.


Inoltre pongo <math>f(x) = x^j</math>, allora per la formula di quadratura deve valere:
<math display="block">\sum_{i=0}^m w_i f(x_i) = \int_a^b f(x) \, dx</math>
che si può riscrivere come
<math display="block">\sum_{i=0}^m w_i x_i^j = \int_a^b x^j \, dx = \frac{b^{j+1}-a^{j+1}}{j+1}, \quad j=0,\dots,k</math>
Ottengo il sistema lineare
<math display="block">\begin{array}{cccc}1 & 1 & \dots & 1  \\x_0 & x_1 & \dots & x_k \\x_0^2 & x_1^2 & \dots & x_k^2 \\\dots & \dots & \dots & \dots \\x_0^k & x_1^k & \dots & x_n^k\end{array}  \times\begin{array}{c}w_0 \\ w_1 \\ \dots \\ w_k\end{array}= \boldsymbol{b}</math>
dove il termine noto ha componenti
<math display="block">b_j = (b^{j+1}-a^{j+1})/(j+1)</math>
La matrice a primo membro è la matrice di Vandermonde, e se i nodi sono distinti <math>\det A \neq 0</math> quindi esiste unica la soluzione.


{{InizioEsempio|titolo= metodo di Simpson|number=6.1|anchor=Esempio6_1}}
Considero un intervallo <math>[0,1]</math> con nodi <math>x_0=0, \, x_1=1/2, \, x_2=1</math>. La distanza <math>h</math> fra i nodi è <math>1/2</math>.


Il sistema corrispondente è:
<math display="block">\begin{array}{ccc} 1 & 1 & 1 \\ 0 & 1/2 & 1 \\ 0 & 1/4 & 1 \end{array}\times\begin{array}{c}w_0 \\ w_1 \\ w_2 \end{array} =\begin{array}{c} 1 \\ 1/2 \\ 1/3 \end{array}</math>
Risolvo il sistema:
<math display="block">\begin{cases}w_1+w_2+w_3 = 1 \\ 1/2 w_2+w_3 = 1/2  \\ 1/4 w_2+w_3 = 1/3 \end{cases}</math><math display="block">\begin{cases}w_1+w_2+w_3 = 1 \\ 1/2 w_2-1/2=1/4 w_2-1/3 \end{cases}</math><math display="block">\begin{cases}w_1+w_2+w_3 = 1 \\ 1/4 w_2=1/6 \end{cases}</math><math display="block">w_2 =2/3, \, w_3 = -1/2*2/3+1/2 = 1/6</math><math display="block">w_1 = 1-w_2-w_3 = 1-2/3-1/6 = 1/6</math>
Siccome <math>h=1/2</math>, si ha
<math display="block">\alpha_0 = \alpha_2 = 1/3, \, \alpha_1 = 4/3</math>
e si ritrova la formula di Simpson.
{{FineEsempio}}

==Grado di precisione==
{{InizioTeorema|titolo=|number=6.2|anchor=Teorema6_2}}
Supponiamo <math>f \in C^4([a,b])</math>, allora
<math display="block">E_2(f) = -h^5/90*f^{(4)}(\xi), \quad h=(b-a)/2, \quad \xi \in [a,b]</math>
{{FineTeorema}}

Questa formula ha grado di precisione 3 (tutti i polinomi di grado 3 hanno derivata quarta nulla, e quindi errore nullo)

Tutte le formule interpolatorie su nodi equispaziati quando il polinomio interpolante ha grado pari hanno un grado di precisione in più di quanto ci si aspetta, Questo deriva dalla simmetria nei nodi.

Vale in generale che dato <math>P_{2,n}</math> come polinomio interpolante, allora il grado di precisione <math>k</math> è <math>2n+1</math>.

==Formula di Simpson composita==
Definisco <math>H = (b-a)/m</math>, e <math>x_k = a+k H/2</math>, <math>k=0,\dots, 2m</math>.
Segue che
<math display="block">w_0 = 1/6, w_1 = 4/6, w_2 = 1/6, w_3 = 4/6, w_4 = 1/6 \dots</math><math display="block">P_{2,m} = H/6*[f(x_0)+2*\sum_{s=1}^{m-1} f(x_{2s})+4 \sum_{s=0}^{m} f(x_{2s-1})+f(x_m)]</math>
(i punti medi vengono contati una volta, gli altri punti due volte)
<math display="block">E_{2,m}(f) = -(b-a)/180*H^4/2*f^{(4)}(\xi), \quad \xi \in [a,b]</math>
Per <math>H \to 0</math>, l'errore decresce.



La formula di quadratura di Simpson e dei trapezi sono dette formule chiuse perché usano gli estremi dell'intervallo <math>[a,b]</math>. Altre formule vengono dette
aperte perché tutti i nodi di quadratura sono interni all'intervallo.

==Formula del punto medio==
Considero la funzione <math>f \colon [a,b] \to \mathbb R</math>, considero il punto medio <math>(a+b)/2</math> e approssimo la funzione con la retta per quel punto. Ottengo quindi:

#''formula semplice'':  <math>I_{0,f} = H*f((a+b)/2)</math>, dove <math>H=b-a</math>.
#''Errore nella formula semplice'': <math>E_0(f) = H^3/3 f^{(2)}(\xi)</math>, con gradi di precisione <math>k=1</math>.
#''Formula composita'': <math>I_{0,n}(f) = H*\sum_{k=0}^{n-1} f(x_k)</math>
#''Errore nella formula composita'': <math>E_{0,m}(f) = (b-a)/24 H^2*f^{(2)}(\xi)</math>, dove <math>H=(b-a)/m</math>, e <math>f \in C^2([a,b])</math>.


Considero <math>n</math> generico.
<math display="block">I_n(f) = \int_a^b p_n(x) \, dx</math>
dove <math>p_n</math> è il polinomio interpolante di grado <math>n</math>, quindi
<math display="block">I_n = \sum_i \int_a^b L_i(x) \, dx*f(x_i)</math>
e <math>w = H \alpha(i)</math>.


{{InizioOsservazione|titolo=|number=6.2|anchor=Osservazione6_2}}
Per grado del polinomio interpolante fino a 7, i coefficienti sono tutti positivi, mentre
per <math>n > 8</math> i pesi possono anche essere negativi (si può quindi avere la possibilità dell'errore di cancellazione).
{{FineOsservazione}}


Vale il seguente teorema:

{{InizioTeorema|titolo=|number=6.3|anchor=Teorema6_3}}
La formula di quadratura <math>I_n(f)=\sum_{i=0}^n w_i f(x_i)</math> con grado di precisione almeno <math>n</math>, coincide con la formula di quadratura interpolatoria che
usa come nodi di interpolazione <math>x_0,x_1, \dots, x_n</math>.
{{FineTeorema}}

{{InizioDimostrazione}}
<math display="block">I_n(f) = \sum_{i=0}^n w_i f(x_i)</math>
è la formula interpolatoria, con <math>p_n</math> polinomio interpolante di grado <math>n</math>, segue che il grado di precisione <math>k</math> è almeno <math>n</math>. Sia <math>\pi(f)</math> il polinomio
interpolante di grado <math>n</math> di <math>f</math>. Allora esiste unica la soluzione del sistema con matrice di Vandermonde, e quindi le due formule di quadratura coincidono.
{{FineDimostrazione}}

==Stime dell'errore a priori==
{{InizioDefinizione|titolo=|number=6.3|anchor=Definizione6_3}}
L'''ordine di infinitsimo'' di una formula di quadratura rispetto a <math>h</math> è  il massimo <math>p \in \mathbb N</math> tale che
<math display="block">|I_n(f)-I(f)| = o(h^p)</math>
{{FineDefinizione}}



{{InizioTeorema|titolo=|number=6.4|anchor=Teorema6_4}}
Sia <math>I_n(f)</math> la formula di quadratura interpolatoria di Newton-Coves (nodi equispaziati e polinomio interpolante di grado <math>n</math>), <math>h = (b-a)/m</math> il passo
di equispaziatura. Se il grado del polinomio interpolante è
pari e <math>f \in C_{n+2}([a,b])</math> allora si dimostra che
<math display="block">E_n(f) = k_m/(n+2)!*h^{n+3} f^{(n+2)}(\xi), \xi \in [a,b]</math>
(grado di precisione <math>n+1</math> e ordine di infinitesimo <math>n+3</math>)
con
<math display="block">K_m =\int_0^m w_{m+1}(t) \, dt</math>
con <math>w_{m+1}</math> polinomio nodale, e questa quantità è negativa.


Se invece <math>n</math> è dispari e <math>f \in C_{n+1}</math> si ha
<math display="block">E_{n,f} = K_m/(n+1)!*H^{n+2} f^{(n+1)}(\xi)</math>
(grado di precisione <math>n</math> e ordine di infinitesimo <math>n+2</math>)
con
<math display="block">K_m = \int_0^m w_{m+1}(t) \, dt</math>
{{FineTeorema}}


Ad esempio, nel caso della formula di Simpson (tre nodi, quindi <math>m=3</math>) dove il polinomio interpolante ha grado 2, il grado di precisione è <math>n+1=3</math>, e l'ordine di
infinitesimo è <math>n+3=5</math>.


{{InizioTeorema|titolo=|number=6.5|anchor=Teorema6_5}}
Per le formule composite, segue che, se <math>n</math> è pari,
<math display="block">E_{n,m}(f) = (b-a)/(n+2)!*\frac{k_m}{n^{n+3}} H^{n+2} f^{(n+2)}(\xi)</math>
e il grado di precisione è <math>n+1</math> e l'ordine di infinitesimo è <math>n+2</math>, mentre per <math>n</math> dispari,
<math display="block">E_{n,m}(f) = (b-a)/(n+1)! \frac{k_m}{n^{n+2}}*H^{n+1} f^{(n+1)}(\xi)</math>
e il grado di precisione è <math>n</math> e l'ordine di infinitesimo è <math>n+1</math>.
In ogni caso
<math display="block">E_{n,m}(f) \to 0 \iff H \to 0</math>
quindi non ci sono problemi di convergenza.
{{FineTeorema}}

==Raddoppio del numero dei sottointervalli==
Consideriamo cosa succede se si raddoppia il numero dei sottointervalli. Nel caso della formula dei trapezi:

#Al primo passo, <math>m=1</math>, i nodi sono <math>a</math> e <math>b</math>, e pongo <math>H_1=b-a</math>.
#Per <math>m=2</math>, considero il punto medio <math>x_1</math> e chiamo <math>H_2 = b-x_1</math>.
#Al passo successivo, <math>m=4</math> (raddoppio il numero dei sottointervalli) e divido i due sottointervalli nuovamente a metà, in modo che<math display="block">a < x_2 < x_1 < x_3 < b</math>dove <math>x_2,x_3</math> sono punti medi dei due sosttointervallini.
#Ad ogni livello vengono introdotti nuovi nodi.

Si ha quindi
<math display="block">I_{1,1}(f) = H_1/2 (f(a)+f(b))</math><math display="block">I_{1,2}(f) = H_2/2 (f(a)+2 f(x_1)+f(b)) = H_2/2 (f(a)+f(b))+H_2 f(x_1) = I_{11}(f)/2+H_2 f(x_1)</math><math display="block">I_{1,2m}(f) = I_{1m}(f)/2 +H_{2m}*\sum_{j=1}^m f(x_j)</math>
dove gli <math>x_j</math> sono i nuovi nodi.


{{InizioEsercizio|titolo=|number=6.1|anchor=Esercizio6_1}}
Nel caso della formula di Simpson, raddoppiano i sottointervalli
<math display="block">I_{1,m}(f) = a_{2m}+b_{2m}+c_{2m}.</math>
Infatti,

#al passo 1,<math display="block">a_1 = H_1/3 (f(a)+f(b))</math><math display="block">b_1 = 0</math><math display="block">c_1 =  4H_1/3 f((a+b)/2) = 4H_1/3*f(x_1)</math>
#Al passo 2, pongo <math>H_2 = H_1/2</math>:<math display="block">E_2 = H_2/3(f(a)+f(b)+2f(x_1)+4f(x_2)+4f(x_3))</math>dove <math>x_2, x_3</math> sono i punti medi di <math>[a,x_1]</math> e <math>[x_1,b]</math>.<math display="block">E_2 = H_2/3(f(a)+f(b))+2/3 H_2 f(x_1)+4H_2/3(f(x_2)+f(x_3)))</math><math display="block">E_2 = H_1/6(f(a)+f(b))+1/3 H_1 f(x_1)+4H_2/3(f(x_2)+f(x_3)))</math><math display="block">E_2 = a_1/2+c_1/4+4 H_2/3*\sum_{i=2}^3 f(x_i)</math>cioè<math display="block">a_2 = a_1/2, \, b_2 = c_1/4, \, c_2 = 4H_2/3*\sum_{i=2}^3 f(x_i)</math>
#al passo <math>n+1</math>, avendo <math>2m=3+2^n</math> nodi, si ha:<math display="block">a_{2m} = a_m/2</math><math display="block">b_{2m} = b_m/2+c_m/4</math><math display="block">c_{2m} = 4 H_{2m}/3 \sum f(\hbox{x nuovi})</math>
{{FineEsercizio}}