Quadratura numerica

19 febbraio 2009

1 Introduzione

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

( ) : _{}( ) _{}^{}( ) ( )

dove è una funzione peso in ( ).

La nostra intenzione è di approssimare ( ) come

( ) _{}( ) : _{1}^{}_{}( _{})  (1)

I termini _{}e _{}[ ] sono detti rispettivamente pesi e nodi.

2 Formule interpolatorie

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

; ; ; _{}^{}( ) ( ) ; ; ; _{}^{}( ) ( ) _{}_{1}

Se

1 ( ) _{1}^{}( _{}) _{}( )

è il polinomio che interpola le coppie ( _{}( _{})) con 1 , dove al solito _{} indica l’i-simo polinomio di Lagrange allora

_{}^{}( ) ( )  _{}^{}_{1} ( ) ( )
_{}^{}_{1}^{}( _{}) _{}( ) ( )

_{1}^{}_{}^{}_{}( ) ( )  ( _{})  (2)

per cui, confrontando con la formula (1) abbiamo

_{}^{}_{}( ) ( ) 1

Se _{1} è un polinomio di grado 1 ovviamente corrisponde col polinomio interpolante _{1} nei nodi a due a due distinti _{1} _{}e quindi la formula risulta esatta per polinomi di grado inferiore o uguale a 1, cioè

( _{1} ) _{1}^{}_{}( _{}) _{}_{}^{}_{}( ) ( ) 1

Si dimostra che se una formula avente quali nodi un insieme _{1} _{}di punti a due a due distinti è esatta per ogni polinomio di grado 1 allora è interpolatoria.

Infatti se è esatta per ogni polinomio di grado 1 allora lo è in particolare per i polinomi di Lagrange _{}_{1} , il che implica che _{}_{}^{}_{}( ) ( ) e quindi i pesi sono proprio quelli della formula interpolatoria corrispondente nei nodi _{1} _{}.

3 Formule di Newton-Cotes

Si supponga [ ] un intervallo chiuso e limitato di � e si ponga 1. Per semplicità di notazione, in questo caso porremo : _{}_{1}. Il primo esempio di formule interpolatorie che consideriamo sono le regole di tipo Newton-Cotes chiuse (cf. [13, p.336]) che si ottengono integrando l’interpolante di in nodi equispaziati

\frac{( 1)( )}{1} 1

Alcune classiche regole sono:

1. regola del trapezio

( ) _{1}( ) : _{1}( ) : \frac{( )( ( ) ( ))}{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 [23] è

1( ) : ( ) _{1}( ) \frac{^{3}}{12} ^{(2)}( )

per qualche ( ); Figura 1: Regola del trapezio e di Cavalieri-Simpson per il calcolo di _{0 5}^{2}sin( ) (rispettivamente area in magenta e in azzurro).

2. regola di Cavalieri-Simpson

( ) _{3}( ) : _{3}( ) : \frac{}{6} ( ) 4 (\frac{}{2} ) ( )

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 [24] è

3( ) : ( ) _{3}( ) \frac{^{5}}{90} ^{(4)}( ) \frac{}{2}

per qualche ( );

Vediamo calcolando i pesi, che in effetti le due formule sono interpolatorie. Partiamo dalla regola del trapezio. Posti _{1} , _{2} abbiamo che

1( ) \frac{}{}_{2}( ) \frac{}{}

e quindi visto che 1 abbiamo

1  _{}^{}_{1}( )  _{}^{}\frac{}{}\frac{1}{}_{}^{}( )
\frac{1}{}\frac{( )^{2}}{2} _{}^{}\frac{1}{}\frac{( )^{2}}{2} _{}^{}\frac{1}{}\frac{( )^{2}}{2} \frac{}{2}  (3)

e

2  _{}^{}_{2}( )  _{}^{}\frac{}{}\frac{1}{}_{}^{}( )
\frac{1}{}\frac{( )^{2}}{2} _{}^{}\frac{1}{}\frac{( )^{2}}{2} _{}^{}\frac{1}{}\frac{( )^{2}}{2} \frac{}{2}  (4)

Figura 2: Cavalieri e Simpson.

Per quanto riguarda la formula di Cavalieri 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], [13, 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 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) [ ] in subintervalli _{}[ _{}_{1}] tali che _{}con ( )/ . Dalle proprietà dell’integrale

_{}^{}( ) _{0}^{1} _{_{}}^{_{1}} ( ) _{0}^{1} ( _{}_{1} )  (5)

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

1. formula composta dei trapezi

_{1}^{( )}: \frac{( _{0})}{2} ( _{1})  ( _{1} ) \frac{( _{})}{2}  (6)

