QUADRATURA NUMERICA^{∗}

A. SOMMARIVA^{†}

Conoscenze richieste. Integrale di Riemann. Teorema di Weierstrass. Polinomi di Lagrange. Derivate di ordine superiore. Operatori lineari limitati. Teorema di Weierstrass sulla densità di uno spazio polinomiale rispetto C([a, b]) con [a, b] chiuso. Programmazione in Matlab/Octave.

Conoscenze ottenute. Formule di quadratura. Grado di precisione. Formule interpolatorie. Formule di Newton-Cotes. Regola del trapezio e di Simpson. Formule composte. Errore di alcune formule di quadratura. Formule a pesi positivi. Errore formule di quadratura di Newton-Cotes. Errore formule di quadratura di Gauss. Stabilità di una formula di quadratura. Teorema di Polya-Steklov.

1. Introduzione. Un classico problema dell’analisi numericaè quello di calcolare l’integrale definito di una funzione f in un intervallo avente estremi di integrazione a, b (non necessariamente finiti) cioè

I_{w}(f) := I_{w}(f, a, b) = ∫_{a}^{b} f(x) w(x)dx

dove w è una funzione peso in (a, b) [1, p.206, p.270].

La nostra intenzione è di approssimare I(f) come

I_{w}(f) ≈ Q_{N}(f) := ∑_{i=1}^{N} w_{i}f(x_{i})  (1.1)

I termini w_{i} e x_{i} ∈ [α, β] sono detti rispettivamente pesi e nodi.

2. Formule interpolatorie. Sia (a, b) l’intervallo di integrazione (non necessariamente limitato), x_{1}, …, x_{N} un insieme di N punti a due a due distinti ed f ∈ C([a, b]) una funzione w-integrabile cioè per cui esista finito I_{w}(f). Se l’intervallo è limitato, per il teorema di Weierstrass e l’integrabilità della funzione peso, questo è vero per qualsiasi funzione continua in quanto

|∫_{a}^{b} f(x) w(x)dx| ≤ ∫_{a}^{b} |f(x)| w(x)dx ≤ ∥f∥_{∞}∥w∥_{1} < +∞

Se

p_{N−1}(x) = ∑_{i=1}^{N} f(x_{i})L_{i}(x)

è il polinomio che interpola le coppie (x_{i}, f(x_{i})) con i = 1, …, N, dove al solito L_{i} indica l’i-simo polinomio di Lagrange allora

∗Ultima revisione: 25 novembre 2011.

†Dipartimento di Matematica Pura ed Applicata, Universitá degli Studi di Padova, stanza 419, via Trieste 63, 35121 Padova, Italia (alvise@euler.math.unipd.it). Telefono: +39-049-8271350.

∫_{a}^{b} f(x)w(x)dx ≈ ∫_{a}^{b} p_{N−1}(x)w(x)dx

= ∫_{a}^{b} ∑_{i=1}^{N} f(x_{i})L_{i}(x)w(x)dx

= ∑_{i=1}^{N} (∫_{a}^{b} L_{i}(x) w(x)dx) f(x_{i})  (2.1)

per cui, confrontando con la formula (1.1) abbiamo

w_{i} = ∫_{a}^{b} L_{i}(x) w(x)dx, i = 1, …, N.

Definizione 2.1. Una formula di quadratura

∫_{a}^{b} f(x)w(x)dx ≈ ∑_{i=1}^{N} w_{i}f(x_{i})  (2.2)

per cui

w_{k} = ∫_{a}^{b} L_{i}(x) w(x)dx, k = 1, …, N  (2.3)

si dice interpolatoria.

Definizione 2.2. Una formula

∫_{a}^{b} f(x)w(x)dx ≈ ∑_{i=1}^{M} w_{i}f(x_{i})

ha grado di precisione almeno N se e solo se è esatta per tutti i polinomi f di grado inferiore o uguale a N. Ha inoltre grado di precisione N se e solo se è esatta per ogni polinomio di grado N ed esiste un polinomio di grado N + 1 per cui non lo sia.

Mostriamo ora il seguente

Teorema 2.3. Una formula

∫_{a}^{b} f(x)w(x)dx ≈ ∑_{i=1}^{N} w_{i}f(x_{i})

è interpolatoria se e solo se ha grado di precisione almeno N − 1.

Dimostrazione. Se f = p_{N−1} è un polinomio di grado N −1 ovviamente corrisponde col polinomio interpolante p_{N−1} nei nodi a due a due distinti x_{1}, …, x_{N} e quindi la formula risulta esatta per polinomi di grado inferiore o uguale a N − 1, cioè

I_{w}(p_{N−1}) = ∑_{i=1}^{N} w_{i}f(x_{i}), w_{i} = ∫_{a}^{b} L_{i}(x) w(x)dx, i = 1, …, N.

Viceversa se è esatta per ogni polinomio di grado N − 1 allora lo è in particolare per i polinomi di Lagrange L_{i} ∈ P_{n−1}, il che implica che w_{i} = ∫_{a}^{b} L_{i}(x) w(x)dx e quindi i pesi sono proprio quelli della formula interpolatoria corrispondente nei nodi x_{1}, …, x_{N}.

Figura 3.1. Regola del trapezio e di Cavalieri-Simpson per il calcolo di ∫_{0.5}^{2} sin (x) dx (rispettivamente area in magenta e in azzurro).

3. Formule di Newton-Cotes. Si supponga [a, b] un intervallo chiuso e limitato di ℝ e si ponga w ≡ 1. Per semplicità di notazione, in questo caso porremo I := I_{w} = I_{1}. Il primo esempio di formule interpolatorie che consideriamo sono le regole di tipo Newton-Cotes chiuse (cf. [9, p.336]) che si ottengono integrando l’interpolante di f in nodi equispaziati

x_{i} = a + \frac{(i − 1) (b − a)}{N − 1}, i = 1, …, N.

Alcune classiche regole sono:

1. regola del trapezio

I(f) ≈ S_{1}(f) := S_{1}(f, a, b) := \frac{(b − a) (f(a) + f(b))}{2}

avente grado di precisione 1, cioè esatta per polinomi di grado inferiore o uguale a 1; si può dimostrare (con un po’ di fatica) dal teorema del resto per l’interpolazione polinomiale (cf. [3, p.132]) che l’errore della regola del trapezio [30] è

E_{1}(f) := I(f) − S_{1}(f) = \frac{−h^{3}}{12} f^{(2)}(ξ)

per qualche ξ ∈ (a, b);

2. regola di Cavalieri-Simpson

I(f) ≈ S_{3}(f) := S_{3}(f, a, b) := \frac{b − a}{6} [f(a) + 4f(\frac{a + b}{2}) + f(b)]

avente grado di precisione 3, cioè esatta per polinomi di grado inferiore o uguale a 3; si può dimostrare (non facile!) che l’errore della regola di Cavalieri-Simpson [31] è

E_{3}(f) := I(f) − S_{3}(f) = \frac{−h^{5}}{90} f^{(4)}(ξ), h = \frac{b − a}{2}

per qualche ξ ∈ (a, b);

Nota 3.1. Vediamo calcolando i pesi, che in effetti le due formule sono interpolatorie. Partiamo dalla regola del trapezio. Posti x_{1} = a, x_{2} = b abbiamo che

L_{1}(x) = \frac{x − b}{a − b}, L_{2}(x) = \frac{x − a}{b − a}

Figura 3.2. Isaac Newton (1642-1727), Roger Cotes (1682-1716), Bonaventura Francesco Cavalieri (1598-1647) e Thomas Simpson (1710-1761).

e quindi visto che w ≡ 1 abbiamo

w_{1} = ∫_{a}^{b} L_{1}(x) dx = ∫_{a}^{b} \frac{x − b}{a − b} dx = \frac{1}{a − b} ∫_{a}^{b} (x − b) dx
= \frac{1}{a − b} \frac{(x − b)^{2}}{2} |_{a}^{b} = \frac{1}{a − b} \frac{(x − b)^{2}}{2} |_{a}^{b} = \frac{1}{a − b} \frac{−(a − b)^{2}}{2} = \frac{b − a}{2} (3.1)

e

w_{2} = ∫_{a}^{b} L_{2}(x) dx = ∫_{a}^{b} \frac{x − a}{b − a} dx = \frac{1}{b − a} ∫_{a}^{b} (x − a) dx
= \frac{1}{b − a} \frac{(x − a)^{2}}{2} |_{a}^{b} = \frac{1}{b − a} \frac{(x − a)^{2}}{2} |_{a}^{b} = \frac{1}{b − a} \frac{(b − a)^{2}}{2} = \frac{b − a}{2} (3.2)

Per quanto riguarda la regola di Cavalieri-Simpson i ragionamenti sono analoghi. D’altra parte essendo quelle dei trapezi e Simpson regole rispettivamente aventi 2 e 3 punti con grado 2 e 4, allora sono entrambe interpolatorie.

Per ulteriori dettagli si confronti [1, p.252-258], [9, p.333-336].

Nota 3.2. Qualora le funzioni da integrare non siano sufficientemente derivabili, una stima dell’errore viene fornita dalle formule dell’errore via nucleo di Peano ([1, p.259]). Ricordiamo che per N ≥ 8 le formule di Newton-Cotes chiuse hanno pesi di segno diverso e sono instabili dal punto di vista della propagazione degli errori (cf. [3, p.196]).

4. Formule di Newton-Cotes composte. Si suddivida l’intervallo (chiuso e limitato) [a, b] in N subintervalli T_{j} = [x_{j}, x_{j+1}] tali che x_{j} = a + jh con h = (b − a)/N. Dalle proprietà dell’integrale

