Quadratura numerica

Alvise Sommariva

Università degli Studi di Padova

Dipartimento di Matematica

Quadratura numerica

4 aprile 2014 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)

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

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

∫_{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)

per cui, confrontando con la formula (1) abbiamo

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

Alvise Sommariva  Quadratura numerica  3/ 116

Quadratura numerica

Definizione

Una formula di quadratura

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

per cui

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

si dice interpolatoria.

Alvise Sommariva  Quadratura numerica  4/ 116

Definizione

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.

Quadratura numerica

Mostriamo ora il seguente

Teorema

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.

Alvise Sommariva  Quadratura numerica  6/ 116

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}.

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. [8, 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.

Alvise Sommariva  Quadratura numerica  8/ 116

Formule di Newton-Cotes

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

Alvise Sommariva  Quadratura numerica  9/ 116

Alcune classiche regole sono: ▶ 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 [14] è

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

per qualche ξ ∈ (a, b); ▶ 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 [15] è

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

per qualche ξ ∈ (a, b); 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}

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}  (5)

Formule di Newton-Cotes

Inoltre

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}  (6)
Alvise Sommariva  Quadratura numerica  13/ 116

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], [8, p.333-336].

▶ 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]).

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}) (7)

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

Formule di Newton-Cotes composte

Due casi particolari sono

▶ 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}] (8)

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);

Alvise Sommariva  Quadratura numerica  17/ 116

Implementazione Matlab di alcune formule composte

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

 f u n c t i o n [ 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 ;

Implementazione Matlab di alcune formule composte

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})  (9)

come pure per (8)

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

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.

Implementazione Matlab di alcune formule composte

Si potrebbe usare il comando Matlab trapz nella sua implementazione

>> h e l p t r a p z

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
f o r 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∗ t r a p z ( fx ) ;
Alvise Sommariva  Quadratura numerica  20/ 116

Implementazione Matlab di alcune formule composte

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 ( 1 0 , 0 , 1 ) ;

>> fx=s i n ( x ) ;

>> I_trapezi_composta=w ’ ∗ fx
I_trapezi_composta =

0.45931454885798

>> h=(1−0) /10;
>> I_trapz=h∗ t r a p z ( fx )

I_trapz =

0.45931454885798 >>

Alvise Sommariva  Quadratura numerica  21/ 116

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).

Implementazione Matlab di alcune formule composte

Per quanto riguarda la formula di Cavalieri-Simpson composta

 f u n c t i o n [ 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 ( l e n g t h ( 3 : 2 : 2 ∗ N−1) , 1 ) ;
w ( 2 : 2 : 2 ∗ N , 1 ) =4∗ones ( l e n g t h ( 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 (9) e (11) 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).

Implementazione Matlab di alcune formule composte

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

>> format long ;

>> [ x , w]= simpson_composta ( 5 , 0 , 1 ) ;

>> fx=s i n ( x ) ;

>> I_simpson=w ’ ∗ fx
I_simpson =

0.45969794982382

>> l e n g t h ( x )
ans =

11 >>

Implementazione Matlab di alcune formule composte

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.

a=−1; b=1; N=11;  %SCEGLIERE N DISPARI .

N_trap=N−1;
[ x_trap , w_trap]= trapezi_composta ( N_trap , a , b ) ;

fx_trap=x_trap . ˆ 2 0 ;  % 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 . ˆ 2 0 ;  % VALUT. FUNZIONE .
I_simp=w_simp ’ ∗ fx_simp ;  % SIMPSON COMPOSTA.
f p r i n t f ( ’ \n\ t [TPZ COMP] [ PTS] : % 4 . 0 f ’ , l e n g t h ( x_trap ) ) ;
f p r i n t f ( ’ \n\ t [TPZ COMP] [ RIS ]:%14.14 f ’ , I_trap ) ;
f p r i n t f ( ’ \n\ t [ SIMPS .COMP] [ PTS] : % 4 . 0 f ’ , l e n g t h ( x_simp ) ) ;
f p r i n t f ( ’ \n\ t [ SIMPS .COMP] [ RIS ]:%14.14 f ’ , I_simp ) ;
Alvise Sommariva  Quadratura numerica  27/ 116

Implementazione Matlab di alcune formule composte

ottenendo

[ TPZ COMP ] [ PTS ] :  11

[ TPZ COMP ] [ RIS ] : 0.20462631505024

[ SIMPS . COMP ] [ PTS ] :  11

[ SIMPS . COMP ] [ RIS ] : 0.13949200364447

Implementazione Matlab di alcune formule gaussiane

Si può vedere che usando formule di tipo gaussiano (cf. [13]) o di tipo Clenshaw-Curtis (cf. [12], [10], [11]) 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.

Implementazione Matlab di alcune formule composte

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

Alvise Sommariva  Quadratura numerica  30/ 116

Formule di Newton-Cotes composte

▶ 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})]; (11)

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).

