Calcolo numerico/Formule di Newton-Cotes

Formula dei trapezi

Consideriamo una funzione positiva crescente. Consideriamo la retta passante per $(a,f(a))$ e $(b,f(b))$ e calcoliamo l'area del trapezio sotteso, ponendo $h=b-a$. Il trapezio ha area

$h/2*(f(a)+f(b))*(b-a)=I_{n}(f)$

Osserviamo che il polinomio interpolante $f$ nei nodi $x_{1}=a$ e $x_{2}=b$ è:

$f(a)*l_{1}(x)+f(b)*l_{2}(x)=f(a)*{\frac {x-a}{b-a}}+f(b)*{\frac {x-b}{a-b}}$

quindi

$\int _{a}^{b}f(x)\,dx$

si approssima con

$I_{n}=\int _{a}^{b}f(a)*{\frac {x-a}{b-a}}+f(b)*{\frac {x-b}{a-b}}\,dx=$

$I_{n}=\int _{a}^{b}f(a)*{\frac {x-a}{b-a}}-f(b)*{\frac {x-b}{b-a}}\,dx=$

$I_{n}={\frac {1}{b-a}}*\int _{a}^{b}f(a)*(x-a)-f(b)*(x-b)\,dx=$

$I_{n}={\frac {1}{b-a}}*\int _{a}^{b}x(f(a)-f(b))-af(a)+bf(b)\,dx=$

$I_{n}={\frac {1}{b-a}}*[x^{2}/2(f(a)-f(b))-af(a)x+bf(b)x]_{a}^{b}=$

$I_{n}={\frac {1}{b-a}}*[b^{2}/2f(a)-a^{2}/2f(a)-b^{2}/2f(b)+a^{2}/2f(b)-abf(a)+a^{2}f(a)+b^{2}f(b)-abf(b)]=$

$I_{n}={\frac {1}{b-a}}*[b^{2}/2f(a)+a^{2}/2f(a)+b^{2}/2f(b)+a^{2}/2f(b)-abf(a)-abf(b)]=$

$I_{n}={\frac {1}{2(b-a)}}*[b^{2}f(a)+a^{2}f(a)+b^{2}f(b)+a^{2}f(b)-2abf(a)-2abf(b)]=$

$I_{n}={\frac {1}{2(b-a)}}*[(b-a)^{2}f(a)+(b-a)^{2}f(b)]$

$I_{n}=1/2[(b-a)f(a)+(b-a)f(b)]=(b-a)/2*(f(b)+f(a))\,\quad {\hbox{formula dei trapezi}}$

Supponiamo che $f\in C^{2}([a,b])$, allora

$f(x)-p_{1}(x)={\frac {f^{(2)}(\psi )}{2!}}*w_{2}(x),\quad \psi =\psi (x)$

(deriva dall'espressione dell'errore analitico per l'interpolazione di polinomi).

Quindi si ha

$E_{an}=I(f)-I_{n}(f)=\int _{a}^{b}f(x)-p_{1}(x)\,dx$

$=1/2*\int _{a}^{b}f^{(2)}(\psi (x))*W_{2}(x)\,dx$

Nel caso della formula dei trapezi, $w_{2}(x)=(x-a)*(x-b)$ ha segno negativo sull'intervallo $[a,b]$. Possiamo quindi applicare il teorema del valor medio.

Teorema 6.1

Supponiamo di avere $G\colon [a,b]\to \mathbb {R}$ continua, $\phi \colon [a,b]\to \mathbb {R}$ integrabile, che non cambia segno in $[a,b]$. Allora esiste $\psi \in [a,b]$ tale che

$\int _{a}^{b}g(x)\phi (x)\,dx=g(\psi )*\int _{a}^{b}\phi (x)\,dx$

Allora, riconsiderando la formula dell'errore, siccome la funzione peso non cambia segno, si applica il teorema del valor medio. Esisterà $\psi$ indipendente da $x$ tale che

$E_{an}=1/2*f^{(2)}(\psi )*\int _{a}^{b}w_{2}(x)\,dx$

$E_{an}=1/2*f^{(2)}(\psi )*\int _{a}^{b}(x-a)(x-b)\,dx$

$E_{an}=1/2*f^{(2)}(\psi )*\int _{a}^{b}x^{2}-(a+b)x+ab\,dx$

$E_{an}=1/2*f^{(2)}(\psi )*[x^{3}/3-(a+b)x^{2}/2+abx]_{a}^{b}=$

$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]=$

$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]=$

$E_{an}=1/2*f^{(2)}(\psi )*[a^{3}/6-ba^{2}/2+ab^{2}/2-b^{3}/6]=$

$E_{an}={\frac {1}{12}}*f^{(2)}(\psi )*[a^{3}-3ba^{2}+3ab^{2}-b^{3}]=$