∫_{a}^{b} f(x) dx = ∑_{j=0}^{N−1} ∫_{x_{j}}^{x_{j+1}} f(x) dx ≈ ∑_{j=0}^{N−1} S(f, x_{j}, x_{j+1})  (4.1)

dove S è una delle regole di quadratura finora esposte (ad esempio S_{3}(f)). Le formule descritte in (4.1) sono dette composte. Due casi particolari sono

1. formula composta dei trapezi

S_{1}^{(c)} := h [\frac{f(x_{0})}{2} + f(x_{1}) + … + f(x_{N−1}) + \frac{f(x_{N})}{2}]  (4.2)

il cui errore è

E_{1}^{(c)}(f) := I(f) − S_{1}^{(c)}(f) = \frac{−(b − a)}{12}h^{2} f^{(2)}(ξ), h = \frac{(b − a)}{N}

per qualche ξ ∈ (a, b);

2. formula composta di Cavalieri-Simpson fissati il numero N di subintervalli e i punti x_{k} = a + kh/2 dove h = \frac{b−a}{N} sia

I(f) ≈ S_{3}^{(c)}(f) := \frac{h}{6} [f(x_{0}) + 2 ∑_{r=1}^{N−1} f(x_{2r}) + 4 ∑_{s=0}^{N−1} f(x_{2s+1}) + f(x_{2N})];

(4.3)

il cui errore è

E_{3}^{(c)}(f) := I(f) − S_{3}^{(c)}(f) = \frac{−(b − a)}{180} (\frac{h}{2})^{4} f^{(4)}(ξ)

per qualche ξ ∈ (a, b).

5. Formule gaussiane. Nelle formule interpolatorie di Newton-Cotes (come ad esempio la regola del Trapezio o di Cavalieri-Simpson) i nodi x_{1}, …, x_{n} sono equispaziati e il grado di precisione δ è generalmente uguale almeno a n − 1 ma in alcuni casi, come per la regola di Cavalieri-Simpson, uguale al numero di nodi n. Vediamo ora formule che a parità di nodi hanno grado di precisione maggiore di n.

Sia w : (a, b) → ℝ (non necessariamente limitato) è una funzione peso, cioè tale che (cf. [1, p.206, p.270])

1. w è nonnegativa in (a, b);

2. w è integrabile in [a, b];

3. esista e sia finito

∫_{a}^{b} |x|^{n}w(x) dx

per ogni n ∈ ℕ;

4. se

∫_{a}^{b} g(x)w(x) dx

per una qualche funzione nonnegativa g allora g ≡ 0 in (a, b).

Tra gli esempi più noti ricordiamo

1. Legendre: w(x) ≡ 1 in [a, b] limitato;

2. Jacobi: w(x) = (1 − x)^{α} (1 + x)^{β} in (−1, 1) per α, β ≥ −1;

3. Chebyshev: w(x) = \frac{1}{\sqrt{1−x^{2}}} in (−1, 1);

4. Laguerre: w(x) = exp (−x) in [0, ∞);

5. Hermite: w(x) = exp (−x^{2}) in (−∞, ∞);

Nota 5.1. Diamo di seguito l’idea sul perchè si utilizzino funzioni peso w per il calcolo di integrali. Quanto detto non deve essere preso letteralmente ma con un po’ di elasticitá. Se dobbiamo approssimare ∫_{a}^{b} g(x) dx con g regolare allora \underline{di solito} non sarà difficile calcolare numericamente il valore dell’integrale, mentre se g è poco regolare sarà molto più oneroso ottenere buoni risultati. A tal proposito introduciamo delle stime d’errore che verranno discusse in una sezione successiva. Se (a, b) è limitato allora

|I(g) − I_{n}(g)| ≤ 2(b − a) min_{q_{n}∈P_{n}} ∥g − q_{n}∥_{∞}.

Essendo min_{q_{n}∈P_{n}} ∥g −q_{n}∥_{∞} l’elemento di miglior approssimazione, ricordiamo il teorema di Jackson

Teorema 5.2. Se f ∈ C^{k}([a, b]), ed f^{(k)} è α hölderiana, cioè

sup_{x,y∈[a,b]} |f^{(k)}(x) − f^{(k)}(y)| ≤ M|x − y|^{α}

per qualche M > 0, 0 < α ≤ 1. Allora esiste una costante M_{k} indipendente da f e n per cui

inf_{p∈P_{n}} ∥f − p∥_{∞} ≤ \frac{M_{k}}{n^{k+α}} , n ≥ 1.

La disuguaglianza del Teorema 5.2 dice che al crescere della regolarità k asintoticamente al crescere di n si ha un errore di migliore approssimazione inferiore.

Introducendo una funzione peso w, e supposto che f = g/w sia regolare, si considera il calcolo di

∫_{a}^{b} g(x) dx = ∫_{a}^{b} \frac{g(x)}{w(x)}w(x) dx = ∫_{a}^{b} f(x)w(x) dx.

La regolarità di f e la particolarità del calcolo dipendentemente dalla funzione peso permette di approssimare meglio ∫_{a}^{b} g(x) dx = ∫_{a}^{b} f(x)w(x) dx. Infatti per I(f) = ∫_{a}^{b} f(x)w(x) dx e I_{n}(f) un’appropriata formula di quadratura (che verrà detta gaussiana rispetto w) abbiamo

|I(f) − I_{n}(f)| ≤ 2∥w∥_{1} min_{q_{n}∈P_{n}} ∥f − q_{n}∥_{∞}.

e quindi, nuovamente dalle stime di Jackson, la facilità di calcolo dell’integrale dipenderà esclusivamente dalla regolarità di f (e non da quella di g = fw che potrebbe essere notevolmente inferiore visto che a priori la funzione g potrebbe essere perfino discontinua).

Si supponga ora di dover calcolare per qualche funzione f : (a, b) → ℝ

I_{w}(f) := ∫_{a}^{b} f(x)w(x) dx.

Il problemaè evidentemente più generale di quello di calcolare un integrale del tipo ∫_{a}^{b} f(x)dx con f ∈ C([a, b]), [a, b] limitato, visto che l’integranda fw noné necessariamente continua in [a, b] (si consideri ad esempio il peso di Chebyshev che ha una singolarità in a = −1, b = 1) oppure può succedere che l’intervallo sia illimitato come nel caso del peso di Laguerre o Hermite.

Esistono nuovamente x_{1}, …, x_{n} e pesi w_{1}, …, w_{n} (detti di Gauss-nome funzione peso) per cui le relative formule di quadratura di tipo interpolatorio abbiano grado di precisione almeno δ > n, cioè calcolino esattamente

∫_{a}^{b} p_{m}(x)w(x) dx

Figura 5.1. Adrien-Marie Legendre (1752-1833), Charles Hermite (1822-1901), Carl Gustav Jacob Jacobi (1804-1851) e Edmond Nicolas Laguerre (1834-1886).

per m > n?

La risposta è affermativa, come si può vedere in [1, p.272].

Teorema 5.3. Per ogni n ≥ 1 esistono e sono unici dei nodi x_{1}, …, x_{n} e pesi w_{1}, …, w_{n} per cui il grado di precisione sia almeno 2n−1. I nodi sono gli zeri del polinomio ortogonale di grado n,

ϕ_{n}(x) = A_{n} · (x − x_{1}) · . . . · (x − x_{n})

e i corrispettivi pesi sono

w_{i} = ∫_{a}^{b} L_{i}(x)w(x)dx = ∫_{a}^{b} L_{i}(x)w(x)dx, i = 1, …, n.

Dimostrazione. Tratta da [3, p.209]. Per prima cosa mostriamo che in effetti con tale scelta dei nodi la formula interpolatoria ha grado di precisione almeno 2n−1, che i pesi sono univocamente determinati e positivi.

Sia p_{2n−1} ∈ P_{2n−1} e q_{n−1}, r_{n−1} ∈ P_{n−1} tali che

p_{2n−1} = q_{n−1}ϕ_{n} + r_{n−1}.

Allora poichè

• ∫_{a}^{b} q_{n−1}(x)ϕ_{n}(x)w(x)dx = (q_{n−1}, ϕ_{n})_{w} = 0, poichè ϕ_{n} è il polinomio ortogonale rispetto w di grado n; infatti essendo

(ϕ_{k}, ϕ_{n})_{w} = 0, k = 0 < n

necessariamente da q_{n−1} = ∑_{k=0}^{n−1} γ_{k}ϕ_{k} abbiamo

(q_{n−1}, ϕ_{n})_{w} = (∑_{k=0}^{n−1} γ_{k}ϕ_{k}, ϕ_{n})_{w}

= ∑_{k=0}^{n−1} γ_{k}(ϕ_{k}, ϕ_{n})_{w} = 0  (5.1)

• la formula è interpolatoria per costruzione (vedere la definizione dei pesi!), per cui esatta per ogni polinomio di grado n − 1 in quanto basata su n punti a due a due distinti;

• se x_{k} è uno zero di ϕ_{n} allora
p_{2n−1}(x_{k}) = q_{n−1}(x_{k})ϕ_{n}(x_{k}) + r_{n−1}(x_{k}) = r_{n−1}(x_{k}).

abbiamo

∫_{a}^{b} p_{2n−1}(x)w(x)dx = ∫_{a}^{b} q_{n−1}(x)ϕ_{n}(x)w(x)dx + ∫_{a}^{b} r_{n−1}(x)w(x)dx
= 0 + ∫_{a}^{b} r_{n−1}(x)w(x)dx
= ∑_{k=1}^{n} w_{k}r_{n−1}(x_{k})

= ∑_{k=1}^{n} w_{k}p_{2n−1}(x_{k})  (5.2)

per cui la formula di Gauss ha grado di precisione almeno 2n − 1. Inoltre, come dimostrato da Stieltjes nel 1884, i pesi sono positivi, poichè in particolare la formula è esatta per ognuno dei quadrati dei polinomi di Lagrange relativo ai punti x_{1}, …, x_{n} per cui