Alvise Sommariva  Quadratura numerica  31/ 116

Esercizio

Usando la formula composta dei trapezi e di Cavalieri-Simpson, con n = 2, 4, 8, …, 512 suddivisioni equispaziate dell’intervallo di integrazione, calcolare

▶ ∫_{0}^{π} exp (x) · cos (4x)dx = \frac{exp (π)−1}{17};

▶ ∫_{0}^{1} x^{5/2}dx = \frac{2}{7};

▶ ∫_{−π}^{π} exp cos (x)dx = 7.95492652101284;

▶ ∫_{0}^{π/4} exp cos (x)dx = 1.93973485062365;

▶ ∫_{0}^{1} x^{1/2}dx = \frac{2}{3}. Descrivere come decresce l’errore in scala semilogaritmica e quale sia l’andamento del rapporto E_{n}(f )/E_{n+1}(f ), dove E_{n} = |I (f ) − I_{n}(f )| con I (f ) integrale esatto e I_{n}(f ) valore fornito dalla formula di quadratura.

Esercizio facoltativo (non facile)

Aiutandosi con quanto fornito online relativamente al metodo di Romberg, fornirne una implementazione e testarla sulle funzioni, utilizzando n = 2, 4, 8, …, 512 suddivisioni equispaziate dell’intervallo di integrazione, calcolando

▶ ∫_{0}^{π} exp (x) · cos (4x)dx = \frac{exp (π)−1}{17};

▶ ∫_{0}^{1} x^{5/2}dx = \frac{2}{7};

▶ ∫_{−π}^{π} exp cos (x)dx = 7.95492652101284;

▶ ∫_{0}^{π/4} exp cos (x)dx = 1.93973485062365;

▶ ∫_{0}^{1} x^{1/2}dx = \frac{2}{3}. Lo si può applicare in tutti i casi considerati?

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.

Consideriamo 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 (−∞, ∞); 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 nodi 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 δ = 2n − 1, cioè calcolino esattamente

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

per ogni polinomio p il cui grado è minore o uguale a 2n − 1? La risposta è affermativa, come si può vedere in [1, p.272].

Teorema

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. Siano 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}.

▶ ∫_{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

Dimostrazione.

▶ 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}).

Formule gaussiane

Dimostrazione.

Quindi, 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}) (12)

per cui la formula di Gauss ha grado di precisione almeno 2n − 1.

Alvise Sommariva  Quadratura numerica  42/ 116

Dimostrazione.

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.

Dimostrazione.

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}.

Dimostrazione.

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})  (13)

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.

Sull’errore di quadratura delle formule di Newton-Cotes e di Gauss

Riguardo gli errori compiuti da alcune formule di quadratura (cf.[1, p.264])

Teorema

Sia la regola di Newton-Cotes I (f ) ≈ I_{n}(f ) = ∑_{i=0}^{n} w_{i,n}f (x_{i,n}).

▶ 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μ;

▶ 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).

Alvise Sommariva  Quadratura numerica  46/ 116

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

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 ) − I_{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.

Alvise Sommariva  Quadratura numerica  47/ 116

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

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).

Alvise Sommariva  Quadratura numerica  48/ 116

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

(14)

e che invece di {f_{j}}_{j} si disponga di una loro approssimazione {f̃_{j}}_{j}.

Stabilità di una formula di quadratura

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}|  (15)

Quindi la quantità

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

è un indice di stabilità della formula di quadratura.

Alvise Sommariva  Quadratura numerica  50/ 116

Supponiamo che alcuni pesi siano negativi e che la formula abbia almeno grado di precisione 0. 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}^{−}  (16)

∑_{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}^{−}(17)

da cui essendo w_{l}^{−} < 0

∑_{k=1}^{η_{n}}|w_{k}| = ∑_{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}^{−} ≥ ∫_{a}^{b} w(x)dx (18)

Stabilità di una formula di quadratura

mentre se i pesi fossero tutti positivi avremmo

∑_{k=1}^{η_{n}}|w_{k}| = ∑_{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.

Alvise Sommariva  Quadratura numerica  53/ 116

Alcune norme di operatori

Se (a, b) è limitato allora 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 ∥_{∞} (19)

Alcune norme di operatori

In particolare, scegliendo opportunamente f si prova che 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}|. Inoltre se (a, b) è limitato

∥I ∥_{∞} = max_{f ∈C([a,b]),f ≠0} \frac{|I (f )|}{∥f ∥_{∞}} = ∫_{a}^{b} w(x)dx = ∥w∥_{∞} in quanto

|I (f )| = |∫_{a}^{b} f (x)w(x)dx| ≤ ∫_{a}^{b} |f (x)|w(x)dx
≤ ∫_{a}^{b} w(x)dx · ∥f ∥_{∞} = ∥w∥_{∞}∥f ∥_{∞}  (20)

e

|I (1)| = ∥w∥_{∞}.
Alvise Sommariva  Quadratura numerica  55/ 116

