Integrazione Approssimata

Formule di quadratura

Sia $f: [a, b] \to \mathbb{R}$ . Vogliamo approssimare

 $S[f] := \int_a^b f(x) \mathrm dx$

§ Definizione. Considereremo formule di quadratura della forma

 $S_{n+1}[f] = \sum_{i=0}^n w_i f(x_i)$

con $x_i \in [a, b]$ detti nodi e $w_i$ detti pesi (più precisamente una formula di quadratura è quindi identificata da $\{(x_i, w_i)\}_{i=0}^n$ ). La quantità

 $r_{n+1} := S[f] - S_{n+1}[f]$

che rappresenta l’errore di troncamento è detto resto.

§ Definizione. Il grado di precisione di una formula $\sum_{i=0}^n w_i f(x_i)$ è l’intero $k$ tale che per ogni funzione:

$f(x) = x^j, j=0, \dots, k \implies r_{n+1} = 0$
$f(x) = x^{k+1} \implies r_{n+1} \neq 0$

Per linearità dell’operatore integrale, se una formula ha grado di precisione $k$ dà un risultato esatto per tutti i polinomi di grado $\le k$ .

Se imponiamo che la formula abbia grado di precisione $k$ per $x_i$ fissati, otteniamo un sistema lineare per i $w_i$ , per ogni $j=0, \dots, k$ :

 $\sum_{i=0}^n w_i x_i^j = \int_a^b x^j \mathrm dx = \left[ \frac{1}{j+1} x^{j+1} \right]_a^b$

Se gli $x_i$ sono tutti distinti e $k = n$ otteniamo una matrice di Vandermonde:

 $\begin{bmatrix} & \vdots & \\ \dots & x_i^j & \dots \\ & \vdots & \end{bmatrix} \begin{bmatrix} w_0 \\ \vdots \\ w_n \end{bmatrix} = \begin{bmatrix} \frac{b^{j+1} - a^{j+1}}{j+1} \\ \vdots \end{bmatrix}$

Osservazione. In particolare un modo di costruire formule di quadratura è approssimare $f(x)$ con il suo polinomio di interpolazione di Lagrange $p_n(x)$ con $\text{deg} \le n$ e approssimare $S[f]$ con $\int_a^b p_n(x) \mathrm dx$ .

Consideriamo il polinomio dell’interpolazione di Lagrange:

 $\begin{aligned} \pi_n(x) &:= (x-x_0) \cdots (x-x_n) \\ p_n(x) &:= \sum_{i=0}^n L_i(x) f(x_i) \end{aligned}$