0 < ∫_{a}^{b} L_{j}^{2}(x)w(x)dx = ∑_{k=1}^{n} w_{k}L_{j}^{2}(x_{k}) = w_{j}.

in quanto deg(L_{i}^{2})=2(n − 1) e la formula ha grado di precisione almeno 2n − 1. Se esistesse un’altra formula interpolatoria con grado di precisione almeno 2n − 1 e avesse nodi {x̃_{j}}_{j=1,...,n} e pesi {w̃_{j}}_{j=1,...,n} per prima cosa i pesi sarebbero positivi poichè il grado di precisione è almeno 2n − 1 e quindi sarebbe esatta per il j-simo polinomio di Lagrange L̃_{j} da cui

0 < ∫_{a}^{b} L̃_{j}^{2}(x)w(x)dx = ∑_{k=1}^{n} w̃_{k}L̃_{j}^{2}(x̃_{k}) = w̃_{j}.

D’altra parte se L̃_{j} è il j-simo polinomio di Lagrange (avente grado n − 1), poichè ϕ_{n} è il polinomio ortogonale di grado n rispetto al peso w, e w̃_{j} > 0 abbiamo che da

0 = (ϕ_{n}, L̃_{j})_{w} = ∫_{a}^{b} ϕ_{n}(x)L̃_{j}(x)w(x)dx = ∑_{k=1}^{n} w̃_{k}L̃_{j}(x̃_{k})ϕ_{n}(x̃_{k}) = w̃_{j} · ϕ_{n}(x̃_{j})

necessariamente x_{j} = x̃_{j} e visto che questo implica L_{j} = L̃_{j} ricaviamo anche

w_{j} = ∫_{a}^{b} L_{j}^{2}(x)w(x)dx = ∫_{a}^{b} L̃_{j}^{2}(x)w(x)dx = w̃_{j}

per cui la formula gaussiana cercata è unica.

8

6. Formule dell’errore. Per quanto riguarda gli errori compiuti da alcune delle formula di quadratura discusse si ha (cf.[1, p.264])

Teorema 6.1. Sia

I(f) ≈ I_{n}(f) = ∑_{i=0}^{n} w_{i,n}f(x_{i,n})

una regola di Newton-Cotes.

1. se n è pari e f ∈ C^{(n+2)}([a, b]) allora

I(f) − I_{n}(f) = C_{n} h^{n+3}f^{(n+2)}(η), η ∈ (a, b)

con

C_{n} = \frac{1}{(n + 2)!} ∫_{0}^{n} μ^{2}(μ − 1) … (μ − n)dμ;

2. se n è dispari e f ∈ C^{(n+1)}([a, b]) allora

I(f) − I_{n}(f) = C_{n} h^{n+2}f^{(n+1)}(η), η ∈ (a, b)

con

C_{n} = \frac{1}{(n + 1)!} ∫_{0}^{n} μ(μ − 1) … (μ − n)dμ;

Si osserva facilmente che quanto visto in precedenza per la regola del trapezio e la regola di Cavalieri-Simpson, è consistente con questi due teoremi. Inoltre, da questi ultimi, si ottengono gli errori di regole composte da un numero di punti minore o uguale a 7 (per motivi di stabilità non si suggeriscono regole con più punti).

Per quanto concerne le formule gaussiane (cf.[1, p.272], cf.[3, p.264], cf.[5, p.344]) ricordiamo il seguente teorema di Markov

Teorema 6.2. Sia f ∈ C^{(2n)}(a, b) con (a, b) compatto e supponiamo

I_{w}(f) = ∫_{a}^{b} f(x)w(x)dx ≈ I_{n}(f) = ∑_{i=1}^{n} w_{i,n}f(x_{i,n})

sia una formula gaussiana rispetto alla funzione peso w. Allora

E_{n}(f) := I_{w}(f) − L_{n}(f) = \frac{γ_{n}}{A_{n}^{2}(2n)!}f^{(2n)}(η), η ∈ (a, b)

dove A_{n} è il coefficiente di grado massimo del polinomio ortogonale ϕ_{n} di grado n, γ_{n} = ∫_{a}^{b} ϕ_{n}^{2}(x)w(x)dx.

In particolare, se w ≡ 1, [a, b] ≡ [−1, 1] allora

E_{n}(f) = \frac{2^{2n+1}(n!)^{4}}{(2n + 1)[(2n)!]^{3}} f^{(2n)}(η), η ∈ (−1, 1).

Figura 6.1. Grafico in scala semilogaritmica della funzione \frac{2^{2n+1}(n!)^{4}}{(2n+1)[(2n)!]^{3}} .

7. Stabilità di una formula di quadratura. Sia (a, b) un intervallo non necessariamente compatto e w una funzione peso in (a, b). Inoltre supponiamo

I_{w}(f) := ∫_{a}^{b} f(x)w(x) dx ≈ I_{n}(f) := ∑_{j=1}^{n} w_{j}f_{j}, con f_{j} = f(x_{j})  (7.1)

e che invece di {f_{j}}_{j} si disponga di una loro approssimazione {f̃_{j}}_{j}. Di conseguenza al posto di I_{n}(f) si calcola

Ĩ_{n}(f) = ∑_{j=1}^{n} w_{j}f̃_{j},

ed è

|I_{n}(f) − Ĩ_{n}(f)| = | ∑_{j=1}^{n} w_{j}(f_{j} − f̃_{j})| ≤ ∑_{j=1}^{n} |w_{j}||f_{j} − f̃_{j}|

≤ (∑_{j=1}^{n} |w_{j}|) · max_{j} |f_{j} − f̃_{j}|  (7.2)

Quindi la quantità

∑_{j=1}^{n} |w_{j}|

è un indice di stabilità della formula di quadratura.

Si dimostra che se (a, b) è limitato allora in effetti l’operatore lineare I_{n} : C([a, b]) → ℝ è continuo in quanto per il teorema di Weierstrass esiste ∥f∥_{∞} ed è

|I_{n}(f)| = |∑_{j=1}^{n} w_{j}f_{j}| ≤ ∑_{j=1}^{n} |w_{j}||f_{j}|

≤ (∑_{j=1}^{n} |w_{j}|) · max_{j} |f_{j}| ≤ (∑_{j=1}^{n} |w_{j}|) · ∥f∥_{∞}  (7.3)

per cui la norma dell’operatore di quadratura

∥I_{n}∥_{∞} = max_{f∈C([a,b]),f≠0} \frac{|I_{n}(f)|}{∥f∥_{∞}}

coincide con ∑_{j=1}^{n} |w_{j}|.

Mostriamo ora il seguente teorema attribuito a Stieltjes, che lega l’errore delle formule di quadratura a quello di miglior approssimazione.

Teorema 7.1. Sia (a, b) un intervallo limitato e si supponga w : (a, b) → ℝ sia una funzione peso e

I_{n}(f) = ∑_{j=1}^{n} w_{j}f_{j}, con f_{j} = f(x_{j})

una formula di quadratura avente grado di precisione almeno n. Allora

|I(f) − I_{n}(f)| ≤ (∥w∥_{1} + ∥I_{n}∥_{∞}) · min_{q_{n}∈P_{n}} ∥f − q_{n}∥_{∞}.

Dimostrazione. Se q_{n} ∈ P_{n} è un polinomio arbitrario di grado n ed I(f) = I_{n}(f) avendo la formula di quadratura grado di precisione almeno n, essendo

|I(f)| = |∫_{a}^{b} f(x)w(x)dx| ≤ max_{x∈[a,b]} |f(x)| · ∫_{a}^{b} w(x)dx = ∥f∥_{∞}∥w∥_{1}

poichè per definizione

|I_{n}(f)| ≤ ∥I_{n}∥_{∞}∥f∥_{∞}

abbiamo

|I(f) − I_{n}(f)| = |I(f) − I_{n}(q) + I_{n}(q) − I_{n}(f)|
≤ |I(f) − I_{n}(q)| + |I_{n}(q) − I_{n}(f)|
≤ |I(f) − I(q)| + |I_{n}(q − f)|
≤ |I(f − q)| + |I_{n}(f − q)|
≤ ∥w∥_{1}∥f − q∥_{∞} + ∥I_{n}∥_{∞}∥f − q∥_{∞}
= (∥w∥_{1} + ∥I_{n}∥_{∞}) · ∥f − q∥_{∞}  (7.4)

L’interesse di questo teorema è il legame col polinomio di miglior approssimazione. Se i pesi sono positivi, allora ∥I_{n}∥_{∞} = ∥w∥_{1} e si ha

|I(f) − I_{n}(f)| ≤ 2∥w∥_{1} min_{q_{n}∈P_{n}} ∥f − q_{n}∥_{∞}.

Per stime delle quantità

min_{q_{n}∈P_{n}} ∥f − q_{n}∥_{∞}

si devono usare i teoremi di Jackson per l’errore del polinomio di miglior approssimazione di una funzione f ∈ C([a, b]) (dotando C([a, b])) della norma infinito).

Così ad esempio, se usiamo la funzione peso di Legendre w ≡ 1 nell’intervallo (−1, 1) si ha che

|I(f) − I_{n}(f)| ≤ 4 · min_{q_{n}∈P_{n}} ∥f − q_{n}∥_{∞}

Vista l’arbitrarietà di q_{n}, scelto quale q_{n} il polinomio di miglior approssimazione in norma infinito, si ottiene l’asserto.

Risulta importante osservare che il teorema precedente che ha quali contributi il prodotto di due termini. Il primo è dovuto alla funzione peso e alla stabilità della formula di quadratura, mentre il secondo è dato esclusivamente dalla miglior approssimazione di f (e non fw)! Quindi se w è una funzione peso con fw non regolare ma f regolare allora l’utilizzo di formule gaussiane rispetto alla funzione peso w, come anticipato prima, offre risultati potenzialmente migliori.