$E_{an}={\frac {1}{12}}*f^{(2)}(\psi )*(b-a)^{3}$

$E_{an}={\frac {H^{3}}{12}}*f^{(2)}(\psi )$

Formula dei trapezi composita

Se l'intervallo è grande, l'errore potrebbe essere elevato, quindi si pensa di passare alla formula composita. Suddividiamo quindi $[a,b]$ in $m$ sottointervalli, consideriamo la partizione

$a=x_{0}<x_{1}<\dots <x_{m}=b$

Poniamo $H=(b-a)/m$ (ipotesi semplificativa: intervalli di uguale ampiezza). Si calcola l'area del trapezio su ogni intervallino, e sommando i vari contributi otteniamo:

$I_{k,m}(f)=H/2*\sum _{k=0}^{m-1}[f(x_{k})+f(x_{k+1})]$

e siccome i nodi interni compaiono due volte possiamo riscriverlo come

$I_{k,m}(f)=H*(f(a)/2+\sum _{i=1}^{m-1}f(x_{i})+f(b)/2)$

$H$ è la distanza tra due nodi successivi, e più la partizione è fine, più l'approssimazione è buona.

Sommando i vari contributi dell'errore analitico, si ha:

$E_{1,m}(f)=-\sum _{k=0}^{m-1}{\frac {f^{(2)}(\psi _{k})}{12}}*H^{3}$

(come dimostrato prima $H^{3}/6$ è l'integrale del polinomio nodale) e moltiplicando e dividendo per $m$:

$=-m{\frac {H^{3}}{12}}*\sum _{k=0}^{m-1}f^{(2)}(\psi (k))/m$

e per il teorema del valor medio discreto, assumendo che $f^{(2)}(x)$ sia continuo,

$=-H^{3}/12*m*f^{(2)}(\psi ),\,\psi \in [a,b]$

$=-H^{2}/12*H*m*f^{(2)}(\psi )$

e siccome $Hm=h$, otteniamo il risultato cercato.

$E_{k,m}(f)=-h/12*H^{2}*f^{(2)}(\psi )$

Formula di Simpson

Consideriamo la funzione $f\colon [a,b]\to \mathbb {R}$, sia $x_{1}$ il punto medio dell'intervallo. Considerando il polinomio di grado 2 interpolante nei tre nodi equispaziati, si ha

$I_{2}(f)=H/3*(f(a)+4f((a+b)/2)+f(b))$

con $h=(b-a)/2$ (passo di equispaziatura dei nodi).

$I_{2}(f)=\int _{a}^{b}p_{2}(x)\,dx$

$=\int _{a}^{b}\sum _{i=0}^{2}f(x_{i})*L_{i}(x)\,dx$

$=\sum _{i}f(x_{i})*\int _{a}^{b}L_{i}(x)\,dx$

Considerando i polinomi di Lagrange di grado 2 otteniamo i pesi $w_{1}=w_{3}=H/3$ e $w_{2}=4H/3$.

Modalità di scrittura dei pesi

Osservazione 6.1

Consideriamo un esempio di formula di quadratura con nodi equispaziati. I pesi della formula di quadratura si possono sempre scrivere come $H\alpha _{i}$ dove gli $\alpha _{i}$ non dipendono dall'intervallo $[a,b]$.

Possiamo scrivere i nodi come

$x_{i}=a+bh,\,h=(b-a)/n$

Il polinomio generico di grado $n$ si scrive come

$\sum _{i=0}^{n}y_{i}*L_{i}(x)$

dove

$L_{i}(x)=\prod _{j\neq i,j=0}^{n}{\frac {x-x_{j}}{x_{i}-x_{j}}}$

Consideriamo un cambio di variabili.

$x=\phi (t)=x_{0}+th$

tale che $\phi (0)=a,\phi (n)=b$.

L'espressione del polinomio di Lagrange diventa

$\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)$

Tenendo conto che $dx=hdt$, il peso si può quindi riscrivere come

$w_{i}=\int _{a}^{b}L_{i}(x)\,dx=\int _{0}^{n}L_{i}(t)hdt$

e definisco

$\alpha _{i}=\int _{0}^{n}L_{i}(t)\,dt$

in modo che $w_{i}=\alpha (t)h$ con $h=(b-a)/n$ nel caso di nodi equispaziati.

Inoltre $\alpha _{u}=\alpha _{n-u}$, cioè i pesi sono simmetrici. Questo deriva dal fatto che anche i nodi equispaziati sull'intervallo $[a,b]$ sono simmetrici.

Metodo dei coefficienti indeterminati

Fissiamo $x_{0},x_{1},\dots ,x_{k}$ nodi distinti, e imponiamo che la formula di quadratura abbia grado di precisione $k$.

Inoltre poniamo $f(x)=x^{j}$, allora per la formula di quadratura deve valere:

