Quadratura numerica.

Alvise Sommariva

Università degli Studi di Padova Dipartimento di Matematica Pura e Applicata

25 novembre 2011

Quadratura numerica

Il problema della quadratura numerica consiste nell’approssimare 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

con

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

dove w è una funzione peso in un intervallo avente estremi a, b. A seconda degli esempi, l’intervallo I può essere aperto o chiuso.

I termini w_{i} e x_{i} ∈ I sono detti rispettivamente pesi e nodi.

Formule interpolatorie

Sia w : (a, b) → ℝ una funzione non negativa, con (a, b) non necessariamente limitato, tale che

1. ∫_{a}^{b} |x|^{n} w(x) dx < +∞ per tutti gli n ∈ ℕ;

2. ∫_{a}^{b} g(x) w(x) dx = 0 per qualche funzione continua e non negativa g implica g ≡ 0 in (a, b).

La funzione w è detta peso. Le funzioni peso più comuni sono 1. w(x) = 1 con x ∈ [−1, 1] (peso di Legendre);

2. w(x) = \frac{1}{\sqrt{1−x^{2}}} con x ∈ (−1, 1) (peso di Chebyshev);

3. w(x) = (1 − x^{2})^{γ−(1/2)} con x ∈ (−1, 1), γ > (−1/2) (peso di Gegenbauer);

4. w(x) = (1 − x)^{α} · (1 + x)^{β} con x ∈ (−1, 1), α > −1, β > −1 (peso di Jacobi);

5. w(x) = exp (−x) con x ∈ (0, +∞) (peso di Laguerre);

6. w(x) = exp (−x^{2}) con x ∈ (−∞, +∞) (peso di Hermite);

Alvise Sommariva  Quadratura numerica.  3/ 37

Supponiamo sia I = (a, b) (o I = [a, b]) con a < b finiti, x_{1}, …, x_{N} un insieme di N punti a 2 a 2 distinti di I ed f ∈ C(I). Allora esiste finito I^{(w)}(f ). Infatti se I è limitato, per il teorema di Weierstrass e l’integrabilità della funzione w, abbiamo

|I^{(w)}(f )| := |∫_{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 pol. che interpola (x_{i}, f (x_{i})) con i = 1, …, N, (L_{i} è l’i-simo polinomio di Lagrange)

I^{(w)}(f ) ≈ ∫_{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})  (1)

per cui Q_{N}^{(w)}(f ) := ∑_{k=1}^{N} w_{i}f (x_{i}) con w_{i} = ∫_{a}^{b} L_{i}(x) w(x)dx.

Alvise Sommariva  Quadratura numerica.  4/ 37

Definizione. Una formula di quadratura

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

per cui

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

si dice interpolatoria.

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.

Formule di Newton-Cotes

Si supponga [a, b] un intervallo chiuso e limitato di ℝ e sia 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 che si ottengono integrando l’interpolante di f in nodi equispaziati

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

Esempi:

▶ N=2, regola del trapezio, g.d.p.=1:

I (f ) ≈ S_{1}(f ) := S_{1}(f , a, b) := \frac{(b−a) (f (a)+f (b))}{2} ▶ N=3, regola di Cavalieri-Simpson, g.d.p.=3:

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

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

Altre formule di Newton-Cotes Errori formule di Newton-Cotes

Denotiamo con ℙ_{k} lo spazio dei polinomi di grado k. Allora valgono le seguenti formule d’errore:

▶ regola del trapezio, g.d.p.=1:

E_{1}(f ) := I (f ) − S_{1}(f ) = \frac{−h^{3}}{12} f^{(2)}(ξ) per qualche ξ ∈ (a, b). ▶ regola di Cavalieri-Simpson, g.d.p.=3:

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

Si noti il legame tra formule dell’errore e g.d.p.:

▶ nel caso della formula dei trapezi, se f ∈ ℙ_{2}, allora potrebbe