Quale esempio si pensi di dover calcolare l’integrale

∫_{−1}^{1} exp (x)\sqrt{1 − x^{2}}dx

con la formula di Gauss-Legendre e una formula di Gauss-Jacobi con esponenti α = 1/2 e β = 0. Quale delle due sarà da usare e perchè ?

Nota 7.2. Supponiamo che alcuni pesi siano negativi. Siano {w_{j}^{+}}_{j=1,...,n_{+}} e {w_{l}^{−}}_{l=1,...,n_{−}} rispettivamente i pesi positivi e negativi. Di conseguenza

∫_{a}^{b} w(x)dx = ∑_{k=1}^{n} w_{k} = ∑_{j=1}^{n_{+}} w_{j}^{+} + ∑_{l=1}^{n_{−}} w_{l}^{−}  (7.5)

∥I_{n}∥_{∞} = ∑_{k=1}^{n} |w_{k}| = ∑_{j=1}^{n_{+}} |w_{j}^{+}| + ∑_{l=1}^{n_{−}} |w_{l}^{−}| = ∑_{j=1}^{n_{+}} w_{j}^{+} − ∑_{l=1}^{n_{−}} w_{l}^{−}  (7.6)

da cui

∥I_{n}∥_{∞} = ∑_{j=1}^{n_{+}} w_{j}^{+}−∑_{l=1}^{n_{−}} w_{l}^{−} = (∫_{a}^{b} w(x)dx − ∑_{j=1}^{n_{+}} w_{j}^{+})−∑_{l=1}^{n_{−}} w_{l}^{−} = ∫_{a}^{b} w(x)dx−2 ∑_{l=1}^{n_{−}} w_{l}^{−},

mentre se i pesi fossero tutti positivi avremmo

∥I_{n}∥_{∞} = ∑_{j=1}^{n_{+}} w_{j}^{+} − ∑_{l=1}^{n_{−}} w_{l}^{−} = ∑_{j=1}^{n_{+}} w_{j}^{+} = ∫_{a}^{b} w(x)dx

per cui maggiore è −2 ∑_{l=1}^{n_{−}} w_{l}^{−} > 0 e peggiore è la stabilità della formula di quadratura.

8. Teoremi di convergenza. Consideriamo una formula di quadratura del tipo

I_{w}(f) := ∫_{a}^{b} f(x)w(x)dx ≈ I_{n}(f) := ∑_{i=0}^{n} w_{i,n}f(x_{i,n})  (8.1)

Figura 8.1. George Polya (1887-1985) e Vladimir Andreevich Steklov (1864-1926).

con al solito w una funzione peso definita nell’intervallo (a, b). Se tale intervallo è limitato ed f continua in [a, b] allora per quanto visto in precedenza si ha fw ∈ L^{1}(a, b). Sia

E_{n}(f) := ∫_{a}^{b} f(x)w(x)dx − ∑_{i=0}^{n} w_{i,n}f(x_{i,n}).

Dimostriamo ora il teorema di Polya-Steklov [3, p.202], lasciando al lettore che ha conoscenze di analisi funzionale la parte di dimostrazione che utilizza il celebrato teorema di Banach-Steinhaus.

Teorema 8.1. Siano x_{i,k}_{i,k} dei punti di un intervallo compatto [a, b]. Condizione necessaria e sufficiente affinchè per ogni funzione continua f ∈ C([a, b]) si abbia

lim_{n→+∞} E_{n}(f) = 0

è che sia 1. esiste M ∈ ℝ tale che per ogni n si abbia

∑_{i=0}^{n} |w_{i,n}| ≤ M;

2. per ogni k ∈ ℕ si abbia

lim_{k→+∞} E_{n}(x^{k}) = 0.

Dimostrazione. Supponiamo che

1. esiste M ∈ ℝ tale che per ogni n si abbia

∑_{i=0}^{n} |w_{i,n}| ≤ M;

2. per ogni k ∈ ℕ si abbia

lim_{k→+∞} E_{n}(x^{k}) = 0.

Per un teorema di densità dovuto a Weierstrass, per ogni ϵ > 0 esiste un polinomio p tale che ∥f − p_{n}∥_{∞} ≤ ϵ e quindi

|E_{n}(f − p_{n})| ≤ |∫_{a}^{b}(f(x) − p_{n}(x))w(x)dx| + |∑_{i=0}^{n} w_{i,n}(f(x_{i,n}) − p_{n}(x_{i,n}))|

≤ ∥f − p_{n}∥_{∞}∥w∥_{1} + ∑_{i=0}^{n} |w_{i,n}||f(x_{i,n}) − p_{n}(x_{i,n})|

≤ ∥f − p_{n}∥_{∞} · ∥w∥_{1} + (∑_{i=0}^{n} |w_{i,n}|) · ∥f − p_{n}∥_{∞}

= (∥w∥_{1} + ∑_{i=0}^{n} |w_{i,n}|) · ϵ.  (8.2)

Di conseguenza

|E_{n}(f)| ≤ |E_{n}(f) − E_{n}(p)| + |E_{n}(p)|
≤ (∥w∥_{1} + ∑_{i=0}^{n} |w_{i,n}|) · ϵ + |E_{n}(p)|
≤ (∥w∥_{1} + M) · ϵ + |E_{n}(p)|

(8.3)

e siccome |E_{n}(p)| → 0 al crescere di n poichè se p(x) = ∑_{k=0}^{m} a_{k}x^{k} allora

E_{n}(p) = ∑_{k=0}^{m} a_{k}E_{n}(x^{k}) → 0

abbiamo

lim_{n} |E_{n}(f)| ≤ (∥w∥_{1} + M) · ϵ + 0 = (∥w∥_{1} + M) · ϵ

da cui per l’arbitrarietà di ϵ abbiamo lim_{n} |E_{n}(f)| = 0 e quindi per la continuità della funzione | · | abbiamo

lim_{n} |E_{n}(f)| = | lim_{n} E_{n}(f)| = 0

che permette di concludere lim_{n} E_{n}(f) = 0.

Viceversa supponiamo lim_{n} E_{n}(f) = 0 per ogni f ∈ C([a, b]). Si vede che

|E_{n}(f)| = | ∫_{a}^{b} f(x)w(x)dx − ∑_{i=0}^{n} w_{i,n}f(x_{i,n})|

≤ ∫_{a}^{b} |f(x)|w(x)dx + | ∑_{i=0}^{n} w_{i,n}f(x_{i,n})|

≤ ∥w∥_{1}∥f∥_{∞} + ∑_{i=0}^{n} |w_{i,n}||f(x_{i,n})|

≤ ∥w∥_{1}∥f∥_{∞} + (∑_{i=0}^{n} |w_{i,n}|)∥f∥_{∞}

= (∥w∥_{1} + ∑_{i=0}^{n} |w_{i,n}|) · ∥f∥_{∞}  (8.4)

In realtà si prova che

∥E_{n}∥_{∞} = max_{f∈C([a,b]),f≠0} \frac{|E_{n}(f)|}{∥f∥} = ∥w∥_{1} + ∑_{i=0}^{n} |w_{i,n}|.

Il teorema di uniforme limitatezza (talvolta citato come di Banach-Steinhaus) [2, p.58] stabilisce che se L_{n} è una sequenza di operatori lineari limitati da uno spazio di Banach V a uno spazio di Banach W e per ogni v ∈ V la sequenza {L_{n}(v)}_{n} è limitata allora

sup_{n} ∥L_{n}∥ < +∞.

Nel nostro caso • V ≡ (C([a, b]), ∥ · ∥_{∞}), W ≡ ℝ sono spazi di Banach,

• posto L_{n} ≡ E_{n}, se f ∈ C([a, b]) abbiamo che E_{n}(f) è limitata in quanto per ipotesi

è infinitesima, quindi per il teorema di Banach-Steinhaus si ha

sup_{n} ∥E_{n}∥_{∞} < +∞.

Essendo

sup_{n} ∑_{i=0}^{n} |w_{i,n}| ≤ sup_{n} (∥w∥_{1} + ∑_{i=0}^{n} |w_{i,n}|) = sup_{n} ∥E_{n}∥_{∞} < +∞,

esiste M finito ed indipendente da n tale che ∑_{i=0}^{n} |w_{i,n}| ≤ M.

Il secondo punto da dimostrare è ovvio in quanto per ogni k, si ha x^{k} ∈ C([a, b])

Notiamo che

1. L’intervallo [a, b] è limitato per cui il teorema di Polya non è applicabile per funzioni peso quali Gauss-Laguerre e Gauss-Hermite.

2. Si consideri una suddivisione Δ_{m} = {τ_{i}}_{i=0,...,m} dell’intervallo (a, b) con τ_{i} < τ_{i+1}, τ_{0} = a, τ_{m} = b; si supponga che h_{i} = τ_{i+1} − τ_{i} e che sia

h_{Δ_{m}} = max_{k=0,...,m−1} |h_{k}|.

Se l’l-simo intervallo (con l = 1, …, m) della suddivisione contiene i punti

τ_{l−1} = x_{(l−1)·n} < x_{(l−1)·n+1} < … < x_{l·n} = τ_{l}

sia s_{Δ_{m},n}(f) la funzione polinomiale a tratti di grado n ottenuta interpolando nell’isimo subintervallo con un polinomio di grado n le coppie

(x_{(l−1)·n}, f(x_{(l−1)·n})), …, (x_{l·n}, f(x_{l·n})), i = 1, …, m.

Come è noto, se f ∈ C^{(n+1)}(a, b), per qualche ξ ∈ (a, b)

