Quadratura numerica

13 febbraio 2008

1 Introduzione

Un classico problema dell’analisi numerica è quello di calcolare l’integrale definito di una funzione nell’intervallo [ ] cioè

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

dove supporremo per il momento che 1. [ ] [ ] � sia un intervallo chiuso e limitato;

2. la funzione : [ ] � sia sufficientemente regolare.

Successivamente indeboliremo queste condizioni. La nostra intenzione è di approssimare ( ) come

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

I termini _{}e _{}[ ] sono detti rispettivamente pesi e nodi. Notiamo che consideriamo [ ] [ ] �: questo perchè in alcune formule si possono considerare nodi esterni al dominio di integrazione.

1.1 Formule di Newton-Cotes

Le prime regole di questo tipo sono dette di Newton-Cotes chiuse (cf. [11, p.336]) e 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}

Figure 1: Cavalieri e Simpson.

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. [2, p.132]) che l’errore della regola del trapezio [20] è

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

per qualche ( );

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

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

per qualche ( );

Per ulteriori dettagli si confronti [1, p.252-258], [11, 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. [2, p.196]).

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

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

1. formula composta dei trapezi

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

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

il cui errore è

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

per qualche ( ).

1.3 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}^{}_{}( _{})  (5)

come pure per (3)

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

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

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

Quadratura numerica

Figure 2: Formula dei trapezi e dei trapezi composta per il calcolo di _{0 5}^{2}sin( ) (area in magenta).

   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 (5) e (4) 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. [18], [19]) o di tipo Clenshaw-Curtis (cf. [16], [13], [14]) 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 [17].

1.4 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 è al massimo uguale al numero di nodi . Ci si domanda se sia possibile scegliere nodi _{1} _{}e pesi _{1} _{}per cui le relative formule di quadratura abbiano grado di precisione (come ad esempio per la regola di Cavalieri-Simpson).

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

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 in (1), visto che l’integranda non é necessariamente continua in [ ] (si consideri ad esempio il peso di Chebyshev) 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 abbiano grado di precisione , cioè calcolino esattamente

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

per ?

Figure 3: Hermite, Jacobi e Laguerre.

La risposta è affermativa, come si può vedere in [1, p.272]. Per ogni 1 esistono e sono unici dei nodi _{1} _{}e pesi _{1} _{}per cui il grado di precisione sia 2 1. 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. Ricordiamo che per una fissata funzione peso , i relativi nodi di Gauss sono le radici del polinomio ortogonale (rispetto ) di grado e che la relativa formula gaussiana ha grado di precisione 2 1 (cf. [1, p. 270]). 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], [2, p.978-985], [11, 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

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

.

  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.

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.

1.5 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 [13] 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  (7)
_{1}^{1} ^{}^{1} 2 3504023872876032  (8)
_{1}^{1} ^{^{2}}  erf(1) \overline{(} ) 1 4936482656248538  (9)
_{1}^{1} 1/(1 16 ^{2})  1/2 atan(4) 0 66290883183401628  (10)
_{1}^{1} ^{1/ ^{2}}  2 ^{1} 2 ^{1} erf(1) ^{1/2} ^{1} 2 ^{1/2} (11)
0 17814771178156086  (12)
_{1}^{1}  ^{3}  1/2  (13)
_{1}^{1} \overline{}2/3  (14)
_{1}^{1} ^{}\overline{(}1 )  2^{1/2} ^{1} 1/2  ^{1/2} erf(2^{1/2})  (15)
1 7791436546919095  (16)

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

1.6 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. [16], [13], [14]). Usiamo una implementazione di Waldvogel [14]

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)  (17)
_{1}^{1} ^{}(a destra, prima riga)  (18)
_{1}^{1} ^{^{2}}  (a sinistra, seconda riga)  (19)
_{1}^{1} 1/(1 16 ^{2})  (a destra, seconda riga)  (20)
_{1}^{1} ^{1/ ^{2}}  (a sinistra, terza riga)  (21)
_{1}^{1}  ^{3}  (a destra, terza riga)  (22)

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

1.7 Un’integrale con funzione peso

Consideriamo ora l’integrale

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

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)

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

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

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

Dal grafico in figura 23, 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,...
      a,b,inline('exp(x).*sqrt(1-x.^2)'));
>> fprintf('\n \t [GAUSS-LEGENDRE]: %15.20f',I_jac);
          [GAUSS-LEGENDRE]: 1.77661983094052260000
>> 1.7791436546919097925911790299941-1.77661983094052260000
ans =
   0.00252382375139
>> 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).

2 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.77914365469190971503octave:31>
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.^2)'));
octave:39> fprintf('\n \t [GAUSS-LEGENDRE]: %15.20f',I_jac)
         [GAUSS-LEGENDRE]: 1.77661983094052411403
octave:40> 1.7791436546919097925911790299941-1.77661983094052411403
ans = 0.0025238
octave:41> length(x_jac)
ans = 10

3 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 [9]

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.

4 Frasi celebri

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

References

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

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

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

[4] W. Gautschi, personal homepage,

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

[5] Mac Tutor, Cavalieri,

http://www-gap.dcs.st-and.ac.uk/ history/Biographies/Cavalieri.html.

[6] Mac Tutor, Cotes,

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

[7] Mac Tutor, Hermite,

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

[8] Mac Tutor, Jacobi,

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

[9] Mac Tutor, Simpson,

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

[10] The MathWorks Inc., Numerical Computing with Matlab,

http://www.mathworks.com/moler.

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

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

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

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

[15] Wikipedia, Cavalieri,

http://it.wikipedia.org/wiki/Bonaventura_Cavalieri.

[16] Wikipedia, Clenshaw-Curtis quadrature,

http://en.wikipedia.org/wiki/Clenshaw-Curtis_quadrature.

[17] Wikipedia, Fast Fourier Transform,

http://en.wikipedia.org/wiki/Fast_Fourier_transform.

[18] Wikipedia, Gaussian quadrature,

http://en.wikipedia.org/wiki/Gaussian_quadrature.

[19] Wikipedia, Quadratura gaussiana,

http://it.wikipedia.org/wiki/Quadratura_di_Gauss.

[20] Wikipedia, Regola del trapezio,

    http://it.wikipedia.org/wiki/Regola_del_trapezio.

[21] Wikipedia, Regola di Cavalieri-Simpson,

    http://it.wikipedia.org/wiki/Regola_di_Cavalieri-Simpson.