Teorema di Stieltjes: legame tra errore di quadratura e miglior approssimazione polinomiale

Mostriamo il seguente teorema (di Stieltjes), che lega l’errore delle formule di quadratura a quello di miglior approssimazione. Con P_{n} denotiamo i polinomi di grado al più n.

Teorema

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 posto E_{n}(f ) = I (f ) − I_{n}(f ) si ha

|E_{n}(f )| ≤ (∥w∥_{1} + ∑_{j=1}^{η_{n}}|w_{j}|) · min_{q_{n}∈P_{n}} ∥f − q_{n}∥_{∞}.  (21)

Dimostrazione.

Se q ∈ P_{n} è un polinomio arbitrario di grado n ed I (q) = I_{n}(q) avendo la formula di quadratura grado di precisione almeno n. Ricordiamo inoltre che gli operatori I ed I_{n} sono lineari e quindi

I_{n}(f − q) = I_{n}(f ) − I_{n}(q),
I (f − q) = I (f ) − I (q).

Dimostrazione.

Per quanto visto

|I_{n}(f )| ≤ ∥I_{n}∥_{∞}∥f ∥_{∞}, |I (f )| ≤ ∥f ∥_{∞}∥w∥_{1}

poichè per ogni polinomio di grado n si ha che I_{n}(q) = I (q),

|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∥_{∞}
= (∥w∥_{1} + ∑_{j=1}^{η_{n}}|w_{j}|) · ∥f − q∥_{∞}  (22)

L’interesse di questo teorema è il legame col polinomio di miglior approssimazione. Risulta importante osservare che in (21) contribuiscono i prodotti di due termini.

1. Il primo è dovuto alla funzione peso e alla stabilità della formula di quadratura.

2. 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, come suggerito dai teoremi di Jackson sulla miglior approssimante polinomiale di una funzione f , che forniscono stime di

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

con f ∈ C([a, b]) (dotando C([a, b])) della norma infinito).

Quale esempio consideriamo una formula a pesi positivi e grado di precisione n ≥ 0. Necessariamente ∥I_{n}∥_{∞} = ∥w∥_{1} = ∫_{a}^{b} w(x)dx, in quanto la formula integra esattamente la costante 1 ed è

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

Quindi, 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}∥_{∞}.

Esercizio. Si calcoli 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è ?

Teorema di Polya-Steklov, sulla convergenza delle formule di quadratura

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})  (23)

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}).

Risulta importante sottolineare che n non è il grado di precisione della formula, bensì l’indice che contraddistingue l’n-simo elemento della sequenza di formule {I_{n}}_{n} = 1, …,.

Alvise Sommariva  Quadratura numerica  63/ 116

Dimostriamo ora il teorema di Polya-Steklov [3, p.202].

Teorema

Sia un intervallo compatto [a, b], e Q_{n}(f ) = ∑_{i=0}^{η_{n}} w_{i,n}f (x_{i,n}) una sequenza di formule di quadratura tale che I_{w}(f ) ≈ Q_{n}(f ). Condizione necessaria e sufficiente affinchè per ogni f ∈ C([a, b])

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

con E_{n}(f ) := I_{w}(f ) − Q_{n}(f ) è che

1. esista M ∈ ℝ indipendente da n per cui si abbia

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

2. per ogni k ∈ ℕ si abbia

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

Alvise Sommariva  Quadratura numerica  64/ 116

Dimostrazione.

⇐ Supponiamo che

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