∥f − s_{Δ_{m},n}∥_{∞} ≤ \frac{h^{(n+1)}f^{(n+1)}(ξ)}{(n + 1)!} ≤ \frac{h^{(n+1)}∥f^{(n+1)}∥_{∞}}{(n + 1)!} .

Se consideriamo una formula composta basata sull’applicare una formula di Newton-Cotes su n + 1 pesi positivi in ognuno degli m subintervalli, si vede facilmente che tutti i suoi pesi sono positivi (perchè ?) e siccome in ogni sottointervallo la formula interpolatoria ha grado di precisione almeno 0, abbiamo

∑_{i=0}^{n} |w_{i,n}| = ∑_{i=0}^{n} w_{i,n} = ∑_{i=0}^{n} w_{i,n} · 1 = ∫_{a}^{b} 1 dx = b − a.

Quindi, relativamente alla prima ipotesi del teorema di Polya-Steklov, essa è verificata per M = b − a. Per quanto riguarda la verifica della seconda ipotesi, tale formula composta integra esattamente ogni funzione polinomiale a tratti di grado n su Δ e quindi dalla formula dell’errore dell’interpolante polinomiale a tratti si ha

∫_{a}^{b} f(x)dx = ∫_{a}^{b}(f(x) − s_{Δ_{m},n})dx + ∫_{a}^{b} s_{Δ_{m},n}dx
= ∫_{a}^{b}(f(x) − s_{Δ_{m},n})dx + ∑_{i=0}^{n} w_{i,n}f(x_{i,n})dx  (8.5)

da cui

|E_{n}(f)| = |∫_{a}^{b} f(x)dx − ∑_{i=0}^{n} w_{i,n}f(x_{i,n})|

= |∫_{a}^{b}(f(x) − s_{Δ_{m},n})dx|

≤ |∫_{a}^{b} 1 · dx| · ∥f − s_{Δ_{m},n}∥_{∞}
≤ (b − a) · ∥f − s_{Δ_{m},n}∥_{∞}
≤ (b − a) · \frac{h^{(n+1)}∥f^{(n+1)}∥_{∞}}{(n + 1)!}

(8.6)

Siccome un polinomio f(x) = x^{k} è infinitamente derivabile, f^{(n+1)} è continua in [a, b] e quindi per il teorema di Weierstrass ∥f^{(n+1)}∥_{∞} è finito. Di conseguenza, se la successione di formule composte è tale che la massima suddivisione h tende a 0 allora E_{n}(x^{k}) → 0. Per il teorema di Polya-Steklov si ha così che la formula composta

∫_{a}^{b} f(x)w(x)dx ≈ ∑_{i} w_{i,n}f(x_{i,n})  (8.7)

è tale che qualsiasi sia la funzione continua f, E_{m}(f) → 0 quando la massima ampiezza dei subintervalli tende a 0.

3. Se consideriamo una formula di Gauss, su un dominio limitato, i nodi w_{i,n} sono positivi ed avendo grado di precisione almeno 1

∑_{i=0}^{n} |w_{i,n}| = ∑_{i=0}^{n} w_{i,n} = ∑_{i=0}^{n} w_{i,n} · 1 = ∫_{a}^{b} 1 · w(x)dx = ∥w∥_{1}.

Quindi posto M = ∥w∥_{1}, la prima ipotesi del teorema di Polya è verificata. Per quanto riguarda la seconda ipotesi, essendo il grado di precisione almeno 2n + 1, fissato k, per n ≥ ceil((k − 1)/2) si ha E_{n}(x^{k}) = 0. Di conseguenza, essendo tutti gli zeri contenuti in (a, b), possiamo applicare il teorema di Polya-Steklov e dedurre che al crescere del numero di punti n + 1 della formula gaussiana si ha che

lim_{n→+∞} E_{n}(f) = 0

qualsiasi sia la funzione continua f ∈ C(a, b).

8.1. Implementazione Matlab di alcune formule composte. Mostreremo di seguito un’implementazione in Matlab/Octave della formula composta dei trapezi e di Cavalieri-Simpson.

function [x,w]=trapezi_composta(N,a,b)
% FORMULA DEI TRAPEZI COMPOSTA.
% INPUT:
% N: NUMERO SUBINTERVALLI.
% a, b: ESTREMI DI INTEGRAZIONE.
% OUTPUT:
% x: NODI INTEGRAZIONE.
% w: PESI INTEGRAZIONE (INCLUDE IL PASSO!).
h=(b-a)/N;               % PASSO INTEGRAZIONE.
x=a:h:b; x=x';          % NODI INTEGRAZIONE.
w=ones(N+1,1);          % PESI INTEGRAZIONE.
w(1)=0.5; w(N+1)=0.5;
w=w*h;

La funzione trapezi composta appena esposta calcola i nodi e i pesi della omonima formula composta. L’unica difficoltà del codice consiste nel calcolo dei pesi w. Essendo per loro definizione

I(f) ≈ S_{1}^{c}(f) := ∑_{i=0}^{N} w_{i}f(x_{i})  (8.8)

Figura 8.2. Formula dei trapezi composta per il calcolo di ∫_{0.5}^{2} sin (x) dx (area in magenta).

come pure per (4.2)

S_{1}^{(c)} := h [\frac{f(x_{0})}{2} + f(x_{1}) + … + f(x_{N−1}) + \frac{f(x_{N})}{2}]  (8.9)

deduciamo che w_{0} = w_{N} = h/2 mentre w_{1} = … = w_{N−1} = h, cosa che giustifica le ultime linee della function trapezi composta.

Si potrebbe usare il comando Matlab trapz nella sua implementazione

>> help trapz
 TRAPZ  Trapezoidal numerical integration.
     Z = TRAPZ(Y) computes an approximation of the integral of Y via
     the trapezoidal method (with unit spacing).  To compute the integral
     for spacing different from one, multiply Z by the spacing increment.
     For vectors, TRAPZ(Y) is the integral of Y.
     ... ... ...

e sostituire la parte relativa al calcolo dei pesi con

I=h*trapz(fx);

Vediamone i dettagli in Matlab (versione 6.1) per il calcolo di

∫_{0}^{1} sin(x)dx = − cos(1) − (− cos(0)) = − cos(1) + 1 ≈ 0.45969769413186.

sia utilizzando la funzione trapz che trapezi composta

>> format long;
>> [x,w]=trapezi_composta(10,0,1);
>> fx=sin(x);
>> I_trapezi_composta=w'*fx
I_trapezi_composta =
   0.45931454885798
>> h=(1-0)/10;
>> I_trapz=h*trapz(fx)

I_trapz =
   0.45931454885798
>>

Di conseguenza per implementare la regolaè del tutto equivalente usare la function trapezi composta o trapz. Si osserva che è sbagliato chiamare trapz senza il passo h (nell’esempio non si dividerebbe per 10, e invece di 0.45931454885798 avremmo 4.5931454885798). Evitiamo il diretto utilizzo di trapz perchè non presente in alcune vecchie versioni di Octave (ma non nella più recente 2.1.73).

Per quanto riguarda la formula di Cavalieri-Simpson composta

function [x,w]=simpson_composta(N,a,b)
% FORMULA DI SIMPSON COMPOSTA.
% INPUT:
% N: NUMERO SUBINTERVALLI.
% a, b: ESTREMI DI INTEGRAZIONE.
% OUTPUT:
% x: INTEGRAZIONE.
% w: PESI INTEGRAZIONE (INCLUDE IL PASSO!).
h=(b-a)/N;               % AMPIEZZA INTERVALLO.
x=a:(h/2):b; x=x';      % NODI INTEGRAZIONE.
w=ones(2*N+1,1);             % PESI INTEGRAZIONE.
w(3:2:2*N-1,1)=2*ones(length(3:2:2*N-1),1);
w(2:2:2*N,1)=4*ones(length(2:2:2*N),1);
w=w*h/6;

Similmente alla routine per il calcolo dei nodi e i pesi di trapezi composta, le ultime righe sono le più difficili da capire, ma un confronto con (8.8) e (4.3) ne spiega il significato.

Una volta noti il vettore (colonna) x dei nodi e w dei pesi di integrazione, se la funzione f è richiamata da un m-file f.m, basta

fx=f(x);                   % VALUT. FUNZIONE.
I=w'*fx;                   % VALORE INTEGRALE.

per calcolare il risultato fornito dalla formula di quadratura composta. Ricordiamo che se w = (w_{k})_{k=0,...,N} ∈ ℝ^{N+1}, fx = (f(x_{k}))_{k=0,...,N} ∈ ℝ^{N+1} sono due vettori colonna allora il prodotto scalare

w ∗ fx := ∑_{k=0}^{N} w_{k} · fx_{k} := ∑_{k=0}^{N} w_{k} · f(x_{k})

si scrive in Matlab/Octave come w'*fx. Osserviamo che dimensionalmente il prodotto di un vettore 1 × (N + 1) con un vettore (N + 1) × 1 dà uno scalare (cioè un vettore 1 × 1).

Applichiamo ora la formula composta di Cavalieri-Simpson all’esempio precedente:

>> format long;

>> [x,w]=simpson_composta(5,0,1);
>> fx=sin(x);
>> I_simpson=w'*fx
I_simpson =
   0.45969794982382
>> length(x)
ans =
     11
>>

Facciamo ora un altro esempio. Calcoliamo numericamente utilizzando le formule composte sopra citate

∫_{−1}^{1} x^{20}dx = (1^{21}/21) − (−1)^{21}/21 = 2/21 ≈ 0.09523809523810.

A tal proposito scriviamo il seguente codice demo composte.