essere essere f^{(2)}(ξ) ≠ 0, mentre se f ∈ ℙ_{1}, allora f^{(2)}(ξ) = 0; ▶ nel caso della formula di Cavalieri-Simpson, se f ∈ ℙ_{4}, allora

potrebbe essere essere f^{(4)}(ξ) ≠ 0, mentre se f ∈ ℙ_{3}, allora f^{(4)}(ξ) = 0.

Formule di Newton-Cotes composte

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

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

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

Formule di Newton-Cotes

Figura: Formula dei trapezi composta ∫_{0.5}^{2} sin (x)dx.

Alvise Sommariva  Quadratura numerica.  11/ 37

Formule dei trapezi composte

Formula composta dei trapezi: fissati il numero N di subintervalli e i punti x_{k} = a + kh dove h = \frac{b−a}{N} è definita da

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

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.  12/ 37

Formule di Cavalieri-Simpson 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})]

(6)

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

Formule gaussiane

Nelle formule interpolatorie di Newton-Cotes (come ad esempio la regola del Trapezio o di Cavalieri-Simpson) i nodi x_{1}, …, x_{n} sono equispaziati e il grado di precisione δ è generalmente uguale almeno a n − 1 ma in alcuni casi, come per la regola di Cavalieri-Simpson, uguale al numero di nodi n. Vediamo ora formule ∑_{k=1}^{n} w_{k}f (x_{k})

▶ hanno grado di precisione δ maggiore di 2n − 1;

▶ per qualche funzione f : (a, b) → ℝ continua approssimino I^{(w)}(f ) := ∫_{a}^{b} f (x)w(x) dx, con (a,b) non necessariamente limitato.

Formule gaussiane

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.

Importante: Si dimostra che gli zeri sono semplici e interni all’intervallo (a,b).

Vantaggi integrali con peso

Diamo di seguito l’idea sul perchè si utilizzino funzioni peso w per il calcolo di integrali. Se dobbiamo approssimare ∫_{a}^{b} g(x) dx con g regolare allora \underline{di solito} non sarà difficile calcolare numericamente il valore dell’integrale, mentre se g è poco regolare sarà molto più oneroso ottenere buoni risultati.

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

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

La regolarità di f e la particolarità del calcolo dipendentemente dalla funzione peso permette di approssimare meglio ∫_{a}^{b} g(x) dx = ∫_{a}^{b} f (x)w(x) dx.

Vantaggi integrali con peso

Teorema. Sia (a, b) un intervallo limitato e si supponga w : (a, b) → ℝ sia una funzione peso, I^{(w)} = ∫_{a}^{b} f (x)w(x)dx e I_{m}^{(w)}(f ) := ∑_{k=1}^{m} w_{i}f (x_{i}) una formula di quadratura avente g.d.p. n. Allora ricordato che ∥I_{n}∥_{∞} = sup_{f ∈C(a,b)\0} |I_{n}(f )|/∥f ∥_{∞}, ∥w∥_{1} = ∫_{a}^{b} w(x)dx, abbiamo

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

Quindi è importante avere formule con grado di precisione n alto e f regolare così da avere minore min_{q_{n}∈ℙ_{n}} ∥f − q_{n}∥_{∞}.

Esempio

Consideriamo il calcolo approssimato

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

(7)

▶ Non è suggerito il calcolo dell’integrale con formule composte

poichè la funzione exp (x)\sqrt{1 − x} è poco regolare. ▶ Non è indicato usare formule gaussiane con peso w(x) = 1

poichè la funzione exp (x)\sqrt{1 − x} è poco regolare. ▶ Si suggerisce di usare formule gaussiane con peso di Jacobi

w(x) = \sqrt{1 − x} = (1 − x)^{1/2} poichè la funzione exp (x) è C^{∞}(−1, 1).

Formula trapezi composta

La funzione trapezi composta appena esposta calcola i nodi e i pesi della omonima formula composta.

 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 ;

Formula Cavalieri-Simpson composta

La funzione simpson composta appena esposta calcola i nodi e i pesi della omonima formula 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 ;