∑_{i=1}^{η_{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 τ_{1} > 0 esiste un polinomio p tale che ∥f − p∥_{∞} ≤ τ_{1}.

Fissato n, ricordato che per la definizione di norma degli operatori, si ha che

∥I ∥_{∞} = sup_{g∈C([a,b]),g≠0} \frac{|I (g)|}{∥g∥_{∞}}

implica, dato che ∥I ∥_{∞} = ∥w∥_{1}, per ogni g ∈ C([a, b])

|I (g)| ≤ ∥I ∥_{∞}∥g∥_{∞} = ∥w∥_{1}∥g∥_{∞}  (24)

Similmente

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

implica, ricordando che ∥I_{n}∥_{∞} = ∑_{i=1}^{η_{n}} |w_{i,n}|,

|I_{n}(g)| ≤ ∥I_{n}∥_{∞}∥g∥_{∞} = ∑_{i=1}^{η_{n}}|w_{i,n}|∥g∥_{∞}.  (25)

Quindi, per g = f − p in (24) e (25)

|E_{n}(f − p)| = |I (f − p) − I_{n}(f − p)| ≤ |I (f − p)| + |I_{n}(f − p)|
≤ ∥w∥_{1}∥f − p∥_{∞} + ∑_{i=1}^{η_{n}}|w_{i,n}|∥f − p∥_{∞}

= (∥w∥_{1} + ∑_{i=1}^{η_{n}}|w_{i,n}|) · ∥f − p∥_{∞}
≤ (∥w∥_{1} + M) · τ_{1}.  (26)

Di conseguenza |E_{n}(f − p)| ≤ (∥w∥_{1} + M) · τ_{1}, per ogni n ∈ ℕ. Si osservi che il secondo membro della precedente disuguaglianza non dipende da n.

Poichè lim_{n→∞} E_{n}(p) = 0, per ogni τ_{2} > 0 esiste Ñ (τ_{2}) ∈ ℕ tale che se n ≥ Ñ (τ_{2}) allora |E_{n}(p)| ≤ τ_{2}. Di conseguenza fissato ϵ > 0 e posto

▶ τ_{1} = ϵ/(2 · (∥w∥_{1} + M)), ▶ τ_{2} = ϵ/2,

abbiamo che esiste N(ϵ) = Ñ (ϵ/2) tale che se n ≥ N(ϵ) allora

|E_{n}(f − p)| ≤ (∥w∥_{1} + M) · τ_{1} = \frac{(∥w∥_{1} + M) · ϵ}{2 · (∥w∥_{1} + M)} = ϵ/2

e

|E_{n}(p)| ≤ τ_{2} = ϵ/2.

Alvise Sommariva  Quadratura numerica  69/ 116

Quindi, dalla linearità di E_{n}, per n ≥ N(ϵ) = Ñ (τ_{2})

|E_{n}(f )| ≤ |E_{n}(f ) − E_{n}(p)| + |E_{n}(p)| = |E_{n}(f − p)| + |E_{n}(p)|
≤ \frac{ϵ}{2} + \frac{ϵ}{2} = ϵ  (27)

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

⇒ Mostriamo il primo punto dell’asserto. Supponiamo che per ogni f ∈ C([a, b]) sia lim_{n} E_{n}(f ) = 0. Per definizione di E_{n}(f ) abbiamo I_{n}(f ) = I (f ) − E_{n}(f ) e quindi per la disuguaglianza triangolare

|I_{n}(f )| ≤ |I (f )| + |E_{n}(f )| ≤ ∥w∥_{1}∥f ∥_{∞} + |E_{n}(f )|.

Poichè lim_{n} E_{n}(f ) = 0 necessariamente lim_{n} |E_{n}(f )| = 0 e quindi, dalla definizione di limite, segue facilmente che esiste M(f ) ∈ ℝ (indipendente da n, ma dipendente da f ) tale che I_{n}(f ) ≤ M(f ) < ∞.

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 , ▶ 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} ≡ I_{n}, operatore lineare limitato con norma

∥I_{n}∥_{∞} = ∑_{i=0}^{η_{n}} |w_{i,n}|, se f ∈ C([a, b]) abbiamo che la sequenza {I_{n}(f )}_{n} è limitata in quanto esiste M(f ) < ∞ indipendente da n tale che I_{n}(f ) ≤ M(f ) < ∞.

Conseguentemente per il teorema di Banach-Steinhaus, ricordando che

∥I_{n}∥_{∞} = ∑_{i=0}^{η_{n}}|w_{i,n}|

si ha

sup_{n}(∑_{i=0}^{η_{n}}|w_{i,n}|) = sup_{n}∥I_{n}∥_{∞} < +∞.

e quindi esiste M finito tale che

∑_{i=0}^{η_{n}}|w_{i,n}| ≤ M < +∞, ∀n ∈ ℕ.

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

Alvise Sommariva  Quadratura numerica  74/ 116

Alcune conseguenze e osservazioni sul Teorema di Polya-Steklov

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

▶ Si osservi che in generale le formule di errore introdotte nei capitoli precedenti, implicavano la convergenza in caso l’integranda f fosse sufficientemente regolare. Nel teorema di Polya-Steklov si chiede esclusivamente che f ∈ C([a, b]), senza però offrire stime dell’errore compiuto.

Teorema

Consideriamo una formula su un dominio limitato, con i pesi w_{i,n} positivi, e g.d.p. almeno 0. Essa è convergente per ogni f ∈ C([a, b]) se e solo se è convergente per ogni polinomio p.

Dimostrazione.

Infatti se è convergente per polinomi in particolare lo è per f (x) = 1, cioè lim_{n} E_{n}(1) = 0. Visto che

E_{n}(1) = ∫_{a}^{b} w(x)dx − ∑_{k=1}^{η_{n}}w_{k} → 0

necessariamente ∑_{k=1}^{η_{n}} w_{k} ≤ M < +∞.

Teorema

Una formula gaussiana su un dominio limitato è convergente.

Dimostrazione.

Sia una formula di Gauss, su un dominio limitato, con n nodi {w_{i,n}}_{i=1,...,n} positivi. Per quanto detto, è convergente per ogni f ∈ C([a, b]) se e solo se è convergente per ogni polinomio p. Ma ciò è verificato banalmente in quanto essendo il grado di precisione almeno 2n − 1, fissato k, per n ≥ ceil((k + 1)/2) si ha E_{n}(x^{k}) = 0.

Quindi, 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 della formula gaussiana si ha che

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

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

Alvise Sommariva  Quadratura numerica  77/ 116

Teorema

Una sequenza di formule composte, basate su regole a pesi positivi e g.d.p. n ≥ 0, risulta convergente qualora l’ampiezza delle suddivisioni tenda a 0.

Dimostrazione.

Si consideri la suddivisione Δ_{m} = {τ_{i}}_{i=0,...,m} dell’intervallo (a, b) con

τ_{i} < τ_{i+1}, τ_{0} = a, τ_{m} = b

nodi della formula in questione.

Visto che la formula composta ha pesi positivi, risulta convergente per ogni f ∈ C([a, b]) se e solo se è convergente per ogni polinomio p.

Alvise Sommariva  Quadratura numerica  78/ 116

Mostriamo che è convergente per qualsiasi polinomio p.

Osserviamo che tale formula composta integra esattamente ogni funzione polinomiale a tratti di grado n su Δ e che se s_{Δ_{m},n} è l’interpolante polinomiale a tratti di grado n della funzione f relativamente alla suddivisione Δ_{m} e ai nodi di quadratura,

∫_{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 (28)

da cui se f ∈ C^{(n+1)}([a, b])

|E_{n}(f )| = |∫_{a}^{b} f (x)dx − ∑_{i=1}^{η_{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)!}

(29)

Alvise Sommariva  Quadratura numerica  80/ 116

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 ampiezza h della suddivisione tende a 0 allora E_{n}(x^{k}) → 0. Per il teorema di Polya-Steklov si ha così che la formula composta a pesi positivi

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

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

Alvise Sommariva  Quadratura numerica  81/ 116

Alcune conseguenze e osservazioni sul Teorema di Polya-Steklov (facoltativo)

[Facoltativo] Una formula I_{n}(f ) ≈ ∫_{a}^{b} f (x)w(x)dx := I (f ) a pesi positivi convergente sui polinomi e avente grado di precisione almeno 0 risulta convergente sulle funzioni continue a tratti in [a, b]. Se f è tale funzione, si dimostra che per ogni ϵ > 0 esistono due funzioni f_{1}, f_{2} ∈ C([a, b]) tali che f_{1} ≤ f ≤ f_{2} e

∥f_{1} − f_{2}∥_{∞} ≤ ϵ.

Osserviamo che

▶ Poichè la formula è a pesi positivi, f_{1} ≤ f ≤ f_{2} implica che

I_{n}(f_{1}) ≤ I_{n}(f ) ≤ I_{n}(f_{2}). ▶ Essendo f_{1} ≤ f ≤ f_{2},

∫_{a}^{b} f_{1}(x)w(x)dx ≤ ∫_{a}^{b} f (x)w(x)dx ≤ ∫_{a}^{b} f_{2}(x)w(x)dx.

▶ Per il teorema di Polya-Steklov abbiamo inoltre che essendo

f_{1}, f_{2} ∈ C([a, b]) allora lim_{n} E_{n}(f_{1}) = 0 e lim_{n} E_{n}(f_{2}) = 0.

Alvise Sommariva  Quadratura numerica  82/ 116

Ora notiamo che

I (f_{1}) − I_{n}(f_{2}) ≤ I (f ) − I_{n}(f ) ≤ I (f_{2}) − I_{n}(f_{1}).  (31)

Per la linearità degli operatori I , I_{n}

E_{n}(f_{1})+I_{n}(f_{1}−f_{2}) = (I (f_{1})−I_{n}(f_{1}))+I_{n}(f_{1}−f_{2}) = I (f_{1})−I_{n}(f_{2}) (32)
E_{n}(f_{2})+I_{n}(f_{2}−f_{1}) = (I (f_{2})−I_{n}(f_{2}))+I_{n}(f_{2}−f_{1}) = I (f_{2})−I_{n}(f_{1}) (33)

Inoltre

I_{n}(f_{2} − f_{1}) = ∑_{k=1}^{η_{n}}w_{k}(f_{2}(x_{k}) − f_{1}(x_{k})) ≤ ∑_{k=1}^{η_{n}}w_{k}∥f_{2} − f_{1}∥_{∞} (34)
I_{n}(f_{1} − f_{2}) = −I_{n}(f_{2} − f_{1}) ≥ − ∑_{k=1}^{η_{n}}w_{k}∥f_{2} − f_{1}∥_{∞}  (35)

Alcune conseguenze e osservazioni sul Teorema di Polya-Steklov (facoltativo)

Quindi da (31), in virtù di (33), (34) E_{n}(f ) ≤ I (f_{2})−I_{n}(f_{1}) = E_{n}(f_{2})+I_{n}(f_{2}−f_{1}) ≤ E_{n}(f_{2})+∑_{k=1}^{η_{n}}w_{k}∥f_{2}−f_{1}∥_{∞} mentre da (31), in virtù di (32), (35) E_{n}(f_{1})−∑_{k=1}^{η_{n}}w_{k}∥f_{2}−f_{1}∥_{∞} ≤ E_{n}(f_{1})+I_{n}(f_{1}−f_{2}) = I (f_{1})−I_{n}(f_{2}) ≤ E_{n}(f ) cioè

E_{n}(f_{1}) − ∑_{k=1}^{η_{n}}w_{k}∥f_{2} − f_{1}∥_{∞} ≤ E_{n}(f ) ≤ E_{n}(f_{2}) + ∑_{k=1}^{η_{n}}w_{k}∥f_{2} − f_{1}∥_{∞}

Dal fatto che ∑_{k=1}^{η_{n}} w_{k} = ∫_{a}^{b} w(x)dx < +∞, lim_{n} E_{n}(f_{1}) = lim_{n} E_{n}(f_{2}) = 0 e ∥f_{2} − f_{1}∥_{∞} ≤ ϵ, dall’arbitrarietà di ϵ deduciamo che E_{n}(f ) → 0.

Alvise Sommariva  Quadratura numerica  85/ 116

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, genera i primi n coeff. di ricorrenza dei polinomi ortogonali \underline{monici} relativamente alla funzione peso di Jacobi

w(x) = (1 − x)^{a} (1 + x)^{b}.

Quindi

▶ 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.  (36)
Alvise Sommariva  Quadratura numerica  88/ 116

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}\}.