N=11;                                   %SCEGLIERE DISPARI.
a=-1; b=1;
N_trap=N-1;
[x_trap,w_trap]=trapezi_composta(N_trap,a,b);
fx_trap=x_trap.^20;                   % VALUT. FUNZIONE.
I_trap=w_trap'*fx_trap;               % TRAPEZI COMPOSTA.
N_simpson=(N-1)/2;
[x_simp,w_simp]=simpson_composta(N_simpson,a,b)
fx_simp=x_simp.^20;                   % VALUT. FUNZIONE.
I_simp=w_simp'*fx_simp;               % SIMPSON COMPOSTA.
fprintf('\n \t [TRAPEZI COMPOSTA] [PTS]: %4.0f', length(x_trap));
fprintf('\n \t [TRAPEZI COMPOSTA] [RIS]: %14.14f', I_trap);
fprintf('\n \t [SIMPSON COMPOSTA] [PTS]: %4.0f', length(x_simp));
fprintf('\n \t [SIMPSON COMPOSTA] [RIS]: %14.14f', I_simp);

ottenendo

           [TRAPEZI COMPOSTA] [PTS]:   11
           [TRAPEZI COMPOSTA] [RIS]: 0.20462631505024
           [SIMPSON COMPOSTA] [PTS]:   11
           [SIMPSON COMPOSTA] [RIS]: 0.13949200364447

Si può vedere che usando formule di tipo gaussiano (cf. [28], [29]) o di tipo Clenshaw-Curtis (cf. [26], [11], [13]) a parità di valutazioni della funzione f avremmo ottenuto

           [GAUSS-LEGENDRE  ]: 0.095238095238095649
           [CLENSHAW-CURTIS ]: 0.094905176204004307

col costo aggiuntivo di dover calcolare tramite complicati algoritmi i pesi e i nodi di Gauss o i nodi di Clenshaw-Curtis via un IFFT [27].

8.2. Formule gaussiane in Matlab. Eccetto che in pochi casi (come ad esempio per la funzione peso di Chebyshev), non esistono formule esplicite per l’individuazione di tali nodi e pesi. Una volta venivano tabulati, oggi si consiglia di applicare del software che si può trovare ad esempio nella pagina di Walter Gautschi

http://www.cs.purdue.edu/people/wxg

Fissato un peso, ad esempio quello di Jacobi, si cercano per prima cosa i coefficienti di ricorrenza, che possono essere calcolati con l’ m-file r jacobi.m nel sito di W. Gautschi. La sintassi è la seguente

ab=r_jacobi(N,a,b)

Come si vede dai commenti al file

% R_JACOBI Recurrence coefficients for monic Jacobi polynomials.
%
%     ab=R_JACOBI(n,a,b) generates the first n recurrence
%     coefficients for monic Jacobi polynomials with parameters
%     a and b. These are orthogonal on [-1,1] relative to the
%     weight function w(t)=(1-t)^a(1+t)^b. The n alpha-coefficients
%     are stored in the first column, the n beta-coefficients in
%     the second column, of the nx2 array ab. The call ab=
%     R_JACOBI(n,a) is the same as ab=R_JACOBI(n,a,a) and
%     ab=R_JACOBI(n) the same as ab=R_JACOBI(n,0,0).
%
%     Supplied by Dirk Laurie, 6-22-1998; edited by Walter
%     Gautschi, 4-4-2002.

a, b corrispondono rispettivamente all’α e β delle formule di Jacobi (e non agli estremi di integrazione!). I coefficienti di ricorrenza sono immagazzinati nella variabile ab. Ricordiamo che se {ϕ_{n}(x)} è una famiglia di polinomi ortogonali su [a, b], rispetto ad una funzione peso w(x), allora per n ≥ 1

ϕ_{n+1}(x) = α_{n}(x − β_{n})ϕ_{n}(x) − γ_{n}ϕ_{n−1}(x)

ove, detto A_{n} il coefficiente di grado n in ϕ_{n}, si ha α_{n} = A_{n+1}/A_{n},

β_{n} = \frac{(xϕ_{n}, ϕ_{n})}{∥ϕ_{n}∥^{2}}, γ_{n} = \frac{(ϕ_{n}, xϕ_{n−1})}{∥ϕ_{n}∥^{2}},

con (f, g) = ∫_{a}^{b} f(x)g(x) dx. Nel caso di r jacobi si usa la formula ricorsiva dei polinomi ortogonali (monici!) [7, p.216]

ϕ_{n+1}(x) = (x − \overline{α}_{n})ϕ_{n}(x) − \overline{β}_{n}ϕ_{n−1}(x), k = 1, 2, …
ϕ_{−1}(t) = 0, ϕ_{0}(t) = 1.  (8.10)

Osserviamo che, essendo i polinomi ortogonali definiti a meno di una costante moltiplicativa non nulla (se (ϕ_{k},ϕ_{j})= 0 per ogni j < k pure (τϕ_{k},ϕ_{j})= 0 per ogni j < k, per ogni τ ≠ 0), si può richiedere che la famiglia triangolare di polinomi ortogonali sia di polinomi monici (cioè con coefficienti di grado massimo uguali a 1). Il vettore ab ha quale prima colonna il vettore dei coefficienti \{\overline{α}_{k}\} e quale seconda colonna il vettore dei coefficienti \{\overline{β}_{k}\}.

Figura 8.3. Grafico che illustra la distribuzione dei 20 nodi e il valore dei 20 pesi di Gauss-Legendre nell’intervallo [1, 1]

.

Per ulteriori dettagli sui polinomi ortogonali si confronti [1, p.270-283], [3, p.978-985], [9, p.375-401]. A questo punto si chiama la funzione gauss.m (reperibile nuovamente presso lo stesso sito di W. Gautschi)

% GAUSS Gauss quadrature rule.
%
%     Given a weight function w encoded by the nx2 array ab of the
%     first n recurrence coefficients for the associated orthogonal
%     polynomials, the first column of ab containing the n alpha-
%     coefficients and the second column the n beta-coefficients,
%     the call xw=GAUSS(n,ab) generates the nodes and weights xw of
%     the n-point Gauss quadrature rule for the weight function w.
%     The nodes, in increasing order, are stored in the first
%     column, the n corresponding weights in the second column, of
%     the nx2 array xw.
%
function xw=gauss(N,ab)
N0=size(ab,1); if N0<N, error('input array ab too short'), end
J=zeros(N);
for n=1:N, J(n,n)=ab(n,1); end
for n=2:N
  J(n,n-1)=sqrt(ab(n,2));
  J(n-1,n)=J(n,n-1);