con

 $L_i(x) := \prod_{j=0, j \neq i}^n \frac{x-x_j}{x_i-x_j} = \frac{\pi_n (x)}{(x-x_i)\pi_n'(x_i)}$

Allora:

 $\begin{gathered} \begin{aligned} S_{n+1}[f] &= \int_a^b p_n(x) \mathrm dx \\ &= \int_a^b \sum_{i=0}^n L_i(x) f(x_i) \mathrm dx \\ &= \sum_{i=0}^n \left( \int_a^b L_i(x) \mathrm dx \right) f(x_i) \\ &= \sum_{i=0}^n w_i f(x_i) \end{aligned} \\ \implies w_i = \int_a^b L_i(x) \mathrm dx = \frac{1}{\pi_n'(x_i)} \int_a^b \frac{\pi_n(x)}{x-x_i} \mathrm dx \end{gathered}$

Queste sono dette formule di quadratura interpolatorie.

Se i nodi sono disposti in modo simmetrico allora anche i $w_i$ lo saranno.
Se $f \in C^{n+1}[a, b]$ e gli $x_i$ sono tutti distinti, allora l’errore $r(x) = f(x) - p_n(x)$ verifica:
```
 $r(x) = \pi_n(x) \frac{f^{(n+1)}(\xi)}{(n+1)!}$ 
```
da cui si ottiene per il resto della formula di quadratura:
```
 $r_{n+1} = \frac{1}{(n+1)!} \int_a^b \pi_n(x) f^{(n+1)}(\xi) \mathrm dx$ 
```

§ Teorema. Una formula di integrazione $S_{n+1}[f]$ è interpolatoria $\iff$ grado di precisione $\ge n$ .

$\boxed{\Rightarrow}$ Se $f(x)$ è un polinomio con $\text{deg} \le n$ , allora $p_n(x)$ (polinomio di interpolazione di grado $\le n$ ) coincide con $f(x) \implies S[f] = S_{n+1}[f]$ e la formula ha grado di precisione $\ge n$ .
$\boxed{\Leftarrow}$ Se la formula ha grado di precisione $\ge n$ . Se poniamo $f(x) = L_k(x)$ che hanno grado $n$ , vale $S[f] = S_{n+1}[f]$ , quindi:
```
 $\int_a^b L_k(x) \mathrm dx = \sum_{i=0}^n w_i L_k(x_i) = w_k$ 
```
poiché $L_k(x_i) = \delta_{ki}$ .

Proprietà di convergenza

§ Teorema. Sia $f \in C[a, b]$ e sia $\{ S_{n+1}(f) \}_n$ una successione di formule di quadratura interpolatorie con

 $S_{n+1}[f] = \sum_{i=0}^n w_i^{(n)} f(x_i^{(n)})$

tali che $\exists H$ costante per cui

 $\forall n \quad \sum_i |w_i^{(n)}| < H$

allora $\lim_n r_{n+1} = 0$ .

Dimostrazione. $\forall \varepsilon > 0$ per il teorema di Weierstrass, esiste un polinomio $p$ tale che

 $\forall x \in [a, b] \quad |f(x) - p(x)| < \varepsilon$

dunque

 $\begin{aligned} r_{n+1} &= \int_a^b f - S_{n+1}[f] = \\ &= \int_a^b (f-p) + \int_a^b p - \sum_i w_i^{(n)} f(x_i^{(n)}) \end{aligned}$

Se $\text{deg } p = m$ , allora $\forall n \ge m$ , $S_{n+1}[f]$ ha grado di precisione $\ge m$ per cui

 $\int_a^b p = \sum_i w_i^{(n)} p(x_i^{(n)})$

da cui

 $\begin{aligned} \implies r_{n+1} &= \int_a^b (f-p) + \sum_i w_i^{(n)} [ p(x_i^{(n)}) - f(x_i^{(n)}) ] \\ |r_{n+1}| &\le \varepsilon \left( \int_a^b 1 \mathrm dx + \sum_i |w_i^{(n)}| \right) \\ &\le \varepsilon \cdot (b - a + H) \end{aligned}$

$\square$

Osservazione. Esistono quindi successioni $\{ x_i^{(n)}, w_i^{(n)} \}_{i, n}$ di nodi distinti e pesi tali che $\forall f \in C[a, b]$

 $\lim_n \sum_{i=0}^n w_i^{(n)} f(x_i^{(n)}) = \int_a^b f(x) \mathrm dx$

(Però non ci sono successioni di un solo indice $\{ x_i, w_i \}_i$ ).

Tutte le formule di quadratura interpolatorie a coefficienti positivi sono convergenti.

§ Teorema. I coefficienti di una formula interpolatoria $S_{n+1}(f)$ con grado di precisione $\ge 2n$ sono tutti positivi.

Siano $x_0, \dots, x_n$ i nodi di una formula di quadratura interpolatoria $S_{n+1}[f]$ . Consideriamo per $i=0, \dots, n$ i polinomi:

 $q_i(x) = \prod_{\substack{r=0 \\ r \neq i}}^n (x - x_r)$

Hanno tutti grado $n$ e poiché $S_{n+1}[f]$ ha grado di precisione $\ge 2n$ , essa integra esattamente i polinomi $q_i(x)^2$ , cioè:

 $\sum_{j=0}^n w_j q_i(x_j)^2 = \int_a^b q_i(x)^2 \mathrm dx$

Dunque per ogni $i$ :

 $w_i \prod_{r \neq i} (x_i - x_r)^2 = \int_a^b q_i(x) \mathrm dx > 0$

poiché il prodotto a sinistra è $> 0$ (essendo un quadrato di termini non nulli) e l’integrale a destra è di una funzione strettamente positiva quasi ovunque. Segue $w_i > 0$ . $\square$

Formule notevoli

Newton-Cotes: si scelgono gli $x_i$ equispaziati su $[a, b]$ . Hanno grado di precisione $n$ o $n+1$ . Per $n \ge 8$ i coefficienti non sono tutti dello stesso segno.
Gaussiane: si scelgono gli $x_i$ in modo da massimizzare il grado di precisione che risulta essere $2n+1$ ed hanno coefficienti sempre positivi.

Formule di Newton-Cotes

§ Definizione (Newton-Cotes). Siano $h = (b - a)/n$ e $x_i = a + ih$ per $i = 0, \dots, n$ su $[a, b]$ . Una formula di quadratura $S_{n+1}[f]$ interpolatoria su questi $n+1$ nodi equidistanti è detta di Newton-Cotes.

Esempio (Newton-Cotes per $n=1$ ). In questo caso $h = b - a$ e i nodi sono $x_0 = a, x_1 = b$ . Si ottiene la regola del trapezio:

 $S_2[f] = \frac{h}{2} (f(x_0) + f(x_1))$

Calcoliamo i pesi:

 $w_0 = \int_{x_0}^{x_1} \frac{x - x_1}{x_0 - x_1} \mathrm dx = \dots = \frac{h}{2} (= w_1)$

Esempio (Newton-Cotes per $n=2$ ). In questo caso $h = (b - a)/2$ e i nodi sono:

 $x_0 = a, \quad x_1 = \frac{a+b}{2}, \quad x_2 = b$

Si ottiene la regola di Simpson:

 $S_3[f] = \frac{h}{3} (f(x_0) + 4f(x_1) + f(x_2))$

Calcoliamo i pesi:

 $w_0 = \int_{x_0}^{x_2} \frac{(x - x_1)(x - x_2)}{(x_0 - x_1)(x_0 - x_2)} \mathrm dx = \dots = \frac{h}{3} (= w_2)$

 $w_1 = \int_{x_0}^{x_2} \frac{(x - x_0)(x - x_2)}{(x_1 - x_0)(x_1 - x_2)} \mathrm dx = \dots = \frac{4}{3}h$

§ Teorema (Errore di Newton-Cotes). Sia $S_{n+1}[f]$ una formula di Newton-Cotes.

Se $n$ è pari e $f \in C^{n+2}[a, b]$ , allora $\exists \xi \in (a, b)$ tale che: $r_{n+1} = \frac{f^{(n+2)}(\xi)}{(n+2)!} \int_a^b x \pi_n(x) \mathrm dx$
Se $n$ è dispari e $f \in C^{n+1}[a, b]$ , allora $\exists \xi \in (a, b)$ tale che: $r_{n+1} = \frac{f^{(n+1)}(\xi)}{(n+1)!} \int_a^b \pi_n(x) \mathrm dx$

Le formule di Newton-Cotes hanno grado di precisione $n+1$ per $n$ pari e $n$ per $n$ dispari.

Osservazione. Le formule di Newton-Cotes per $n \ge 8$ presentano coefficienti non tutti positivi, il che implica che non soddisfano le ipotesi del teorema di convergenza.

Formule di Newton-Cotes composte

L’idea delle formule composte è quella di dividere l’intervallo $[a, b]$ in $N$ sottointervalli $[z_k, z_{k+1}]$ di uguale ampiezza $H = (b-a)/N$ .

Si ha quindi:

 $\int_a^b f(x) \mathrm dx = \sum_{k=0}^{N-1} \int_{z_k}^{z_{k+1}} f(x) \mathrm dx \approx \sum_{k=0}^{N-1} S_{n+1}^{(k)}[f] := J_{n+1}^{(N)}[f]$

Esempio (Formula dei trapezi composta, $n=1$ ). Dividendo $[a, b]$ in $N$ intervalli con nodi $z_k = a + k \frac{b-a}{N}$ , si ottiene:

 $J_2^{(N)}[f] = \frac{b-a}{2N} \left( f(z_0) + 2 \sum_{k=1}^{N-1} f(z_k) + f(z_N) \right)$

Esempio (Formula di Simpson composta, $n=2$ ). Dividendo $[a, b]$ in $N$ intervalli, si ottiene:

 $J_3^{(N)}[f] = \frac{b-a}{6N} \left( f(z_0) + 2 \sum_{k=1}^{N-1} f(z_k) + 4 \sum_{k=0}^{N-1} f\left(\frac{z_k + z_{k+1}}{2}\right) + f(z_N) \right)$

Osservazione. Per le formule di Newton-Cotes composte, il resto totale si analizza sommando i resti sui singoli sottointervalli. Utilizzando il teorema della media, è possibile esprimere il resto globale in termini della derivata della funzione in un punto intermedio dell’intervallo $[a, b]$ .

Osservazione. Nel caso della formula dei trapezi composta ( $n=1$ ), il resto su un singolo intervallo è:

 $r_2^{(k)} = -\frac{h^3}{12} f''(\xi_k) \quad \text{con } h = z_{k+1} - z_k \text{ e } \xi_k \in (z_k, z_{k+1})$

Per la formula composta su $N$ intervalli si ha:

 $\begin{aligned} R_2^{(N)} &= \sum_{k=0}^{N-1} r_2^{(k)} = -\frac{h^3}{12} \sum_{k=0}^{N-1} f''(\xi_k) \\ &= -\frac{h^3}{12} N f''(\xi) \quad \text{(usando il teorema della media)} \\ &= -\frac{(b-a)h^2}{12} f''(\xi), \quad \xi \in (a, b) \end{aligned}$

Osservazione. Per la formula di Simpson composta ( $n=2$ ):

 $\begin{aligned} R_3^{(N)} &= \sum_{k=0}^{N-1} r_3^{(k)} = -\frac{1}{90} \left(\frac{h}{2}\right)^5 \sum_{k=0}^{N-1} f^{(4)}(\xi_k) \\ &= -\frac{(b-a)h^4}{2880} f^{(4)}(\xi), \quad \xi \in (a, b) \end{aligned}$

Stima dell’errore e estrapolazione

Per il metodo dei trapezi:

 $S[f] - J_2^{(N)}[f] = \frac{\delta_N}{N^2}, \quad S[f] - J_2^{(2N)}[f] = \frac{\delta_{2N}}{4N^2} \approx \frac{\delta_N}{4N^2}$

Studiamo cosa succede raddoppiando il numero di intervalli da $N$ a $2N$ , ottenendo (supponendo $\delta_N \approx \delta_{2N} = \delta$ ):

 $J_2^{(2N)}[f] - J_2^{(N)}[f] \approx \frac{3\delta}{4N^2}$

da cui si ottiene la seguente stima del resto:

 $S[f] - J_2^{(2N)}[f] \approx \frac{J_2^{(2N)}[f] - J_2^{(N)}[f]}{3}$

Dunque si può procedere con successivi raddoppi (riciclando le valutazioni di $f$ già fatte) fino a quando

 $\left| \frac{J_2^{(2N)}[f] - J_2^{(N)}[f]}{3} \right| < \varepsilon$

e estrapolare il valore dell’integrale finale usando l’ultimo resto:

 $J_2^{(2N)}[f] + \frac{J_2^{(2N)}[f] - J_2^{(N)}[f]}{3} = \frac{4 J_2^{(2N)}[f] - J_2^{(N)}[f]}{3}$

e viene detta estrapolazione di Richardson.

Vediamo ora una tecnica che applica iterativamente l’estrapolazione di Richardson, usando il metodo dei trapezi.

§ Teorema (Eulero-Maclaurin). Sia $f \in C^{2r+2}[a, b]$ e sia $x_i = a + ih$ , $i=0, \dots, N$ , $h = (b-a)/N$ e sia

 $T(a:h:b)[f] = \frac{h}{2} \left( f(a) + 2 \sum_{i=1}^{N-1} f(x_i) + f(b) \right)$

Allora

 $T(a:h:b)[f] - \int_a^b f(x) \mathrm dx = \sum_{j=1}^r \frac{B_{2j} h^{2j}}{(2j)!} \left( f^{(2j-1)}(b) - f^{(2j-1)}(a) \right) + R_{2r+2}(f, a, b, h)$

con $R_{2r+2}(f, a, b, h) = O(h^{2r+2})$ e i coefficienti $B_j$ sono i numeri di Bernoulli definiti da

 $B_0 = 1, \quad \sum_{i=0}^{n-1} \binom{n}{i} B_i = 0, \quad n \ge 2$

e vale

 $\begin{aligned} |R_{2r+2}(f, a, b, h)| &\le 2 \frac{|B_{2r+2}| h^{2r+2}}{(2r+2)!} \int_a^b |f^{(2r+2)}(x)| \mathrm dx \\ &\approx 2 \left( \frac{h}{2\pi} \right)^{2r+2} \int_a^b |f^{(2r+2)}(x)| \mathrm dx \end{aligned}$

Per la formula precedente, per $f$ sufficientemente regolare esistono costanti $c_1, \dots, c_r$ tali che:

 $T(a:h:b) - \int_a^b f(x) \mathrm dx = \sum_{j=1}^r c_j h^{2j} + R_{2r+2}(f, a, b, h)$

(non ci interessano le costanti $c_j$ …) Per due passi diversi otteniamo $h_1, h_2$ :

 $\begin{aligned} &\sum_{j=1}^r c_j h_1^{2j} + R_{2r+2}(f, a, b, h_1) \\ \text{e } &\sum_{j=1}^r c_j h_2^{2j} + R_{2r+2}(f, a, b, h_2) \end{aligned}$

Quindi moltiplicando la prima per $h_2^2 / h_1^2$ e sottraendo otteniamo:

 $T(a:h_2:b)[f] - \frac{h_2^2}{h_1^2} T(a:h_1:b)[f] - \left(1-\frac{h_2^2}{h_1^2}\right) S[f] = O(h_2^4)$

Possiamo quindi definire

 $\tilde{T} = \frac{T(a:h_2:b)[f] - h_2^2/h_1^2 T(a:h_1:b)}{1 - h_2^2/h_1^2}$

ed avremo $\tilde{T} - S[f] = O(h_2^4)$ che è superiore.

Posto $h_1 = 2h$ e $h_2 = h$ e $f_i := f(a+ih)$ per $i=0, \dots, N$ vale

 $\begin{aligned} \tilde{T} &= \frac{T(a:h:b)[f] - 1/4 T(a:2h:b)[f]}{1 - 1/4} = \\ &= \frac{h}{3} (f_0 + 4f_1 + 2f_2 + \dots + 2f_{N-2} + 4f_{N-1} + f_N) \end{aligned}$

detto metodo di Cavalieri-Simpson.

§ Definizione (Metodo di Romberg). Il metodo di Romberg consiste nell’applicare questa tecnica con passi successivi:

 $h_1 = \frac{b-a}{n_1}, \quad h_2 = \frac{h_1}{n_2}, \quad \dots, \quad h_q = \frac{h_1}{n_q}$

Con $n_1 < n_2 < \dots < n_q$ interi, poniamo

 $T_{m,1} := T(a:h_m:b)[f] \quad \forall m=1, \dots, q$

 $T_{m, k+1} = T_{m, k} + \frac{T_{m, k} - T_{m-1, k}}{(h_{m-k}/h_m)^2 - 1} \quad \forall k$

che ci permette di calcolare i valori di estrapolazione. Di solito per riciclare il lavoro già fatto si sceglie $h_m / h_{m-1} = 2$ , in questo modo il denominatore $(h_{m-k}/h_m)^2 - 1$ diventa $2^{2k} - 1$ .

 $\begin{bmatrix} T_{1,1} \\ T_{2,1} & T_{2,2} \\ T_{3,1} & T_{3,2} & T_{3,3} \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix}$

Calcolato iterativamente per righe e ci si ferma quando due termini successivi in una riga sono minori di una certa tolleranza ovvero quando $|T_{m, k} - T_{m, k-1}| < \varepsilon$ .

§ Teorema. In questo caso l’errore del metodo di Romberg è

 $T_{m, k} - \int_a^b f = r_k h^{2k} (b-a) f^{(2k)}(\xi)$

per $\xi \in (a, b)$ con

 $r_k = 2^{k(k-1)} |B_{2k}| / (2k)! \quad h = \frac{b-a}{2^{m-1}}$

Superconvergenza metodo trapezi

Osservazione. Il metodo dei trapezi integra esattamente solo i polinomi di grado 1 ma si comporta molto bene con i polinomi trigonometrici.

§ Teorema. Sia $t_m$ un polinomio trigonometrico di grado $m$

 $t_m(x) = \sum_{j=0}^m (a_j \cos(jx) + b_j \sin(jx))$

e si consideri $S[t_m] = \int_0^{2\pi} t_m(x) \mathrm dx$ , allora la formula dei trapezi con $h = 2\pi/N$ e $N > m$ calcola esattamente $S[t_m]$ .

Formule di Fejer e Clenshaw-Curtis

I nodi equispaziati, come abbiamo già visto, funzionano per $n < 8$ , altrimenti non è garantita la convergenza a zero.

Formule in cui i nodi di interpolazione si addensano agli estremi hanno migliori proprietà per $n \gg 0$ .

Formule di Fejer: sono formule aperte e come nodi si usano gli zeri dei polinomi di Chebyshev di I specie o II specie.
Formule di Clenshaw-Curtis: sono formule chiuse ed usano gli zeri dei polinomi di Chebyshev di II specie più gli estremi $\{-1, 1\}$ .

§ Vogliamo approssimare $\int_{-1}^{+1} f(x) \mathrm dx$ . Come nodi prendiamo:

 $\begin{gathered} x_k := \cos \theta_k, \quad \theta_k := \frac{(2k-1)\pi}{2n}, \quad k = 1, \dots, n \\ w_k^{F1} = \frac{2}{n} \left( 1 - 2 \sum_{j=1}^{\lfloor n/2 \rfloor} \frac{\cos(2j\theta_k)}{4j^2-1} \right) \end{gathered}$

§ Le formule di Fejer del II tipo usano:

 $\begin{gathered} x_k := \cos \theta_k, \quad \theta_k = \frac{k\pi}{2n}, \quad k = 1, \dots, n-1 \\ w_k^{F2} := \frac{4 \sin \theta_k}{n} \sum_{j=1}^{\lfloor n/2 \rfloor} \frac{\sin(2j-1)\theta_k}{2j-1} \end{gathered}$

Osservazione. La formula di II specie è la più pratica perché ci permette di riciclare computazione.

§ Le quadrature di Clenshaw-Curtis invece sono formule chiuse su $[-1, 1]$ , con $x_0 = -1, x_n = 1$ .

 $\begin{gathered} x_k = \cos \theta_k, \quad \theta_k = \frac{k\pi}{n}, \quad k = 0, \dots, n \\ w_k^{cc} = \frac{c_k}{n} \left( 1 - \sum_{j=1}^{\lfloor n/2 \rfloor} \frac{b_j}{4j^2 - 1} \cos(2j\theta_k) \right) \end{gathered}$

con

 $b_j = \begin{cases} 1 & j = n/2 \\ 2 & j < n/2 \end{cases}, \quad c_k = \begin{cases} 1 & k = 0, n \\ 2 & \text{altrimenti} \end{cases}$

Osservazione. Il calcolo esplicito dei pesi può essere evitato esprimendo $p_n(x)$ per $f(x)$ nella base dei polinomi di Chebyshev:

 $p_n(x) = \frac{1}{2} c_0 T_0(x) + \sum_{k=1}^{n-1} c_k T_k(x) + \frac{1}{2} c_n T_n(x)$

con

 $c_k = \frac{2}{n} \left( \frac{f(x_0)}{2} + \sum_{r=1}^{n-1} f(x_r) \cos\left(\frac{kr\pi}{n}\right) + \frac{f(x_n)}{2} \right)$

Questa può essere vista come la parte reale di una IDFT (calcolabile usando una FFT) dunque

 $\begin{aligned} S_{n+1}[f] &= \int_{-1}^1 p_n(x) \mathrm dx \\ &= \frac{1}{2} c_0 \mu_0 + \sum_{k=1}^{n-1} c_k \mu_k + \frac{1}{2} c_n \mu_n \end{aligned}$

con $\mu_k$ i momenti dei polinomi di Chebyshev:

 $\mu_k = \int_{-1}^1 T_k(x) \mathrm dx = \begin{cases} 0 & k \text{ dispari} \\[0.5em] \dfrac{2}{1-k^2} & k \text{ pari} \end{cases}$

Formule Gaussiane

Vorremmo ora cercare formule di quadratura interpolatorie con grado di precisione il più alto possibile.

Generalizziamo il problema di approssimazione a:

 $S[f] = \int_a^b f(x) \omega(x) \mathrm dx$

§ Teorema. Sia $n \geq 0$ e $1 \leq k \leq n + 1$ e sia $S_{n+1}[f]$ una formula di integrazione interpolatoria, allora:

$S_{n+1}[f]$ ha grado di precisione $\leq n + k$ $\iff$ per ogni $p(x)$ con $\deg p \lt n$ vale:

 $\int_a^b p(x) \pi_n(x) \omega(x) \mathrm dx = \langle p, \pi_n \rangle = 0$

dove $\omega(x)$ è una funzione peso positiva quasi ovunque su $[a, b]$ .

Dimostrazione.

$\boxed{\Rightarrow}$ Sia $p$ un polinomio con $\deg p \lt k$ . Consideriamo il polinomio $p(x) \pi_n(x)$ , abbiamo anche che $\deg p(x) \lt k \Leftrightarrow \deg p(x) \leq k - 1$ e che $\deg \pi_n(x) = n + 1$ . Quindi $\deg p(x) \pi_n(x) \leq n + k$ e per ipotesi la formula di quadratura lo integra esattamente ovvero

 $\int_a^b p(x) \pi_n(x) \omega(x) \mathrm dx = S_{n+1}[p(x) \pi_n(x)] = \sum_{i=0}^n w_i p(x_i) \underbrace{\pi_n(x_i)}_{=0} = 0$

$\boxed{\Leftarrow}$ Sia $p(x)$ un polinomio con $\deg p(x) \leq n + k$ , facciamo la divisione euclidea di $p(x)$ per $\pi_n(x)$

 $p(x) = q(x) \pi_n(x) + r(x)$

osserviamo che $\deg q(x) \lt k$ e $\deg r(x) \lt n + 1$ . Passiamo ora all’integrale

 $\int_a^b p(x) \omega(x) \mathrm dx \\ = \underbrace{\int_a^b q(x) \pi_n(x) \omega(x) \mathrm dx}_{=0 \text{ per ipotesi}} + \int_a^b r(x) \omega(x) \mathrm dx$

D’altra il termine in $r(x)$ viene integrato esattamente dalla formula di quadratura

 $\int_a^b r(x) \omega(x) \mathrm dx = \sum_{i=0}^n w_i r(x_i)$

ma osserviamo che valutando $p = q \pi_n + r$ nei nodi $x_i$ si ha

 $p(x_i) = q(x_i) \underbrace{\pi_n(x_i)}_{=0} + r(x_i) = r(x_i)$

$\Rightarrow r(x_i) = p(x_i)$ . Infine otteniamo che

 $\int_a^b p(x) \omega(x) \mathrm dx = \sum_{i=0}^n w_i p(x_i) = S_{n+1}[p]$

ovvero $S_{n+1}$ integra esattamente $p(x)$ di grado $\leq n + k$ . $\square$

Questo ci suggerisce come scegliere i nodi per massimizzare il grado di precisione. Si usano gli zeri dell’ $(n+1)$ -esimo polinomio ortogonale $p_{n+1}$ su $[a, b]$ rispetto al peso $\omega(x) \equiv 1$ .

§ Teorema (Formula per i pesi gaussiani). Sia $p_n(x)$ l’n-esimo polinomio ortogonale monico su $[a, b]$ rispetto a $\omega(x)$ , $h_n := \langle p_n, p_n \rangle$ , allora i coefficienti di $S_{n+1}[f]$ sono dati da:

 $w_i := \frac{h_n}{p_{n+1}'(x_i) p_n(x_i)}$

Consideriamo la formula di Christoffel-Darboux con $y = x_i$ :

 $\begin{aligned} (x - x_i) \sum_{j=0}^n \frac{p_j(x) p_j(x_i)}{h_j} &= \frac{1}{h_n} [p_{n+1}(x) p_n(x_i) - \underbrace{p_{n+1}(x_i)}_{=0} p_n(x)] \\ \implies \sum_{j=0}^n \frac{p_j(x) p_j(x_i)}{h_j} &= \frac{p_{n+1}(x) p_n(x_i)}{h_n (x - x_i)} \end{aligned}$

Integrando rispetto a $\omega(x)$ su $[a, b]$ :

 $\int_a^b \sum_{j=0}^n \frac{1}{h_j} p_j(x) p_j(x_i) \cdot \omega(x) \mathrm dx = \sum_{j=0}^n \frac{1}{h_j} p_j(x_i) \int_a^b p_j(x) \omega(x) \mathrm dx = \frac{h_0}{h_0} = 1$

poiché $p_0(x) = p_0(x_i) = 1$ e gli integrali per $i > 0$ sono nulli per ortogonalità. Per l’altro termine osserviamo che:

 $\begin{aligned} \int_a^b \frac{p_{n+1}(x) p_n(x_i)}{h_n (x - x_i)} \omega(x) \mathrm dx & = \\ = \frac{1}{h_n} p_n(x_i) \int_a^b \frac{p_{n+1}(x)}{x - x_i} \omega(x) \mathrm{d}x & = 1 \\ \implies \int_a^b \frac{p_{n+1}(x)}{x-x_i} \omega(x) \mathrm dx & = \frac{h_n}{p_n(x_i)} \end{aligned}$

Ora ricordiamo che per la definizione di polinomio di Lagrange:

 $L_i(x) = \prod_{\substack{j=0 \\ j \neq i}}^n \frac{x - x_j}{x_i - x_j} = \frac{\pi_{n+1}(x)}{(x - x_i) \pi_{n+1}'(x_i)}$

E nel nostro caso $\pi_{n+1}(x) = p_{n+1}(x)$ visto che i nodi sono gli zeri di $p_{n+1}(x)$ . Dalla definizione di pesi per formule interpolatorie:

 $\begin{aligned} w_i &= \int_a^b L_i(x) \omega(x) \mathrm dx \\ &= \frac{1}{p_{n+1}'(x_i)} \int_a^b \frac{p_{n+1}(x)}{x-x_i} \omega(x) \mathrm dx \\ &= \frac{h_n}{p_{n+1}'(x_i) p_n(x_i)} \end{aligned}$

$\square$

Calcolo di nodi e pesi

A parte il caso dei polinomi di Chebyshev, il calcolo dei nodi e dei pesi delle formule gaussiane risulta computazionalmente più efficiente se si sfruttano certe relazioni.

I polinomi ortogonali rispettano la relazione a 3 termini:

 $p_0 = 1, \quad p_{j+1} = (x + B_j) p_j - C_j p_{j-1}$

Introduciamo la matrice di Jacobi $J_{n+1}$ (tridiagonale):

 $J_{n+1} = \begin{bmatrix} -B_0 & 1 & & \\ C_1 & -B_1 & 1 & \\ & \ddots & \ddots & \ddots \\ & & C_n & -B_n \end{bmatrix}, \quad q(x) = \begin{bmatrix} p_0 \\ \vdots \\ p_n \end{bmatrix}$

allora la stessa relazione in forma matriciale diventa:

 $J_{n+1} q(x) = x q(x) - \begin{bmatrix} 0 \\ \vdots \\ 0 \\ p_{n+1}(x) \end{bmatrix}$

Osserviamo che $x_i$ sono gli zeri di $p_{n+1}$ :

 $J_{n+1} q(x_i) = x_i q(x_i)$

ovvero $x_i$ è un autovalore di $J_{n+1}$ e $q(x_i)$ è il corrispondente autovettore con prima componente 1.

E sempre per la formula di Christoffel-Darboux:

 $\sum_{j=0}^n \frac{p_j(x) p_j(x_i)}{h_j} = \frac{1}{h_n} \frac{p_{n+1}(x) p_n(x_i) - p_{n+1}(x_i) p_n(x)}{x - x_i}$

e per $x \to x_i$ si ha:

 $\sum_{j=0}^n \frac{p_j^2(x_i)}{h_j} = \frac{1}{h_n} p_{n+1}'(x_i) p_n(x_i) = \frac{1}{w_i}$