Formule gaussiane in Matlab

A questo punto si chiama la funzione gauss.m (reperibile nuovamente presso lo stesso sito di W. Gautschi)

f u n c t i o n xw=gauss ( N , ab )
N0=s i z e ( ab , 1 ) ; i f N0<N , e r r o r ( ’ i n p u t a r r a y ab too s h o r t
’ ) , end
J=z e r o s ( N ) ;
f o r n=1:N , J ( n , n )=ab ( n , 1 ) ; end
f o r n=2:N
J ( n , n−1)=s q r t ( ab ( n , 2 ) ) ;
J ( n−1,n )=J ( n , n−1) ;

end

[ V , D]= e i g ( J ) ;
[ D , I]= s o r t ( d i a g ( D ) ) ;
V=V ( : , I ) ;
xw=[D ab ( 1 , 2 ) ∗V ( 1 , : ) ’ . ˆ 2 ] ;

▶ Tale routine, che per un certo N, uguale al massimo al numero di righe della matrice di coefficienti di ricorrenza ab, fornisce N nodi x e i corrispettivi pesi w immagazzinati in una matrice xw che ha quale prima colonna x e quale seconda colonna w. Osserviamo che in tale codice 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.

Formule gaussiane in Matlab

Figura: Grafico che illustra la distribuzione dei 20 nodi di Gauss-Legendre nell’intervallo [1, 1]