end
[V,D]=eig(J);
[D,I]=sort(diag(D));
V=V(:,I);
xw=[D ab(1,2)*V(1,:)'.^2];

che per un certo N, uguale al massimo al numero di righe della matrice di coefficienti di ricorrenza ab, fornisce nodi x e pesi w immagazzinati in una matrice xw che ha quale prima colonna x e quale seconda colonna w. Osserviamo che dall’help si evince che i polinomi ortogonali sono monici e quindi compatibili con quelli forniti da r jacobi.

La descrizione di perchè tale software fornisca il risultato desideratoè complicata ma può essere trovata nella monografia di W. Gautschi sui polinomi ortogonali, per Acta Numerica.

Conseguentemente per trovare i nodi e pesi relativi alla formula di Gauss-Legendre in (a, b), che è una formula di tipo Gauss-Jacobi per α = 0 e β = 0, nonchè calcolare con essi degli integrali di una funzione f, si procede come segue.

8.2.1. Una demo di Gauss-Jacobi. Salviamo nel file integrazione gauss jacobi.m

function [I_jac,x_jac,w_jac]=integrazione_gauss_jacobi(N,ajac,bjac,a,b,f)
% INPUT:
%
% N: GRADO DI GAUSS-JACOBI. SE SI VUOLE CALCOLARE I(f) CON
%                   I(f)=int_a^b f(x) dx
%     PORRE ajac=0, bjac=0 (FUNZIONE PESO DI LEGENDRE w(x)=1).
% ajac, bjac: PARAMETRI DI GAUSS-JACOBI, CIOE' SI INTEGRA:
%          I(fw)=int_a^b f(x) w(x) dx
%     CON w(x)=(1-x)^ajac (1+x)^bjac.
% a,b: ESTREMI INTERVALLO INTEGRAZIONE.
% f: FUNZIONE DA INTEGRARE.
%
% OUTPUT:
%
% I_jac: VALORE DI I(fw) APPROSSIMATO DA gauss_jacobi. PER ajac=0,
%         bjac=0 CORRISPONDE A I(f).
% x_jac: NODI GAUSS-JACOBI.
% w_jac: PESI GAUSS-JACOBI.
% ROUTINES ESTERNE: r_jacobi, gauss.

ab_jac=r_jacobi(N,ajac,bjac);         % TERM. RICORSIVI.
xw_jac=gauss(N,ab_jac);                % NODI E PESI IN MATRICE.
x_jac=xw_jac(:,1);                      % NODI GAUSS-LEGENDRE [-1,1].
x_jac_ab=((a+b)/2)+((b-a)/2)*x_jac;  % NODI GAUSS-LEGENDRE [a,b].
w_jac=xw_jac(:,2);                      % PESI GAUSS-LEGENDRE [-1,1].
w_jac_ab=((b-a)/2)*w_jac;              % PESI GAUSS-LEGENDRE [a,b].
fx_jac_ab=feval(f,x_jac_ab);          % VALUTAZIONE FUNZIONE.
I_jac=w_jac_ab'*fx_jac_ab;             % VALORE INTEGRALE.

e nel file f.m delle funzioni su cui effettueremo dei test. Un esempio è

function fx=f(x)
fx=x.^20;
% ALCUNE FUNZIONI CHE FANNO PARTE DEL SET STUDIATO NELL'ARTICOLO:
% "IS GAUSS QUADRATURE BETTER THAN CLENSHAW-CURTIS?"
% DI L.N. TREFETHEN.
% fx=exp(x);
% fx=exp(-x.^2);
% fx=1./(1+16*(x.^2));
% fx=exp(-x.^(-2));
% fx=abs(x); fx=fx.^3;

% fx=x.^(0.5);
% fx=exp(x).*(sqrt(1-x));

Altre funzioni test [11] sono

fx = x.^20
fx = exp(x)
fx = exp(-x.^2)
fx = 1./(1+16*(x.^2))
fx = exp(-x.^(-2))
fx = (abs(x)).^3;
fx=x.^(0.5);
fx=exp(x).*(sqrt(1-x));

con a = −1, b = 1 e per cambiare il tipo di funzione, basta modificare la posizione dei caratteri %. Scaricando dalla pagina web del corso i files integrazione gauss jacobi.m, f.m e dalla pagina di W. Gautschi r jacobi.m e gauss.m possiamo fare alcuni esperimenti.

1. Si osservi che per ottenere i nodi e pesi di Gauss-Legendre in (a, b) da quelli di Gauss-Jacobi in (−1, 1) abbiamo effettuato uno scaling: se x_{i,[−1,1]} sono i nodi di Gauss-Jacobi in [−1, 1] allora i nodi di Gauss-Jacobi x_{i,[a,b]} in [a, b] sono

x_{i,[a,b]} = ((a + b)/2) + ((b − a)/2) · x_{i,[−1,−1]};

2. se w_{i,[−1,1]} sono i pesi di Gauss-Jacobi in [−1, 1] allora i pesi di Gauss-Jacobi w_{i,[a,b]} in [a, b] sono

w_{i,[a,b]} = ((b − a)/2) · w_{i,[−1,1]}.

L’idea è la seguente. Dovendo le formule essere esatte per le costanti come la funzione f ≡ 1, nel caso della funzione peso di Legendre w ≡ 1

∑_{i} w_{i,[a,b]} = ∫_{a}^{b} w(x) dx = b − a

mentre nel caso di a = −1, b = 1 derivante dalla funzione peso di Jacobi abbiamo

∑_{i} w_{i,[−1,1]} = ∫_{−1}^{1} w(x) dx = 2;

si ha quindi l’intuizione che i pesi in (a, b) siano quelli in (−1, 1) moltiplicati per \frac{b−a}{2}.

3. Nonostante l’introduzione riguardante nodi e pesi in [a, b] non necessariamente uguale a [−1, 1], alla fine eseguiremo test esclusivamente in quest’ultimo intervallo. Qualora necessario, basta aggiungere nuove funzioni matematiche al file f.m, modificare adeguatamente a, b, fornire il relativo risultato esatto in exact results.m e testare la formula gaussiana. Per avere il risultato con alta precisione si usi la fun-

zione quadl.m o quad8.m di Matlab oppure integrazione gauss jacobi.m

con N = 500;

4. l’intervallo tipico di Gauss-Legendre è [a, b] (chiuso), mentre per Gauss-Jacobi l’intervallo è tipicamente (a, b) (aperto),poichè in generale gli esponenti della funzione peso di Jacobi possono essere negativi; per Gauss-Legendre il problema non sussiste, visto che la funzione peso è continua, e i punti a, b sono un insieme trascurabile.

Se ora testiamo il primo esempio, cioè il calcolo di

∫_{−1}^{1} x^{20} dx ≈ 0.09523809523809523300

otteniamo

>> format long;
>> N=11; ajac=0; bjac=0; a=-1; b=1;
>> [I_jac,x_jac,w_jac]=integrazione_gauss_jacobi(N,ajac,bjac,a,b,@f);
>> I_jac
I_jac =
0.09523809523810
>> length(x_jac)
ans =
 11
>> fprintf('\n \t [GAUSS-LEGENDRE]: \%15.20f',I_jac);
        [GAUSS-LEGENDRE]: 0.09523809523809564900
>> 0.09523809523809523300-0.09523809523809564900
ans =
-4.163336342344337e-016
>>

8.2.2. Un’integrale con funzione peso. Consideriamo ora l’integrale

∫_{−1}^{1} exp (x)\sqrt{1 − x} dx = 1.7791436546919097925911790299941.  (8.11)

Tale risultato è stato ottenuto usando il comando simbolico di Matlab 6.5 (non funziona in Octave, vedere in alternativa il programma Maxima!!)

>> syms x
>> int('(exp(x)) * ( (1-x)^(0.5) )',-1,1)
ans =
1.7791436546919097925911790299941

Si capisce che

1. syms x rende la variabile x di tipo simbolico (e non numerico!);

2. il termine

    int('(exp(x)) * ( (1-x)^(0.5) )',-1,1)

calcola simbolicamente l’integrale

∫_{−1}^{1} exp (x)\sqrt{1 − x} dx.

E’ immediato osservare che w(x) = \sqrt{1 − x} è un peso di Gauss-Jacobi

w(t) = (1 − t)^{α}(1 + t)^{β}

per α = 1/2 e β = 0.

Figura 8.4. Grafico che illustra l’errore delle formule di quadratura dei trapezi composta, Cavalieri-Simpson composta, Gauss-Legendre e Clenshaw-Curtis [11] sulla funzione test (8.11).

Infatti se w(x) = \sqrt{1 − x}, allora g(x) = exp (x) w(x) il che corrisponde a usare Gauss-Jacobi con f(x) = exp (x). Quindi paragoniamo le formule gaussiane (il codice funziona in Matlab 6.1 come pure nella release 2.1.73 di Octave

>> a=-1; b=1;
>> [I_jac,x_jac,w_jac]=integrazione_gauss_jacobi(10,1/2,0,a,b,@exp);
>> fprintf('\n \t [GAUSS-JACOBI]: %15.20f',I_jac);
           [GAUSS-JACOBI]: 1.77914365469190930000
>> 1.7791436546919097925911790299941-1.77914365469190930000
ans =
     4.440892098500626e-016
>> length(x_jac)
ans =
     10
>> [I_jac,x_jac,w_jac]=...
 integrazione_gauss_jacobi(10,0,0,-1,1,inline('exp(x).*sqrt(1-x)'));
>> fprintf('\n \t [GAUSS-LEGENDRE]: %15.20f',I_jac)
           [GAUSS-LEGENDRE]: 1.77984112101478020000
>> 1.7791436546919097925911790299941-1.77984112101478020000
ans =
 -6.9747e-004
>> length(x_jac)
ans =
     10
>>

Entrambe le formule hanno lo stesso numero di nodi (e pesi), come si vede dalla riga di comando

>> length(x_jac)
ans =
     10

ma offrono risultati diversi, con un errore assoluto di circa 4.44 · 10^{−16} per Gauss-Jacobi con a = 1/2, b = 0 e di 2.52 · 10^{−3} per Gauss-Legendre (cioè Gauss-Jacobi con a = 0, b = 0).

8.3. Esercizio. Si calcolino per N = 10, 20 con la formula composta dei trapezi, di Cavalieri-Simpson e un’appropriata formula gaussiana i seguenti integrali

∫_{−1}^{1} x^{20} dx = 2/21 ≈ 0.095238095238096801  (8.12)
∫_{−1}^{1} e^{x} dx = e − e^{−1} ≈ 2.3504023872876032  (8.13)
∫_{−1}^{1} e^{−x^{2}} dx = erf(1)·\sqrt{π} ≈ 1.4936482656248538  (8.14)
∫_{−1}^{1} 1/(1 + 16x^{2}) dx = 1/2 · atan(4) ≈ 0.66290883183401628  (8.15)
∫_{−1}^{1} e^{−1/x^{2}} dx ≈ 0.17814771178156086  (8.16)
∫_{−1}^{1} |x|^{3} dx = 1/2  (8.17)
∫_{0}^{1} \sqrt{x} dx = 2/3  (8.18)
∫_{−1}^{1} e^{x}\sqrt{1 − x} dx ≈ 1.7791436546919095  (8.19)

1. Quali delle due formule ha errori relativi inferiori?

2. Quali funzioni risultano più difficili da integrare numericamente?

3. La scelta N = 10, 20 nei singoli codici a quanti nodi corrisponde?

4. Suggerimento: ricordarsi che il peso w(x) ≡ 1 corrisponde ad un’opportuna scelta del peso di Jacobi.

9. Facoltativo: Equazioni di Fredholm di seconda specie. Si consideri l’equazione integrale di Fredholm di seconda specie

λu(x) = ∫_{a}^{b} k(x, y)u(y)dy + f(x), x ∈ [a, b]  (9.1)

dove k ∈ C([a, b] × [a, b]) e f ∈ C([a, b]), λ ∈ ℝ.

Esistono vari teoremi di esistenza ed unicità della soluzione u^{∗}. Un primo esempio è che sia (cf.[2, p.139])

max_{a≤x≤b} ∫ |k(x, y)|dy < |λ|.

Per un’analisi più accurata si consulti [12].

Introduciamo ora il metodo di Nyström (cf.[2, p.373]). Sia una sequenza di formule di quadratura

∫_{a}^{b} g(y)dy ≈ ∑_{j=1}^{n} w_{j,n}g(x_{j,n}), g ∈ C([a, b])  (9.2)

Figura 9.1. Erik Ivar Fredholm (1866-1927).

e supponiamo che tale sequenza converga all’integrale esatto qualsiasi sia g ∈ C([a, b]). Possiamo così considerare quale appossimazione del problema originario, la determinazione della funzione u_{n} soluzione del problema

λu_{n}(x) = ∑_{j=1}^{n} w_{j,n}k(x, x_{j,n})u_{n}(x_{j,n}) + f(x), x ∈ [a, b].  (9.3)

dove {u_{n}(x_{j,n})}_{j=1,...,q_{n}} verifica

λu(x_{i,n}) = ∑_{j=1}^{n} w_{j,n}k(x_{i,n}, x_{j,n})u_{n}(x_{j,n}) + f(x), x ∈ [a, b].  (9.4)

Posti A_{i,j}^{(n)} = w_{j,n}k(x_{i,n}, x_{j,n}), u_{n} = {u_{n}(x_{j,n})}_{j}, f_{n} = {f(x_{j,n})}_{j} si ha così che

λu_{n} = A^{(n)}u_{n} + f_{n}

per cui se B^{(n)} = λ · I − A^{(n)} è non singolare, allora

B^{(n)}u_{n} = f_{n},

sistema lineare che può essere risolto con il comando \ di Matlab.

Una volta ottenuto u_{n}(x_{j,n}), tramite (9.3), la funzione soluzione u_{n} può essere valutata in qualsiasi punto x ∈ [a, b]. Per le proprietà di convergenza di questo metodo si veda [2, p.376]

9.1. Facoltativo. Esercizio. Si consideri il problema

2 · u(x) = ∫_{0}^{1} exp (xy)u(y)dy + (2 exp(x) − \frac{1}{x + 1} · (exp(x + 1) − 1)), x ∈ [a, b]. (9.5)

avente unica soluzione exp (x) (perchè ?). Usando una formula di Gauss-Legendre, scalata nell’intervallo [0, 1], si valutino per n = 2, 4, 8, 16, 32 le soluzioni approssimate u_{n} nei nodi test x_{k} = kh dove h = 1/100, k = 0, …, 100 e si calcoli

E_{n} = max_{k=0,...,100} |u_{n}(x_{k}) − exp (x_{k})|.

Come diminuisce l’errore (farne un plot in scala semilogaritmica)? Suggerimento: in una versione di base, si può descrivere la valutazione dell’operatore integrale tramite due cicli for innestati

for i_index=1:length(x_pts)
for j_index=1:length(x_gl_ab)
 kern_eval=feval(kern,x_pts(i_index),x_gl_ab(j_index));
 B_gl_grid(i_index,j_index)=kern_eval*w_gl_ab(j_index);
end
end

10. Online. Quale nota storica osserviamo che da [22] Simpson is best remembered for his work on interpolation and numerical methods of integration. However the numerical method known today as Simpson’s rule, although it did appear in his work, was in fact due to Newton as Simpson himself acknowledged. By way of compensation, however, the Newton-Raphson method for solving the equation f(x) = 0 is, in its present form, due to Simpson. Newton described an algebraic process for solving polynomial equations which Raphson later improved. The method of approximating the roots did not use the differential calculus. The modern iterative form x_{n+1} = x_{n} − f(x_{n})/f′(x_{n}) is due to Simpson, who published it in 1740.

11. Frasi celebri.

1. Cayley, Hermite and I constitute the Invariantive Trinity. (Sylvester)

2. The traditional mathematics professor of the popular legend is absentminded. He usually appears in public with a lost umbrella in each hand. He prefers to face a blackboard and to turn his back on the class. He writes a, he says b, he means c, but it should be d. Some of his sayings are handed down from generation to generation:

(a) In order to solve this differential equation you look at it till a solution occurs to you.

(b) This principle is so perfectly general that no particular application of it is possible.

(c) Geometry is the science of correct reasoning on incorrect figures.

(d) My method to overcome a difficulty is to go round it.

(e) What is the difference between method and device? A method is a device

which you used twice.

3. Even fairly good students, when they have obtained the solution of the problem and written down neatly the argument, shut their books and look for something else. Doing so, they miss an important and instructive phase of the work. ... A good teacher should understand and impress on his students the view that no problem whatever is completely exhausted. (Polya)

4. One of the first and foremost duties of the teacher is not to give his students the impression that mathematical problems have little connection with each other, and no connection at all with anything else. We have a natural opportunity to investigate the connections of a problem when looking back at its solution. (Polya)

5. If there is a problem you can’t solve, then there is an easier problem you can solve: find it. (Polya)

6. A great discovery solves a great problem, but there is a grain of discovery in the solution of any problem. Your problem may be modest, but if it challenges your curiosity and brings into play your inventive faculties, and if you solve it by your own means, you may experience the tension and enjoy the triumph of discovery. (Polya)

7. The first rule of discovery is to have brains and good luck. The second rule of discovery is to sit tight and wait till you get a bright idea. (Polya)

8. If you have to prove a theorem, do not rush. First of all, understand fully what the theorem says, try to see clearly what it means. Then check the theorem; it could be false. Examine the consequences, verify as many particular instances as are needed to convince yourself of the truth. When you have satisfied yourself that the theorem is true, you can start proving it. (Polya)

9. Mathematics consists of proving the most obvious thing in the least obvious way. (Polya)

10. Mathematics is the cheapest science. Unlike physics or chemistry, it does not require any expensive equipment. All one needs for mathematics is a pencil and paper. (Polya)

11. There are many questions which fools can ask that wise men cannot answer. (Polya)

12. A mathematician who can only generalise is like a monkey who can only climb up a tree, and a mathematician who can only specialise is like a monkey who can only climb down a tree. In fact neither the up monkey nor the down monkey is a viable creature. A real monkey must find food and escape his enemies and so must be able to incessantly climb up and down. A real mathematician must be able to generalise and specialise. (Polya)

13. Look around when you have got your first mushroom or made your first discovery: they grow in clusters. (Polya)

14. John von Neumann was the only student I was ever afraid of. (Polya)

15. The apex and culmination of modern mathematics is a theorem so perfectly general that no particular application of it is feasible. (Polya)

16. I am too good for philosophy and not good enough for physics. Mathematics is in between. (Polya)

17. The elegance of a mathematical theorem is directly proportional to the number of independent ideas one can see in the theorem and inversely proportional to the effort it takes to see them. (Polya)

RIFERIMENTI BIBLIOGRAFICI

[1] K. Atkinson, Introduction to Numerical Analysis, Wiley, 1989.

[2] K. Atkinson e W. Han, Theoretical Numerical Analysis, Springer, 2001.

[3] V. Comincioli, Analisi Numerica, metodi modelli applicazioni, Mc Graw-Hill, 1990.

[4] S.D. Conte e C. de Boor, Elementary Numerical Analysis, 3rd Edition, Mc Graw-Hill, 1980.

[5] P.J. Davis, Interpolation and Approximation, Dover, 1975.

[6] W. Gautschi, personal homepage, http://www.cs.purdue.edu/people/wxg .

[7] W. Gautschi, , Orthogonal polynomials (in Matlab), Journal of Computational and Applied Mathematics 178 (2005) 215234

[8] The MathWorks Inc., Numerical Computing with Matlab, http://www.mathworks.com/moler .

[9] A. Quarteroni e F. Saleri, Introduzione al calcolo scientifico, Springer Verlag, 2006.

[10] A. Suli e D. Mayers, An Introduction to Numerical Analysis, Cambridge University Press, 2003.

[11] L.N. Trefethen, Is Gauss quadrature better than Clenshaw-Curtis?, SIAM Reviews, to appear (2007), http://web.comlab.ox.ac.uk/oucl/work/nick.trefethen/CC trefethen revised2.pdf .

[12] F.G. Tricomi, Integral Equations, Dover, 1985.

[13] J. Waldvogel, Fast Construction of the Fejer and Clenshaw-Curtis Quadrature Rules, BIT 46 (2006), no. 1, 195–202, http://www.math.ethz.ch/w̃aldvoge/Papers/fejer.pdf .

[14] Mac Tutor, Cavalieri, http://www-gap.dcs.st-and.ac.uk/ history/Biographies/Cavalieri.html .

[15] Mac Tutor, Cotes,

http://www-history.mcs.st-andrews.ac.uk/Biographies/Cotes.html .

[16] Mac Tutor, Fredholm, http://www-history.mcs.st-andrews.ac.uk/Biographies/Fredholm.html .

[17] Mac Tutor, Hermite, http://www-history.mcs.st-andrews.ac.uk/Biographies/Hermite.html .

[18] Mac Tutor, Jacobi, http://www-history.mcs.st-andrews.ac.uk/Biographies/Jacobi.html .

[19] Mac Tutor, Laguerre, http://www-history.mcs.st-andrews.ac.uk/Biographies/Laguerre.html .

[20] Mac Tutor, Legendre, http://www-history.mcs.st-andrews.ac.uk/Biographies/Legendre.html .

[21] Mac Tutor, Polya, http://www-history.mcs.st-andrews.ac.uk/Biographies/Polya.html .

[22] Mac Tutor, Simpson, http://www-history.mcs.st-andrews.ac.uk/Biographies/Simpson.html .

[23] Mac Tutor, Steklov, http://www-history.mcs.st-andrews.ac.uk/Biographies/Steklov.html .

[24] Mac Tutor, Stieltjes, http://www-history.mcs.st-andrews.ac.uk/Biographies/Stieltjes.html .

[25] Wikipedia, Cavalieri, http://it.wikipedia.org/wiki/Bonaventura Cavalieri .

[26] Wikipedia, Clenshaw-Curtis quadrature, http://en.wikipedia.org/wiki/Clenshaw-Curtis quadrature .

[27] Wikipedia, Fast Fourier Transform, http://en.wikipedia.org/wiki/Fast Fourier transform .

[28] Wikipedia, Gaussian quadrature, http://en.wikipedia.org/wiki/Gaussian quadrature .

[29] Wikipedia, Quadratura gaussiana, http://it.wikipedia.org/wiki/Quadratura di Gauss .

[30] Wikipedia, Regola del trapezio, http://it.wikipedia.org/wiki/Regola del trapezio .

[31] Wikipedia, Regola di Cavalieri-Simpson, http://it.wikipedia.org/wiki/Regola di Cavalieri-Simpson .