Vogliamo approssimare l’integrale ∫_{−1}^{1} x^{20}dx = 2/21 ≈ 0.09523809523810 tramite le formule composte dei trapezi e Cavalieri-Simpson. Perchè il paragone sia veritiero, decidiamo che il numero di valutazioni della funzione integranda sia uguale. Scriviamo nel file demo composte.m

N=11; %SCEGLIERE DISPARI .
a=−1; b=1; f=inline ( ’ x .ˆ 20 ’ ) ;
% a=0; b=1; f=i n l i n e ( ’ exp ( x ) ’ ) ;
N_trap=N−1; % TRAPEZI COMPOSTA.

[ x_trap , w_trap]= trapezi_composta ( N_trap , a , b ) ;
fx_trap=f e v a l ( f , x_trap ) ; I_trap=w_trap ’ ∗ fx_trap ;

N_simpson=(N−1) / 2 ; % CAV. SIMPSON COMPOSTA.

[ x_simp , w_simp]= simpson_composta ( N_simpson , a , b ) ;
fx_simp=f e v a l ( f , x_simp ) ; I_simp=w_simp ’ ∗ fx_simp ;

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 [ CS .COMP . ] [ PTS ] : %4.0 f ’ , l e n g t h ( x_simp ) ) ;
f p r i n t f ( ’ \n\ t [ CS .COMP . ] [ RIS ] : %14.14 f \n\n ’ , I_simp ) ;

Ricordando che il risultato è 0.09523809523810, otteniamo i poco soddisfacenti

[ TRAPEZI COMPOSTA ] [ PTS ] :  11

[ TRAPEZI COMPOSTA ] [ RIS ] : 0.20462631505024

[ SIMPSON COMPOSTA ] [ PTS ] :  11

[ SIMPSON COMPOSTA ] [ RIS ] : 0.13949200364447

Posto N = 51 nella prima riga di demo composte.m

[ TPZ . COMP . ] [ PTS ] :  51

[ TPZ . COMP . ] [ RIS ] : 0.10052328836742

[ CS . COMP . ] [ PTS ] :  51

[ CS . COMP . ] [ RIS ] : 0.09542292188917

L’errore rel. della formula comp. dei trapezi è solo di 5.55 · 10^{−2} mentre per quella di Cav. Simpson composta 1.94 · 10^{−3}, mentre gli err. assoluti sono risp. circa 5.28 · 10^{−3} e 1.85 · 10^{−4}.

Alvise Sommariva  Quadratura numerica.  22/ 37

Approssimiamo ∫_{0}^{1} exp (x)dx = exp(1) − 1 ≈ 1.718281828459046 tramite le formule composte dei trapezi e Cavalieri-Simpson. Perchè il paragone sia veritiero, decidiamo che il numero di valutazioni della funzione integranda sia uguale. Modifichiamo la funzione f nel file demo composte.m e ricaviamo

[ TPZ . COMP . ] [ PTS ] :  11

[ TPZ . COMP . ] [ RIS ] : 1.71971349138931

[ CS . COMP . ] [ PTS ] :  11

[ CS . COMP . ] [ RIS ] : 1.71828278192482

Gli errori relativi sono rispettivamente circa 8.33 · 10^{−4} e 5.55 · 10^{−7}, quelli assoluti circa 1.43 · 10^{−3} e 9.55 · 10^{−7}.

Esempio 2: errore

Nel caso della formula dei trapezi composta, abbiamo

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) e quindi

|E_{1}^{(c)}(f )| ≤ |\frac{−(b − a)}{12}h^{2} max_{x∈(0,1)}exp(x)| ≤ \frac{1}{12}h^{2}exp(1)

Nel nostro esempio h = 0.2 e quindi

|E_{1}^{(c)}(f )| ≤ \frac{1}{12}h^{2}exp(1) = \frac{1}{12}0.2^{2}exp(1) < 9.07 · 10^{−3}, mentre l’errore assoluto effettivo era circa 1.43 · 10^{−3}.

Alvise Sommariva  Quadratura numerica.  24/ 37