.

Alvise Sommariva  Quadratura numerica  92/ 116

Formule gaussiane in Matlab

Figura: Grafico che illustra la distribuzione dei 20 pesi di Gauss-Legendre nell’intervallo [1, 1]

.

Alvise Sommariva  Quadratura numerica  93/ 116

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.

Una demo di Gauss-Jacobi

Salviamo nel file integrazione gauss jacobi.m

 f u n c t i o n [ I , x , w]= integrazione_gauss_jacobi ( N , alpha , beta

, a , b , f )
% INPUT :
% N: NUMERO NODI .
% alpha , beta : w( x ) =((1−x ) ˆ a l p h a ) ∗((1+x ) ˆ beta )
% a , b : SE a l p h a =0, beta =0, SONO ESTREMI DI INTEGRAZIONE

,

%  ALTRIMENTI a=−1, b=1.
% f : INTEGRANDA (RISPETTO w) .
% OUTPUT:
% I : APPROSSIMAZIONE DI \ i n t a ˆ b f ( x ) w( x ) dx
% x ,w: NODI/PESI UTILIZZATI .

Una demo di Gauss-Jacobi

Salviamo nel file gauss jacobi.m

ab=r_jacobi ( N , alpha , beta ) ; %TERM. RICORSIVI
xw=gauss ( N , ab ) ; %NODI/PESI IN MATRICE .
x=xw ( : , 1 ) ; w=xw ( : , 2 ) ;%NODI/PESI GAUSS JAC . [ − 1 , 1 ] .
i f ( alpha = 0) & ( beta = 0) & ( ( a˜ =−1) | ( b˜ =1) )
x=((a+b ) /2) +((b−a ) /2) ∗x ;%NODI GAUSS−LEG . [ a , b ] .
w=((b−a ) /2) ∗w ;%PESI GAUSS−LEG . [ a , b ] .

end

fx=f e v a l ( f , x ) ;%VALUTAZIONE FUNZIONE .
I=w ’ ∗ fx ;%VALORE INTEGRALE .

Una demo di Gauss-Jacobi

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