il cui errore è

_{1}^{( )}( ) : ( ) _{1}^{( )}( ) \frac{( )}{12} ^{2} ^{(2)}( ) \frac{( )}{}

per qualche ( );

2. formula composta di Cavalieri-Simpson fissati il numero di subintervalli e i punti _{}/2 dove

\frac{}{}

sia

( ) _{3}^{( )}( ) : \frac{}{6} ( _{0}) 2_{1}^{1} ( _{2}) 4_{0}^{1} ( _{2 1}) ( _{2} ) ; (7)

il cui errore è

_{3}^{( )}( ) : ( ) _{3}^{( )}( ) \frac{( )}{180} \frac{}{2} ^{4} ^{(4)}( )

per qualche ( ).

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

( ) _{1}^{}( ) : _{0}^{}_{}( _{})  (8)

come pure per (6)

_{1}^{( )}: \frac{( _{0})}{2} ( _{1})  ( _{1} ) \frac{( _{})}{2}  (9)

deduciamo che _{0} _{}/2 mentre _{1} _{1} , cosa che giustifica le ultime linee della function trapezi_composta.

Si potrebbe usare il comando Matlab trapz nella sua implementazione Figura 3: Formula dei trapezi composta per il calcolo di _{0 5}^{2}sin( ) (area in magenta).

>> 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( ) cos(1) ( cos(0)) cos(1) 1 0 45969769413186

sia utilizzando la funzione trapezi_composta che trapz

>> 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 (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) e (7) ne spiega il significato.

Una volta noti il vettore (colonna) dei nodi e dei pesi di integrazione, se la funzione è 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_{})_{0} �^{1} , fx ( ( _{}))_{0} �^{1} sono due vettori colonna allora il prodotto scalare

w fx: _{0}^{}w_{}fx_{}: _{0}^{}w_{}( _{})

si scrive in Matlab/Octave come w'*fx. Osserviamo che dimensionalmente il prodotto di un vettore 1 ( 1) con un vettore ( 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} ^{20} (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. [21], [22]) o di tipo Clenshaw-Curtis (cf. [19], [15], [17]) a parità di valutazioni della funzione 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 [20].

4.2 Esercizio

Si calcolino con la formula composta dei trapezi e di Cavalieri-Simpson, con 5 7 9 11 13 15 i seguenti integrali

_{1}^{1} ^{20}  2/21 0 095238095238096801  (10)
_{1}^{1} ^{}^{1} 2 3504023872876032  (11)
_{1}^{1} ^{^{2}}  erf(1) \overline{(} ) 1 4936482656248538  (12)
_{1}^{1} 1/(1 16 ^{2})  1/2 atan(4) 0 66290883183401628  (13)
_{1}^{1} ^{1/ ^{2}}  2 ^{1} 2 ^{1} erf(1) ^{1/2} ^{1} 2 ^{1/2} (14)
0 17814771178156086  (15)
_{1}^{1}  ^{3}  1/2  (16)
_{0}^{1} \overline{}2/3  (17)
_{1}^{1} ^{}\overline{(}1 )  2^{1/2} ^{1} 1/2  ^{1/2} erf(2^{1/2})  (18)
1 7791436546919095  (19)

Quali delle due formule ha errori relativi inferiori? Quali funzioni risultano più difficili da integrare numericamente?

5 Formule gaussiane

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

Sia : ( ) � (non necessariamente limitato) è una funzione peso, cioè tale che (cf. [1, p.206, p.270])

1. è nonnegativa in ( );

2. è integrabile in [ ];

3. esista e sia finito

_{}^{}^{}( )

per ogni �; 4. se

_{}^{}( ) ( )

per una qualche funzione nonnegativa allora 0 in ( ).

Tra gli esempi più noti ricordiamo

1. Legendre: ( ) 1 in [ ] limitato;

2. Jacobi: ( ) (1 )^{}(1 )^{}in ( 1 1) per , 1;

3. Chebyshev: ( ) \frac{1}{\overline{1 ^{2}}} in ( 1 1);

4. Laguerre: ( ) exp( ) in [0 );

5. Hermite: ( ) exp( ^{2}) in ( );

Si supponga ora di dover calcolare per qualche funzione : ( ) �

( ) : _{}^{}( ) ( )

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

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

_{}^{}_{}( ) ( )

per ?

Figura 4: Hermite, Jacobi e Laguerre.

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

Teorema 5.1 Per ogni 1 esistono e sono unici dei nodi _{1}  _{}e pesi _{1}  _{}

per cui il grado di precisione sia 2 1. I nodi sono gli zeri del polinomio ortogonale di grado ,