Esempio 2: errore

Nel caso della formula di Cavalieri-Simpson composta, abbiamo

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

|E_{3}^{(c)}(f )| ≤ |\frac{−(b − a)}{180} (\frac{h}{2})^{4} max_{x∈(0,1)}exp(x)| ≤ \frac{1}{180} (\frac{h}{2})^{4} exp(1)

Nel nostro esempio h = 0.4 e quindi

|E_{3}^{(c)}(f )| ≤ \frac{1}{180} (\frac{h}{2})^{4} exp(1) = \frac{1}{180} (\frac{0.4}{2})^{4} exp(1) < 2.42 · 10^{−5},

mentre l’errore assoluto effettivo era circa 9.55 · 10^{−7}.

Alvise Sommariva  Quadratura numerica.  25/ 37

Formule gaussiane in Matlab

Il calcolo di nodi e pesi delle formule gaussiane, per una fissata funzione peso w non è semplice. Di seguito ci interesseremo al caso particolare delle funzioni peso di Jacobi. A tal proposito aiutano le routines gauss.m e r jacobi.m di W. Gautschi e D. Laurie che possono essere scaricate dalla homepage del primo autore. La funzione gauss jacobi.m, se registrata nella stessa cartella di gauss.m e r jacobi.m, calcola i nodi x = {x_{i}}_{i=1,...,n} e i pesi w = {w_{i}}_{i=1,...,n} della formula gaussiana (avente g.d.p. 2n − 1)

∫_{−1}^{1} f (x) · (1 − x)^{α}(1 + x)^{β} ≈ ∑_{i=1}^{n}w_{i}f (x_{i})

f u n c t i o n [ x , w]=gauss_jacobi ( n , alpha , beta )
ab=r_jacobi ( n , alpha , beta ) ;
xw=gauss ( n , ab ) ;

x=xw ( : , 1 ) ;
w=xw ( : , 2 ) ;

Formule gaussiane in Matlab: scalatura

Specialmente per le formule con peso w(x) ≡ 1 si ha spesso da effettuare l’integrale ∫_{a}^{b} f (x)dx con (a, b) diverso da (−1, 1). In tal caso si osserva che se γ(t) = (a(1 − t)/2) + (b(1 + t)/2), per {t_{k}} e {w_{k}} nodi e pesi della formula di Gauss-Legendre,

∫_{a}^{b} f (x)dx = ∫_{−1}^{1} \frac{b − a}{2}f (γ(t))dt ≈ ∑_{k=1}^{n} \frac{b − a}{2}w_{k}f (γ(t_{k}))

e quindi

∫_{a}^{b} f (x)dx ≈ ∑_{k=1}^{n}w_{k}^{∗}f (x_{k}^{∗})

con w_{k}^{∗} = \frac{b−a}{2}w_{k}, x_{k} = γ(t_{k}) = (a(1 − t_{k})/2) + (b(1 + t_{k})/2). Questo processo si chiama scalatura dei nodi e dei pesi.

Formule gaussiane in Matlab: demo

Scriviamo quindi nel file demo gauss jacobi

N=11;
alpha =0; beta =0;
a=−1; b=1; f=inline ( ’ x .ˆ 20 ’ ) ;
% a=0; b=1; f=i n l i n e ( ’ exp ( x ) ’ ) ;

[ x , w]=gauss_jacobi ( N , alpha , beta ) ;

xx=a∗(1−x )/2+b∗(1+x ) / 2 ; ww=(b−a ) ∗w / 2 ;
fxx=f e v a l ( f , xx ) ;

I=ww ’ ∗ fxx ;

f p r i n t f ( ’ \n\ t [N] : % 3 . 0 f [ I ] : %1.15 e \n\n ’ , N , I )

calcolando così con una formula gaussiana di 11 nodi l’integrale ∫_{0}^{1} exp (x)dx = exp(1) − 1 ≈ 1.718281828459046.

Alvise Sommariva  Quadratura numerica.  28/ 37

Formule gaussiane in Matlab: demo ed esempio 1