f u n c t i o n fx=f ( x )
fx=x . ˆ 2 0 ;
% ALCUNE FUNZIONI TRATTE DA
% ” IS GAUSS QUADRATURE BETTER THAN CLENSHAW−CURTIS?”
% DI L .N. TREFETHEN.
% f x=exp ( x ) ;
% f x=exp(−x . ˆ 2 ) ;
% f x =1./(1+16∗(x . ˆ 2 ) ) ;
% f x=exp(−x .ˆ ( −2) ) ;
% f x=abs ( x ) ; f x=f x . ˆ 3 ;
% f x=x . ˆ ( 0 . 5 ) ;
% f x=exp ( x ) . ∗ ( s q r t (1−x ) ) ;

Una demo di Gauss-Jacobi

Le funzioni test [10]

1. fx = x.^{2}0;

2. fx = exp(x);

3. fx = exp(−x.^{2});

4. fx = 1./(1 + 16 ∗ (x.^{2}));

5. fx = exp(−x.^{(} − 2));

6. fx = (abs(x)).^{3};

7. fx = x.^{(}0.5);

8. fx = exp(x). ∗ (sqrt(1 − x)); con a = −1, b = 1 possono essere valutate a scelta cencellando % riponenolo altrove. 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 quindi alcuni esperimenti.

Alvise Sommariva  Quadratura numerica  98/ 116

Una demo di Gauss-Jacobi: scaling

L’idea dello scaling consiste nel modificare la formula di quadratura relativa ad un intervallo [a, b] di riferimento così da poterla utilizzare in un intervallo generico [α, β], cambiando esclusivamente nodi e pesi.

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

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

▶ 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]}.

Una demo di Gauss-Jacobi

Nota.

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.

Esempio 1

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

>> l e n g t h ( x_jac )
ans =

11

>> f p r i n t f ( ’ \n \ t [ GAUSS−LEGENDRE ] : \%15.20 f ’ , I_jac ) ;
[ GAUSS−LEGENDRE ] : 0.09523809523809564900
>> 0.09523809523809523300 −0.09523809523809564900
ans =
−4.163336342344337e−016

>>

Alvise Sommariva  Quadratura numerica  101/ 116

Esempio 2

Consideriamo ora l’integrale

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

(37)

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.

Figura: Grafico che illustra l’errore delle formule di quad. dei trap. composta, Cav.-Simpson comp., Gauss-Legendre e Clenshaw-Curtis [10] sulla funzione test (37).

Alvise Sommariva  Quadratura numerica  104/ 116

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. 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

( 1 0 , 1 / 2 , 0 , a , b , @exp ) ;
>> f p r i n t f ( ’ \n \ t [ GAUSS−JACOBI ] : %15.20 f ’ , I_jac ) ;
[ GAUSS−JACOBI ] : 1.77914365469190930000

> >1.77914365469190979259117902999 −1.779143654691909300

ans = 4.440892098500626 e−016
>> l e n g t h ( x_jac )
ans = 10

>> [ I_jac , x_jac , w_jac ] = . . .
  integrazione_gauss_jacobi (10 ,0 ,0 , −1 ,1 , inline ( ’ exp ( x ) . ∗

s q r t (1−x ) ’ ) ) ;
>> f p r i n t f ( ’ \n \ t [ GAUSS−LEGENDRE ] : %15.20 f ’ , I_jac )
[ GAUSS−LEGENDRE ] : 1.77984112101478020000

> >1.77914365469190979259117902999 −1.779841121014780200

ans =−6.9747e−004
>> l e n g t h ( x_jac )
ans = 10

>>

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

>> l e n g t h ( 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).

Esercizio

Usando luna opportuna formula gaussiana, con n = 10, 20, 30 punti, calcolare

▶ ∫_{0}^{π} exp (x) · cos (4x)dx = \frac{exp (π)−1}{17};

▶ ∫_{0}^{1} x^{5/2}dx = \frac{2}{7};

▶ ∫_{−π}^{π} exp cos (x)dx = 7.95492652101284;

▶ ∫_{0}^{π/4} exp cos (x)dx = 1.93973485062365;

▶ ∫_{0}^{1} x^{1/2}dx = \frac{2}{3}. Descrivere come decresce l’errore in scala semilogaritmica. \underline{Non serve vedere} quale sia l’andamento del rapporto E_{n}(f )/E_{n+1}(f ), dove E_{n} = |I (f ) − I_{n}(f )| con I (f ) integrale esatto e I_{n}(f ) valore fornito dalla formula di quadratura.

Esercizio (facoltativo)

Si calcolino per N = 10, 20 con la formula composta dei trapezi, di Cav.-Simpson e un’appropriata formula gaussiana