( ) _{}( _{1}) ( _{})

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

Sia _{2 1}  _{2 1} e _{1} _{1}  _{1} tali che
2 1 _{1} _{}_{1}

Allora poichè

• _{}^{}_{1} ( ) _{}( ) ( ) 0, poichè _{}è il polinomio ortogonale rispetto di grado ;

• la formula è interpolatoria, per cui esatta per ogni polinomio di grado 1 in quanto basata su punti a due a due distinti;

• se _{}è uno zero di _{}allora

2 1 ( _{}) _{1} ( _{}) _{}( _{}) _{1} ( _{}) _{1} ( _{})

abbiamo

_{}^{}_{2 1} ( ) ( )  _{}^{}_{1} ( ) _{}( ) ( )  _{}^{}_{1} ( ) ( )
0 _{}^{}_{1} ( ) ( )
_{1}^{}_{}_{1} ( _{})
_{1}^{}_{}_{2 1} ( _{})  (20)

per cui la formula di Gauss ha grado di precisione 2 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 _{1} _{}per cui

0 _{}^{}_{}^{2}( ) ( ) _{1}^{}_{}_{}^{2}( _{}) _{}

Se esistesse un’altra formula interpolatoria con grado di precisione 2 1 e avesse nodi _{}_{ ˜ 1} e pesi _{}_{ ˜ 1} per prima cosa i pesi sarebbero positivi poichè il grado di precisione è 2 1 e quindi sarebbe esatta per il -simo polinomio di Lagrange _{} ˜ da cui

0 _{}^{}˜_{}^{2}( ) ( ) _{1}^{}_{}˜˜_{}^{2}( _{} ˜ ) _{} ˜

D’altra parte se _{} ˜ è il -simo polinomio di Lagrange, poichè _{}è il polinomio ortogonale di grado rispetto al peso , e _{} ˜ 0 abbiamo che da

0 ( _{}_{}˜ )_{}_{1}^{}_{}˜_{}˜ ( _{} ˜ ) _{}( _{} ˜ ) _{}_{} ˜ ( _{} ˜ )

necessariamente _{}_{}˜ e visto che quindi _{}_{}˜ ricaviamo anche

_{}^{}_{}^{2}( ) ( )  _{}^{}˜_{}^{2}( ) ( ) _{}  ˜

per cui la formula gaussiana cercata è unica.

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

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

1 ( ) _{}( _{}) _{}( ) _{}_{1} ( )

ove, detto _{}il coefficiente di grado in _{}, si ha _{}_{1} / _{},

\frac{( _{}_{})}{_{}^{2}} _{}\frac{( _{}_{1} )}{_{}^{2}}