Dalla shell di Matlab/Octave

>> demo_gauss_jacobi

[ N ] : 11 [ I ] : 9.523809523809551 e−02
>> I
I = 9.523809523809551 e−02
>> Iex =2/21; abs ( I−Iex )
ans = 2.775557561562891 e−16

che va confrontato, a parità di nodi, con 5.28 · 10^{−3} e 1.85 · 10^{−4} rispettivamente delle formule dei trapezi e Cavalieri-Simpson composte. Si osservi che il grado di precisione della formula gaussiana con n = 11 è 2n − 1 = 21 e quindi l’integrale ∫_{−1}^{1} x^{20}dx è calcolato, a meno di errori macchina, esattamente.

Dalla shell di Matlab/Octave, rimaneggiando i commenti del codice,

>> demo_gauss_jacobi

[ N ] : 11 [ I ] : 1.718281828459045 e+00
>> Iex=exp ( 1 )−exp ( 0 ) ; I =1.718281828459045 e+00;
>> abs ( Iex−I )
ans = 4.440892098500626 e−16

che va confrontato, a parità di nodi, con 1.43 · 10^{−3} e 9.55 · 10^{−7} rispettivamente delle formule dei trapezi e Cavalieri-Simpson composte.

Approssimiamo in un file Matlab/Octave demo comparison.m ∫_{−1}^{1} exp (x)\sqrt{1 − x} dx = 1.77914365469190979259 …. mediante la formula dei trapezi composta, Cav.-Simpson composta, Gauss-Legendre e Gauss-Jacobi con α = 1/2, β = 0.

N=11; a=−1;b=1; Iex =1.7791436546919097925911790299941;
g=inline ( ’ exp ( x ) . ∗ s q r t (1−x ) ’ ) ; f=inline ( ’ exp ( x ) ’ ) ;

[ xt , wt]= trapezi_composta ( N , a , b ) ; % TPZ .

fxt=f e v a l ( g , xt ) ; It=wt ’ ∗ fxt ; et=abs ( Iex−It ) ;

[ xcs , wcs]= simpson_composta ( N , a , b ) ; % C . S .
fxcs=f e v a l ( g , xcs ) ; Ics=wcs ’ ∗ fxcs ; ecs=abs ( Iex−Ics ) ;
[ xgl , wgl]=gauss_jacobi ( N , 0 , 0 ) ; % G . L .
fxgl=f e v a l ( g , xgl ) ; Igl=wgl ’ ∗ fxgl ; egl=abs ( Iex−Igl ) ;

[ xj , wj]=gauss_jacobi ( N , 0 . 5 , 0 ) ; % G . J .
fxj=f e v a l ( f , xj ) ; Ij=wj ’ ∗ fxj ; ej=abs ( Iex−Ij ) ;
f p r i n t f ( ’ \n TPZ:%1.1 e CS:%1.1 e ’ , et , ecs ) ;
f p r i n t f ( ’GL:%1.1 e GJ:%1.1 e \n\n ’ , egl , ej ) ;

Con ovvia notazione, a parità di valutazioni della funzione integranda, abbiamo

>> demo_comparison

[ TPZ ] : 4 . 4 e−02 [ CS ] : 6 . 2 e−03 [ GL ] : 5 . 3 e−04 [ GJ ] : 6 . 7 e
−16

>>

che evidenzia la superiorità delle formule gaussiane, una volta che il peso è scelto opportunamente.

Cubatura e Metodo Montecarlo

Si parla di cubatura quando l’integrale

I^{w}(f ) = ∫_{Ω} f (x)w(x)dx  (8)

è da calcolarsi in un dominio Ω ⊂ ℝ^{m} con m > 1. Il problema è particolarmente difficile per dimensioni m ≫ 1. In tal caso, un classico approccio consiste nell’applicare il Metodo Montecarlo. Per semplicità tratteremo il caso del cubo m dimensionale Ω = [0, 1]^{m}, e w(x) ≡ 1. Nella sua versione più semplice, il Metodo Montecarlo consiste nell’approssimare (8) con