∫_{−1}^{1} x^{20} dx = 2/21 ≈ 0.095238095238096801  (38)
∫_{−1}^{1} e^{x} dx = e − e^{−1} ≈ 2.3504023872876032  (39)
∫_{−1}^{1} e^{−x2} dx = erf(1)·\sqrt{π} ≈ 1.4936482656248538  (40)
∫_{−1}^{1} 1/(1 + 16x^{2}) dx = 1/2 · atan(4) ≈ 0.66290883183401628  (41)
∫_{−1}^{1} e^{−1/x2} dx ≈ 0.17814771178156086  (42)
∫_{−1}^{1} |x|^{3} dx = 1/2  (43)
∫_{0}^{1} \sqrt{x} dx = 2/3  (44)
∫_{−1}^{1} e^{x} \sqrt{1 − x} dx ≈ 1.7791436546919095  (45)

Esercizio (facoltativo)

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.

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.

Calcolare con una opportuna formula gaussiana basata su 11 nodi: 1.

∫_{−1}^{1} x^{20}\sqrt{1 − x^{2}}dx

2.

∫_{−1}^{1} exp (x^{2}) · (1 − x)^{0.5} · (1 + x)^{1.5}dx

Il primo integrale, può essere calcolato esattamente con Mathematica, Maple (da Matlab, aiutarsi con Google) o quadgk di Matlab, una toolbox di Matlab specifica per l’integrazione numerica, con un errore assoluto di 10^{−15}

>> Q = quadgk ( @ ( x ) ( x . ˆ 2 0 ) .∗(1 −x . ˆ 2 ) . ˆ ( 0 . 5 ) , −1 ,1 , ’
AbsTol ’ ,10ˆ ( −15) )

Q =

0.025160880188796 >>

Alvise Sommariva  Quadratura numerica  112/ 116

In alternativa usando online il Wolfram Integrator, e ricordando che devo integrare in (−1, 1), dal teorema fondamentale del calcolo si vede che si tratta di valutare

\frac{14549535}{1816657920} · (arcsin (1) − arcsin (−1)).

Utilizzando la shell di Matlab

>> c =14549535/1816657920;
>> s=a s i n ( 1 )−a s i n (−1) ;
>> I=s∗c ;
>> format long ;
>> I

I =

0.025160880188796 >>

Esercizio (facoltativo)

Il secondo integrale non è integrabile col Wolfram integrator e nemmeno col calcolo simbolico di Matlab.

>> syms x

>> int ( ( exp(−x . ˆ 2 ) ) .∗((1 −x ) . ˆ ( 0 . 5 ) ) .∗((1+x ) . ˆ ( 1 . 5 ) )
, −1 ,1)

Warning : Explicit integral could not be found .
> In sym . int at 58

ans =
int ( exp(−x ˆ 2) ∗(1−x ) ˆ (1/2) ∗(1+x ) ˆ (3/2) , x = −1 . . 1)

Esercizio (facoltativo)

Utilizzando la shell di Matlab e quadgk con un errore assoluto di 10^{−15}

>> format long ;
>> I=quadgk ( @ ( x ) ( exp ( x . ˆ 2 ) ) .∗((1 −x ) . ˆ ( 0 . 5 ) ) .∗((1+x )
. ˆ ( 1 . 5 ) ) , − 1 , 1 , …
’ AbsTol ’ ,10ˆ ( −15) )

I =

2.086318843895086 >>

Bibliografia

K. Atkinson, Introduction to Numerical Analysis, Wiley, 1989. K. Atkinson e W. Han, Theoretical Numerical Analysis, Springer, 2001. V. Comincioli, Analisi Numerica, metodi modelli applicazioni, Mc Graw-Hill, 1990. S.D. Conte e C. de Boor, Elementary Numerical Analysis, 3rd Edition, Mc Graw-Hill, 1980. P.J. Davis, Interpolation and Approximation, Dover, 1975. W. Gautschi, personal homepage, http://www.cs.purdue.edu/people/wxg. W. Gautschi, , Orthogonal polynomials (in Matlab), JCAM 178 (2005) 215–234 A. Quarteroni e F. Saleri, Introduzione al calcolo scientifico, Springer Verlag, 2006. A. Suli e D. Mayers, An Introduction to Numerical Analysis, Cambridge University Press, 2003. L.N. Trefethen, Is Gauss quadrature better than Clenshaw-Curtis?, SIAM Reviews, (2007). J. Waldvogel, Fast Construction of the Fejer and Clenshaw-Curtis Quadrature Rules,BIT 46 (2006),no. 1,195–202. Wikipedia, Clenshaw-Curtis quadrature, http://en.wikipedia.org/wiki/Clenshaw-Curtis quadrature. Wikipedia, Quadratura gaussiana, http://it.wikipedia.org/wiki/Quadratura di Gauss. Wikipedia, Regola del trapezio,http://it.wikipedia.org/wiki/Regola del trapezio. Wikipedia, Regola di Cavalieri-Simpson, http://it.wikipedia.org/wiki/Regola di Cavalieri-Simpson.

Alvise Sommariva  Quadratura numerica  116/ 116