con ( ) _{}^{}( ) ( ) . Per ulteriori dettagli sui polinomi ortogonali si confronti [1, p.270-283], [3, p.978-985], [13, 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 , uguale al massimo al numero di righe della matrice di coefficienti di ricorrenza , fornisce nodi e pesi immagazzinati in una matrice che ha quale prima colonna e quale seconda colonna . Osserviamo che dall’help si evince che i polinomi ortogonali sono monici e quindi _{}1. Inoltre ( )_{} ( _{}_{}).

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.

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

.

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

5.2 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 [15] 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 1, 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. Una lista di possibili tests è la seguente:

_{1}^{1} ^{20}  2/21 0 095238095238096801  (21)
_{1}^{1} ^{}^{1} 2 3504023872876032  (22)
_{1}^{1} ^{^{2}}  erf(1) \overline{(} ) 1 4936482656248538  (23)
_{1}^{1} 1/(1 16 ^{2})  1/2 atan(4) 0 66290883183401628  (24)
_{1}^{1} ^{1/ ^{2}}  2 ^{1} 2 ^{1} erf(1) ^{1/2} ^{1} 2 ^{1/2} (25)
0 17814771178156086  (26)
_{1}^{1}  ^{3}  1/2  (27)
_{1}^{1} \overline{}2/3  (28)
_{1}^{1} ^{}\overline{(}1 )  2^{1/2} ^{1} 1/2  ^{1/2} erf(2^{1/2})  (29)
1 7791436546919095  (30)

Alcune note prima di cominciare.

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

[ ] (( )/2) (( )/2) _{[ 1 1]};

2. se _{[ 1 1]}sono i pesi di Gauss-Jacobi in [ 1 1] allora i pesi di Gauss-Jacobi [ ]in [ ] sono

[ ] (( )/2) _{[ 1 1]}

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

_{}_{[ ]} _{}^{}( )

mentre nel caso di 1, 1 derivante dalla funzione peso di Jacobi abbiamo

_{}_{[ 1 1]} _{1}^{1} ( ) 2;

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

3. Nonostante l’introduzione riguardante nodi e pesi in [ ] 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 , , fornire il relativo risultato esatto in exact_results.m e testare la formula gaussiana. Per avere il risultato con alta precisione si usi la funzione quadl.m o quad8.m di Matlab oppure

integrazione_gauss_jacobi.m con 500;

4. l’intervallo tipico di Gauss-Legendre è [ ] (chiuso), mentre per Gauss-Jacobi l’intervallo è tipicamente ( ) (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 , sono un insieme trascurabile.

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

_{1}^{1} ^{20} 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
>>

5.3 Facoltativo. Una implementazione di Clenshaw-Curtis

La quadrature di tipo Clenshaw-Curtis si basa su un’espansione dell’integranda in termini di polinomi ortogonali di Chebyshev (cf. [19], [15], [17]). Usiamo una implementazione di Waldvogel [17]

function [x,w]=clenshaw_curtis(n,a,b)
N=[1:2:n-1]'; l=length(N); m=n-l; K=[0:m-1]';

g0=-ones(n,1); g0(1+l)=g0(1+l)+n; g0(1+m)=g0(1+m)+n;
g=g0/(n^2-1+mod(n,2));
end_N=length(N);
v0=[2./N./(N-2); 1/N(end_N); zeros(m,1)];
end_v0=length(v0);
v2=-v0(1:end_v0-1)-v0(end_v0:-1:2);
wcc=ifft(v2+g); weights=[wcc;wcc(1,1)];
k=0:n; nodes=(cos(k*pi/n))';
x=(a+b)/2+((b-a)/2)*nodes;
w=weights*((b-a)/2);

Noti i nodi e i pesi di Clenshaw-Curtis, definiamo la function integrazione_clenshaw_curtis

function [I_cc,x_cc,w_cc]=integrazione_clenshaw_curtis(N_cc,a,b,f)
% INTEGRAZIONE DI CLENSHAW-CURTIS.
% INPUT:
% N_cc: GRADO DI CLENSHAW-CURTIS, CIOE' NUMERO DI NODI MENO UNO.
% a, b: ESTREMI DI INTEGRAZIONE.
% f   : FUNZIONE DA INTEGRARE.
% OUTPUT:
% I_cc: APPROSSIMAZIONE DELL'INTEGRALE CON FORMULA DI CLENSHAW-CURTIS.
% x_cc: NODI DI CLENSHAW-CURTIS.
% w_cc: PESI DI CLENSHAW-CURTIS.
% ROUTINES ESTERNE: clenshaw_curtis.
[x_cc,w_cc]=clenshaw_curtis(N_cc,a,b);  % NODI E PESI DI CL.-CURTIS.
fx_cc=feval(f,x_cc);                    % VALUTAZIONE FUNZIONE.
I_cc=w_cc'*fx_cc;                       % VALORE INTEGRALE.

Applichiamola ora al calcolo di

_{1}^{1} ^{20} 2/21 0 09523809523809523300

>> N_cc=10; a=-1; b=1;
>> [I_cc,x_cc,w_cc]=integrazione_clenshaw_curtis(N_cc,a,b,@f);
>> fprintf('\n \t [CLENSHAW-CURTIS]: %15.20f',I_cc);
          [CLENSHAW-CURTIS]: 0.09490517620400430700
>> 0.09490517620400430700-2/21
ans =
   -3.329190340909255e-004
>> length(x_cc)

ans =
    11
>>

La formula ha 11 nodi e pesi, ed ha un’errore assoluto di circa 3 32 10^{4} . Osserviamo il passaggio della funzione attraverso @f.

Dalla figura con i 6 grafici di errore, si vedono le performances delle formule sopracitate quando applicate a calcolare

_{1}^{1} ^{20}  (a sinistra, prima riga)  (31)
_{1}^{1} ^{}(a destra, prima riga)  (32)
_{1}^{1} ^{^{2}}  (a sinistra, seconda riga)  (33)
_{1}^{1} 1/(1 16 ^{2})  (a destra, seconda riga)  (34)
_{1}^{1} ^{1/ ^{2}}  (a sinistra, terza riga)  (35)
_{1}^{1}  ^{3}  (a destra, terza riga)  (36)

Da osservare che fissato , le formule composte hanno nodi mentre le gaussiane e Clenshaw-Curtis ne hanno 1.

5.4 Un’integrale con funzione peso

Consideriamo ora l’integrale

_{1}^{1} exp( ) \overline{1}  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)

Figura 6: Grafico che illustra l’errore delle formule di quadratura dei trapezi composta, Cavalieri-Simpson composta, Gauss-Legendre e Clenshaw-Curtis su 6 funzioni test.

Figura 7: Grafico che illustra l’errore delle formule di quadratura dei trapezi composta, Cavalieri-Simpson composta, Gauss-Legendre e Clenshaw-Curtis sulla funzione test (37).

calcola simbolicamente l’integrale

_{1}^{1} exp( ) \overline{1}

Dal grafico in figura 37, si vede come pure la formula di Gauss-Legendre non abbia una grande performance.

E’ immediato osservare che ( ) \overline{1} è un peso di Gauss-Jacobi

( ) (1 )^{}(1 )^{}

per 1/2 e 0.

Infatti se ( ) \overline{1} , allora ( ) exp( ) ( ) il che corrisponde a usare Gauss-Jacobi con ( ) exp( ). 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 1/2, 0 e di 2 52 10^{3} per Gauss-Legendre (cioè Gauss-Jacobi con 0, 0).

5.5 Esperimenti in Octave

Per via della compatibilità Matlab/Octave, ricapitoliamo quanto visto in questa lezione in Octave (versione 2.1.73). Osserviamo che per Octave, octave:8>, octave:9>, corrisponde in Matlab a >>.

octave:8> format long;
octave:9> [x,w]=trapezi_composta(10,0,1);
octave:10> fx=sin(x);
octave:11> I_trapezi_composta=w'*fx
I_trapezi_composta = 0.459314548857976
octave:12> h=(1-0)/10;
octave:13> I_trapz=h*trapz(fx)
I_trapz = 0.459314548857976
octave:14> fx=f(x);
octave:15> I=w'*fx
I = 0.0633941318588326
octave:16> demo_composte
x_simp =
   -1.00000000000000e+00
   -8.00000000000000e-01
   -6.00000000000000e-01
   -4.00000000000000e-01
   -2.00000000000000e-01
    5.55111512312578e-17
    2.00000000000000e-01
    4.00000000000000e-01
    6.00000000000000e-01
    8.00000000000000e-01
    1.00000000000000e+00

w_simp =
0.0666666666666667
0.2666666666666667
0.1333333333333333
0.2666666666666667
0.1333333333333333
0.2666666666666667
0.1333333333333333
0.2666666666666667
0.1333333333333333
0.2666666666666667
0.0666666666666667

       [TRAPEZI COMPOSTA] [PTS]:   11
       [TRAPEZI COMPOSTA] [RIS]: 0.20462631505024
       [SIMPSON COMPOSTA] [PTS]:   11
       [SIMPSON COMPOSTA] [RIS]: 0.13949200364447
octave:17> % ESEMPIO GAUSS JACOBI.
octave:18> format long;
octave:19> [I_jac_x,x_jac,w_jac]=integrazione_gauss_jacobi(N,ajac,bjac,a,b,@f)
I_jac_x = 0.0952380952380953
x_jac =
 -9.78228658146057e-01
 -8.87062599768095e-01
 -7.30152005574049e-01
 -5.19096129206812e-01
 -2.69543155952345e-01
  1.25162426693849e-17
  2.69543155952345e-01
  5.19096129206812e-01
  7.30152005574049e-01
  8.87062599768095e-01
  9.78228658146057e-01
w_jac =
0.0556685671161738
0.1255803694649046
0.1862902109277342
0.2331937645919905
0.2628045445102464
0.2729250867779003
0.2628045445102470
0.2331937645919902
0.1862902109277343
0.1255803694649046

  0.0556685671161738
octave:20> % ESEMPIO CLENSHAW-CURTIS.
octave:21> N_cc=10; a=-1; b=1;
octave:22> [I_cc,x_cc,w_cc]=integrazione_clenshaw_curtis(N_cc,a,b,@f);
octave:23> fprintf('\n \t [CLENSHAW-CURTIS]: %15.20f',I_cc)
         [CLENSHAW-CURTIS]: 0.09490517620400433507
octave:24>  0.09490517620400433507-2/21
ans =  -3.32919034090898e-04
octave:25> length(x_cc)
ans = 11
octave:26> % ESEMPIO GAUSS-JACOBI (-1/2,0)
octave:27> % PER exp(x).*sqrt(1-x.^2)
octave:28> a=-1; b=1;
octave:29> [I_jac,x_jac,w_jac]=integrazione_gauss_jacobi(10,1/2,0,a,b,@exp);
octave:30> fprintf('\n\t [GAUSS-JACOBI]: %15.20f',I_jac)
         [GAUSS-JACOBI]: 1.77914365469190971503
octave:31> 1.7791436546919097925911790299941-1.77914365469190971503
ans = 0
octave:32> format long e
octave:33> 1.7791436546919097925911790299941-1.77914365469190971503
ans =  0.00000000000000e+00
octave:34> format short
octave:35> length(x_jac)
ans = 10
octave:36> % ESEMPIO GAUSS-LEGENDRE (CIOE' GAUSS-JACOBI (0,0)).
octave:37> % PER exp(x).*sqrt(1-x.^2)
octave:38> [I_jac,x_jac,w_jac]=integrazione_gauss_jacobi(10,0,0,...
            a,b,inline('exp(x).*sqrt(1-x)'));
octave:39> [I_jac,x_jac, w_jac]=integrazione_gauss_jacobi(10,0,0,...
        -1,1,inline('exp(x).*sqrt(1-x)'));
octave:40> fprintf('\n \t [GAUSS-LEGENDRE: %15.20f]',I_jac)
         [GAUSS-LEGENDRE: 1.77984112101478153534
octave:41> 1.7791436546919097925911790299941-I_jac
ans =  -6.9747e-04
octave:42> length(x_jac)
ans = 10

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

( ) _{}( ) _{0}^{}_{}( _{})

una regola di Newton-Cotes.

1. se n è pari e ^{( 2)} ([ ]) allora

( ) _{}( ) _{}^{3} ^{( 2)} ( ) ( )

con

\frac{1}{( 2)!} _{0}^{}^{2}( 1) ( ) ;

2. se n è dispari e ^{( 1)} ([ ]) allora

( ) _{}( ) _{}^{2} ^{( 1)} ( ) ( )

con

\frac{1}{( 1)!} _{0}^{}( 1) ( ) ;

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 ^{(2 )}( ) con ( ) compatto e supponiamo

( ) _{}^{}( ) ( ) _{}( ) _{1}^{}_{}( _{})

sia una formula gaussiana rispetto alla funzione peso . Allora

( ) : _{}( ) _{}( ) \frac{_{}}{_{}^{2}(2 )!} ^{(2 )}( ) ( )

dove _{}è il coefficiente di grado massimo del polinomio ortogonale _{}di grado , _{}^{}_{}^{2}( ) ( ) .

In particolare, se 1, [ ] [ 1 1] allora

( ) \frac{2^{2 1} ( !)^{4}}{(2 1)[(2 )!]^{3}} ^{(2 )}( ) ( 1 1)

7 Stabilità di una formula di quadratura

Sia [ ] un intervallo compatto e una funzione peso in ( ). Inoltre supponiamo

( ) _{1}^{}_{}_{}con _{}( _{})  (38)

e che invece di _{}_{}si disponga di una loro approssimazione _{}_{} ˜ . Di conseguenza al posto di _{}( ) si calcola

˜ ( ) _{1}^{}_{}_{}˜

ed è

( ) _{} ˜ ( )  _{1}^{}_{}( _{}_{}˜ ) _{1}^{}_{}_{}_{}˜

_{1}^{}_{}max_{}_{}_{}˜  (39)

Quindi la quantità

_{1}^{}_{}

è un indice di stabilità della formula di quadratura.

Si dimostra che in effetti l’operatore lineare _{}: ([ ]) � è continuo in quanto

( )   ;  ;  ;  ; _{1}^{}_{}_{} ;  ;  ;  ; _{1}^{}_{}_{}

_{1}^{}_{}max_{}_{}_{1}^{}_{}_{}(40)

e che la norma dell’operatore di quadratura

max_{([ ]) 0}\frac{_{}( )}{_{}}

coincide con _{1}^{}_{}.

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 ( ) un intervallo limitato e si supponga : ( ) � sia una funzione peso e

( ) _{1}^{}_{}_{}con _{}( _{})

una formula di quadratura avente grado di precisione . Allora

( ) _{}( ) ( _{1} _{}_{}) min_{_{}_{}}_{}_{}

Dimostrazione. Se _{}_{}è un polinomio arbitrario di grado ed ( ) _{}( ) avendo la formula di quadratura grado di precisione , essendo

( ) ; ; ; _{}^{}( ) ( ) ; ; ; _{}_{1}

poichè per definizione ( ) _{}_{}_{}

( ) _{}( )  ( ) _{}( ) _{}( ) _{}( )
( ) _{}( )  _{}( ) _{}( )
( ) ( )  _{}( )
( )  _{}( )
1  _{}_{}_{}_{}
( _{1}  _{}_{})  _{}(41)

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

( ) _{}( ) 2 _{1} min_{_{}_{}}_{}_{}

Per stime delle quantità min_{_{}_{}}_{}_{}

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

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

( ) _{}( ) 4 min_{_{}_{}}_{}_{}

Nota. Supponiamo che alcuni pesi siano negativi. Di conseguenza

_{}^{}( )  _{}_{}^{}_{}_{}^{}(42)
_{}_{}^{}_{}_{}^{}_{}_{}^{}_{}_{}^{}(43)

da cui

_{}_{}^{}_{}_{}^{}_{}^{}( ) _{}_{}^{}_{}_{}^{}_{}^{}( ) 2_{}_{}^{}

mentre se i pesi fossero tutti positivi avremmo

_{}_{}^{}_{}_{}^{}_{}_{}^{}_{}^{}( )

per cui 2 _{}_{}^{}0 peggiora la stabilità della formula di quadratura.

Figura 8: G. Polya (1887-1985) e V.A. Steklov (1864-1926).

8 Teoremi di convergenza

Consideriamo una formula di quadratura del tipo

_{}^{}( ) ( )  _{}( ) : _{0}^{}_{}( _{})  (44)

con al solito una funzione peso definita nell’intervallo ( ). Se tale intervallo è limitato ed continua in [ ], allora ^{1}( ). Sia

( ) : _{}^{}( ) ( ) _{0}^{}_{}( _{})

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

Teorema 8.1 Siano _{}_{}dei punti di un intervallo fissato e chiuso [ ]. Condizione necessaria e sufficiente affinchè per ogni funzione continua ( ) si abbia

lim_{}_{}( ) 0

è che sia

1. esiste � tale che per ogni si abbia

_{0}^{}_{};

2. per ogni � si abbia lim_{}_{}( ^{}) 0

Dimostrazione. Supponiamo che

1. esiste � tale che per ogni si abbia

_{0}^{}_{};

2. per ogni � si abbia lim_{}_{}( ^{}) 0

Per un teorema di densità dovuto a Weierstrass, per ogni 0 esiste un polinomio tale che _{}_{}e quindi

( _{})   ;  ;  ; _{}^{}( ( ) _{}( )) ( )  ;  ;  ;  ;  ;  ;  ; _{0}^{}_{}( ( _{}) _{}( _{})) ;  ;  ;  ;
1 _{0}^{}_{}( _{}) _{}( _{})

1 _{0}^{}_{}_{}_{}

1 _{0}^{}_{}(45)

Di conseguenza

( )  _{}( ) _{}( )  _{}( )
( _{1} _{0}^{}_{})  _{}( )
( _{1} )  _{}( )

(46)

e siccome _{}( ) 0 al crescere di poichè se ( ) _{0}^{}_{}^{}allora

( ) _{0}^{}_{}_{}( ^{}) 0

abbiamo

lim_{}_{}( ) ( _{1} ) 0

da cui per l’arbitrarietà di abbiamo lim_{}_{}( ) 0.

Viceversa supponiamo lim_{}_{}( ) 0 per ogni ([ ]). Si vede che

( )  _{}^{}( ) ( ) _{0}^{}_{}( _{})
_{}^{}( ) ( )  _{0}^{}_{}( _{})
1  _{}_{0}^{}_{}( _{})
1  _{}(_{0}^{}_{}) _{}

1 _{0}^{}_{}_{}(47)

In realtà si prova che

max_{([ ]) 0}\frac{_{}( )}{}_{1} _{0}^{}_{}

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

sup_{}_{}

Nel nostro caso

• ( ([ ]) _{}), � sono spazi di Banach,

• posto _{}_{}, se ([ ]) abbiamo che _{}( ) è limitata in quanto per ipotesi è infinitesima,

e di conseguenza

sup_{}_{0}^{}_{}sup_{}_{1} _{0}^{}_{}sup_{}_{}_{}

per cui esiste finito ed indipendente da tale che _{0}^{}_{}.

Il secondo punto da dimostrare è ovvio in quanto per ogni , si ha ^{}([ ]) . Notiamo che

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

2. Si consideri una suddivisione _{}_{}_{0} dell’intervallo ( ) con _{} 1, _{0} , _{}; si supponga che _{}_{1} _{}e che sia

max_{0 1} _{}

Se l’ -simo intervallo (con 1 ) della suddivisione contiene i punti

1 _{( 1)} _{( 1) 1} _{}_{}

sia _{_{}}( ) la funzione polinomiale a tratti di grado ottenuta interpolando nell’ -simo subintervallo con un polinomio di grado le coppie

( _{( 1)} ( _{( 1)}))  ( _{}( _{}))  1

Come è noto, se ^{( 1)} ( ), per qualche ( )

\frac{^{( 1)} ^{( 1)} ( )}{( 1)!} \frac{^{( 1)} ^{( 1)} _{}}{( 1)!}

Se consideriamo una formula composta basata sull’applicare una formula di Newton-Cotes su 1 pesi positivi in ognuno degli 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

_{0}^{}_{}_{0}^{}_{}_{0}^{}_{}1 _{}^{}1

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

_{}^{}( )  _{}^{}( ( ) _{_{}})  _{}^{}_{_{}}
_{}^{}( ( ) _{_{}}) _{0}^{}_{}( _{})  (48)

da cui

( )   ;  ;  ;  ; _{}^{}( ) _{0}^{}_{}( _{}) ;  ;  ;  ;

 ;  ;  ; _{}^{}( ( ) _{_{}})  ;  ;  ;
 ;  ;  ; _{}^{}1  ;  ;  ; _{_{}}_{}
( )  _{_{}}_{}
( ) \frac{^{( 1)} ^{( 1)} _{}}{( 1)!}

(49)

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

_{}^{}( ) ( ) _{}_{}( _{})  (50)

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

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

_{0}^{}_{}_{0}^{}_{}_{0}^{}_{}1 _{}^{}1 ( ) _{1}

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

lim_{}_{}( ) 0

qualsiasi sia la funzione continua ( ).

9 Equazioni di Fredholm di seconda specie

Si consideri l’equazione integrale di Fredholm di seconda specie

( ) _{}^{}( ) ( )  ( )  [ ]  (51)

dove ([ ] [ ]) e ([ ]), �.

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

max_{}( )

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

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

_{}^{}( ) _{1}^{}_{}( _{})  ([ ])  (52)

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

( ) _{1}^{}_{}( _{}) _{}( _{}) ( )  [ ]  (53)

dove _{}( _{}) _{1 _{}}verifica

( _{}) _{1}^{}_{}( _{}_{}) _{}( _{}) ( )  [ ]  (54)

Posti _{}^{( )} _{}( _{}_{}), u_{}_{}( _{}) _{}, f_{}( _{}) _{}si ha così che

u_{}^{( )}u_{}f_{}

per cui se ^{( )} ^{( )} è non singolare, allora

( )u_{}f_{}

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

Una volta ottenuto _{}( _{}), tramite (53), la funzione soluzione _{}può essere valutata in qualsiasi punto [ ]. Per le proprietà di convergenza di questo metodo si veda [2, p.376]

9.1 Esercizio

Si consideri il problema

2 ( ) _{0}^{1}exp( ) ( ) (2exp( ) \frac{1}{1} (exp( 1) 1))  [ ] (55)

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

max_{0 100} _{}( _{}) exp( _{})

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

Alcuni link interessanti sono i seguenti:

http://www.cs.purdue.edu/people/wxg
http://en.wikipedia.org/wiki/Clenshaw-Curtis_quadrature
http://en.wikipedia.org/wiki/Fast_Fourier_transform
http://en.wikipedia.org/wiki/Gaussian_quadrature
http://it.wikipedia.org/wiki/Quadratura_di_Gauss
http://it.wikipedia.org/wiki/Regola_del_trapezio
http://it.wikipedia.org/wiki/Regola_di_Cavalieri-Simpson

Per le biografie di alcuni matematici menzionati in questa lezioni si considerino

http://it.wikipedia.org/wiki/Bonaventura_Cavalieri
http://it.wikipedia.org/wiki/Roger_Cotes
http://www-gap.dcs.st-and.ac.uk/~history/Biographies/Cavalieri.html
http://www-history.mcs.st-andrews.ac.uk/Biographies/Cotes.html
http://www-history.mcs.st-andrews.ac.uk/Biographies/Hermite.html
http://www-history.mcs.st-andrews.ac.uk/Biographies/Jacobi.html
http://www-history.mcs.st-andrews.ac.uk/Biographies/Laguerre.html
http://www-history.mcs.st-andrews.ac.uk/Biographies/Simpson.html

Quale nota storica osserviamo che da [11]

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 _{1} _{}( _{})/ ^{}( _{}) 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)

References

[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] Mac Tutor, Cavalieri, http://www-gap.dcs.st-and.ac.uk/ history/Biographies/Cavalieri.html .

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

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

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

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

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

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

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

[15] 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 .

[16] F.G. Tricomi, Integral Euations, Dover, 1985.

[17] 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 .

[18] Wikipedia, Cavalieri, http://it.wikipedia.org/wiki/Bonaventura_Cavalieri .

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

[20] Wikipedia, Fast Fourier Transform, http://en.wikipedia.org/wiki/Fast_Fourier_transform .

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

[22] Wikipedia, Quadratura gaussiana, http://it.wikipedia.org/wiki/Quadratura_di_Gauss .

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

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