Q_{N}(f ) = ∑_{k=1}^{N}w_{k}f (x_{k})

dove x_{1}, …, x_{N} sono generati casualmente da una distribuzione uniforme su Ω e w_{k} = 1/N.

Alvise Sommariva  Quadratura numerica.  33/ 37

Si conosce n stima probabilistica dell’errore

|I (f ) − I_{N}(f )| ≈ \sqrt{I (f^{2}) − (I (f ))^{2}}/\sqrt{N}

Si osservi che la convergenza, seppure erratica e lenta, non dipende dalla dimensione m.

I metodiquasi Montecarlo scelgono in modo opportuno i nodi, al fine di avere una performance migliore. Esistono varie distribuzioni di nodi come quelle di Sobol, di Niederreiter, di Fibonacci nonchè le cosidette lattice rules. Per semplicità , paragoneremo il metodo Montecarlo con una di queste distribuzioni (quella di Sobol), per il calcolo di un certo integrale.

Calcoliamo quale test

∫_{0}^{1} ∫_{0}^{1} ∫_{0}^{1} ∫_{0}^{1} x · cos (y) · exp (w) · z^{1}0dxdydwdz

con il Metodo di Montecarlo avente 4^{k} nodi e quello di Sobol, basato su una sequenza con la stessa cardinalità . Dagli esperimenti vedremo la miglior performance della sequenza di Sobol.

Cubatura e Metodo Montecarlo

Salviamo in demo sobol:

Npow =10; N=4ˆ Npow ;

f=inline ( ’ x . ∗ cos ( y ) . ∗ exp (w) . ∗ ( z . ˆ 1 0 ) ’ , ’ x ’ , ’ y ’ , ’w ’ , ’ z ’ ) ;
Iex =(1/2) ∗( s i n ( 1 )−s i n ( 0 ) ) ∗( exp ( 1 ) −1) ∗(1/11) ;
% MONTECARLO.
X=rand ( N , 4 ) ; fX=f e v a l ( f , X ( : , 1 ) , X ( : , 2 ) , X ( : , 3 ) , X ( : , 4 ) ) ;
Imc=sum ( fX ) /N ;
% SOBOL.
P=sobolset ( 4 ) ; X= net ( P , 4 ˆ Npow ) ; fX=f e v a l ( f , X ( : , 1 ) , X
( : , 2 ) , X ( : , 3 ) , X ( : , 4 ) ) ;
Isb=sum ( fX ) /N ;
% STATS .
f p r i n t f ( ’ \n \ t [N pt ] : %1.16 e \n ’ , N ) ;
f p r i n t f ( ’ \n \ t [ I . EX ] : %1.16 e ’ , Iex ) ;
f p r i n t f ( ’ \n \ t [ I .MC] : %1.16 e ’ , Imc ) ;
f p r i n t f ( ’ \n \ t [ I . Sb ] : %1.16 e \n ’ , Isb ) ;
f p r i n t f ( ’ \n \ t [ E .MC] : %1.16 e ’ , abs ( Iex−Imc ) ) ;
f p r i n t f ( ’ \n \ t [ E . Sb ] : %1.16 e \n \n ’ , abs ( Iex−Isb ) ) ;
Alvise Sommariva  Quadratura numerica.  36/ 37

In forma tabulare, posto e_{MC} l’errore assoluto del Metodo Montecarlo, e con e_{QMC} quello di Quasi-Montecarlo con la sequenza di Sobol:

N  e_{MC}  e_{QMC}
64 1.8836e − 02 4.3077e − 03
256 8.0959e − 03 1.6564e − 03
1024 8.2839e − 03 3.7884e − 04
4096 3.9135e − 03 9.6362e − 05
16384 4.6608e − 04 1.7776e − 05
65536 3.2112e − 04 6.8086e − 06
262144 3.2929e − 04 1.6914e − 06
1048576 7.4100e − 05 4.3227e − 07

Alvise Sommariva  Quadratura numerica.  37/ 37