$\sum _{i=0}^{m}w_{i}f(x_{i})=\int _{a}^{b}f(x)\,dx$

che si può riscrivere come

$\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$

Otteniamo il sistema lineare

${\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}}$

dove il termine noto ha componenti

$b_{j}=(b^{j+1}-a^{j+1})/(j+1)$

La matrice a primo membro è la matrice di Vandermonde, e se i nodi sono distinti $\det A\neq 0$ esiste unica la soluzione.

Esempio 6.1 (metodo di Simpson)

Consideriamo un intervallo $[0,1]$ con nodi $x_{0}=0,\,x_{1}=1/2,\,x_{2}=1$. La distanza $h$ fra i nodi è $1/2$.

Il sistema corrispondente è:

${\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}}$

Risolviamo il sistema:

${\begin{cases}w_{1}+w_{2}+w_{3}=1\\1/2w_{2}+w_{3}=1/2\\1/4w_{2}+w_{3}=1/3\end{cases}}$

${\begin{cases}w_{1}+w_{2}+w_{3}=1\\1/2w_{2}-1/2=1/4w_{2}-1/3\end{cases}}$

${\begin{cases}w_{1}+w_{2}+w_{3}=1\\1/4w_{2}=1/6\end{cases}}$

$w_{2}=2/3,\,w_{3}=-1/2*2/3+1/2=1/6$

$w_{1}=1-w_{2}-w_{3}=1-2/3-1/6=1/6$

Siccome $h=1/2$, si ha

$\alpha _{0}=\alpha _{2}=1/3,\,\alpha _{1}=4/3$

e si ritrova la formula di Simpson.

Grado di precisione

Teorema 6.2

Supponiamo $f\in C^{4}([a,b])$, allora

$E_{2}(f)=-h^{5}/90*f^{(4)}(\xi ),\quad h=(b-a)/2,\quad \xi \in [a,b]$

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 $P_{2,n}$ come polinomio interpolante, allora il grado di precisione $k$ è $2n+1$.

Formula di Simpson composita

Definiamo $H=(b-a)/m$, e $x_{k}=a+kH/2$, $k=0,\dots ,2m$. Segue che

$w_{0}=1/6,w_{1}=4/6,w_{2}=1/6,w_{3}=4/6,w_{4}=1/6\dots$

$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})]$

(i punti medi vengono contati una volta, gli altri punti due volte).

$E_{2,m}(f)=-(b-a)/180*H^{4}/2*f^{(4)}(\xi ),\quad \xi \in [a,b]$

Per $H\to 0$, l'errore decresce.

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

Formula del punto medio

Consideriamo la funzione $f\colon [a,b]\to \mathbb {R}$, consideriamo il punto medio $(a+b)/2$ e approssimiamo la funzione con la retta per quel punto. Otteniamo quindi:

1. formula semplice: $I_{0,f}=H*f((a+b)/2)$, dove $H=b-a$;
2. errore nella formula semplice: $E_{0}(f)=H^{3}/3f^{(2)}(\xi )$, con gradi di precisione $k=1$;
3. formula composita: $I_{0,n}(f)=H*\sum _{k=0}^{n-1}f(x_{k})$;
4. errore nella formula composita: $E_{0,m}(f)=(b-a)/24H^{2}*f^{(2)}(\xi )$, dove $H=(b-a)/m$, e $f\in C^{2}([a,b])$.

Consideriamo $n$ generico.

$I_{n}(f)=\int _{a}^{b}p_{n}(x)\,dx$

dove $p_{n}$ è il polinomio interpolante di grado $n$, quindi

$I_{n}=\sum _{i}\int _{a}^{b}L_{i}(x)\,dx*f(x_{i})$

e $w=H\alpha (i)$.

Osservazione 6.2

Per grado del polinomio interpolante fino a 7, i coefficienti sono tutti positivi, mentre per $n>8$ i pesi possono anche essere negativi (si può quindi avere la possibilità dell'errore di cancellazione).

Vale il seguente teorema.

Teorema 6.3

La formula di quadratura $I_{n}(f)=\sum _{i=0}^{n}w_{i}f(x_{i})$ con grado di precisione almeno $n$, coincide con la formula di quadratura interpolatoria che usa come nodi di interpolazione $x_{0},x_{1},\dots ,x_{n}$.

Dimostrazione

$I_{n}(f)=\sum _{i=0}^{n}w_{i}f(x_{i})$

è la formula interpolatoria, con $p_{n}$ polinomio interpolante di grado $n$, segue che il grado di precisione $k$ è almeno $n$. Sia $\pi (f)$ il polinomio interpolante di grado $n$ di $f$. Allora esiste unica la soluzione del sistema con matrice di Vandermonde, e quindi le due formule di quadratura coincidono.

Stime dell'errore a priori

Definizione

L'ordine di infinitesimo di una formula di quadratura rispetto a $h$ è il massimo $p\in \mathbb {N}$ tale che

$|I_{n}(f)-I(f)|=o(h^{p})$

Teorema 6.4

Sia $I_{n}(f)$ la formula di quadratura interpolatoria di Newton-Coves (nodi equispaziati e polinomio interpolante di grado $n$), $h=(b-a)/m$ il passo di equispaziatura. Se il grado del polinomio interpolante è pari a $f\in C_{n+2}([a,b])$ allora si dimostra che

$E_{n}(f)=k_{m}/(n+2)!*h^{n+3}f^{(n+2)}(\xi ),\xi \in [a,b]$

(grado di precisione $n+1$ e ordine di infinitesimo $n+3$) con

$K_{m}=\int _{0}^{m}w_{m+1}(t)\,dt$

con $w_{m+1}$ polinomio nodale, e questa quantità è negativa.

Se invece $n$ è dispari $f\in C_{n+1}$ si ha

$E_{n,f}=K_{m}/(n+1)!*H^{n+2}f^{(n+1)}(\xi )$

(grado di precisione $n$ e ordine di infinitesimo $n+2$) con

$K_{m}=\int _{0}^{m}w_{m+1}(t)\,dt$

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

Teorema 6.5

Per le formule composite, segue che, se $n$ è pari,

$E_{n,m}(f)=(b-a)/(n+2)!*{\frac {k_{m}}{n^{n+3}}}H^{n+2}f^{(n+2)}(\xi )$

e il grado di precisione è $n+1$ e l'ordine di infinitesimo è $n+2$, mentre per $n$ dispari,

$E_{n,m}(f)=(b-a)/(n+1)!{\frac {k_{m}}{n^{n+2}}}*H^{n+1}f^{(n+1)}(\xi )$

e il grado di precisione è $n$ e l'ordine di infinitesimo è $n+1$.

In ogni caso

$E_{n,m}(f)\to 0\iff H\to 0$

quindi non ci sono problemi di convergenza.

Raddoppio del numero dei sottointervalli

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

1. Al primo passo, $m=1$, i nodi sono $a$ e $b$, e poniamo $H_{1}=b-a$.

2. Per $m=2$, consideriamo il punto medio $x_{1}$ e chiamiamo $H_{2}=b-x_{1}$.

3. Al passo successivo, $m=4$ (raddoppiamo il numero dei sottointervalli) e dividiamo i due sottointervalli nuovamente a metà, in modo che

$a<x_{2}<x_{1}<x_{3}<b$

dove $x_{2},x_{3}$ sono punti medi dei due sosttointervallini.

4. Ad ogni livello vengono introdotti nuovi nodi.

Si ha quindi

$I_{1,1}(f)=H_{1}/2(f(a)+f(b))$

$I_{1,2}(f)=H_{2}/2(f(a)+2f(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})$

$I_{1,2m}(f)=I_{1m}(f)/2+H_{2m}*\sum _{j=1}^{m}f(x_{j})$

dove gli $x_{j}$ sono i nuovi nodi.

Esercizio 6.1

Nel caso della formula di Simpson, raddoppiano i sottointervalli

$I_{1,m}(f)=a_{2m}+b_{2m}+c_{2m}.$

Infatti al passo 1,

$a_{1}=H_{1}/3(f(a)+f(b))$

$b_{1}=0$

$c_{1}=4H_{1}/3f((a+b)/2)=4H_{1}/3*f(x_{1})$

Al passo 2, poniamo $H_{2}=H_{1}/2$:

$E_{2}=H_{2}/3(f(a)+f(b)+2f(x_{1})+4f(x_{2})+4f(x_{3}))$

dove $x_{2},x_{3}$ sono i punti medi di $[a,x_{1}]$ e $[x_{1},b]$.

$E_{2}=H_{2}/3(f(a)+f(b))+2/3H_{2}f(x_{1})+4H_{2}/3(f(x_{2})+f(x_{3})))$

$E_{2}=H_{1}/6(f(a)+f(b))+1/3H_{1}f(x_{1})+4H_{2}/3(f(x_{2})+f(x_{3})))$

$E_{2}=a_{1}/2+c_{1}/4+4H_{2}/3*\sum _{i=2}^{3}f(x_{i})$

cioè

$a_{2}=a_{1}/2,\,b_{2}=c_{1}/4,\,c_{2}=4H_{2}/3*\sum _{i=2}^{3}f(x_{i})$

Al passo $n+1$, avendo $2m=3+2^{n}$ nodi, si ha:

$a_{2m}=a_{m}/2$

$b_{2m}=b_{m}/2+c_{m}/4$

$c_{2m}=4H_{2m}/3\sum f({\hbox{x nuovi}})$