Matlab: complessità e stabilità degli algoritmi. Alcuni esempi.

Alvise Sommariva

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

21 novembre 2011

Stabilità : radici secondo grado

Dato x^{2} + 2 px − q, con \sqrt{p^{2} + q} ≥ 0 eseguiamo un primo algoritmo Matlab che valuta la radice via:

y = −p + \sqrt{p^{2} + q}.  (1)

▶ \sqrt{p^{2} + q} ≥ 0 implica radici reali.

▶ Potenzialmente instabile per p ≫ q a causa della sottrazione tra p e \sqrt{p^{2} + q} (cancellazione).

Valutiamo radice con un secondo algoritmo stabile via razionalizzazione di (1):

y = −p + \sqrt{p^{2} + q} = \frac{(−p + \sqrt{p^{2} + q})(p + \sqrt{p^{2} + q})}{(p + \sqrt{p^{2} + q})}
= \frac{q}{(p + \sqrt{p^{2} + q})}  (2)

Alvise Sommariva  Matlab: complessità e stabilità degli algoritmi.  2/ 24

Codice stabilità : algoritmo 1

Salviamo il seguente codice in radicesecgrado.m.

p=1000; q =0.018000000081; sol =0.9∗10ˆ ( −5) ;
% ALGORITMO 1
s=p ˆ 2 ;

t=s+q ;

i f t >=0
u=s q r t ( t ) ;

e l s e

f p r i n t f ( ’ \n \ t [ RADICI COMPLESSE] ’ ) ;

end

s1=−p+u ;

Codice stabilità : algoritmo 2

% ALGORITMO 2
s=p ˆ 2 ;

t=s+q ;

i f t >=0
u=s q r t ( t ) ;

e l s e

f p r i n t f ( ’ \n \ t [ RADICI COMPLESSE] ’ ) ;

end

v=p+u ;
t1=q/v ;

Codice stabilità : stampa risultati

f p r i n t f ( ’ \n \ t [ ALG . 1 ] : %10.19 f ’ , s1 ) ;
f p r i n t f ( ’ \n \ t [ ALG . 2 ] : %10.19 f ’ , t1 ) ;
i f l e n g t h ( sol ) > 0 & ( sol ˜ = 0)
rerr1 =abs ( s1−sol ) / abs ( sol ) ;
rerr2=abs ( t1−sol ) / abs ( sol ) ;
f p r i n t f ( ’ \n \ t [ REL . ERR . ALG . 1 ] : %2.2 e ’ , rerr1 ) ;
f p r i n t f ( ’ \n \ t [ REL . ERR . ALG . 2 ] : %2.2 e ’ , rerr2 ) ;

end

Test.

Come previsto, il secondo algoritmo si comporta notevolmente meglio del primo, che compie un errore relativo dell’ordine di circa 10^{−9}. Infatti:

>> radicesecgrado

[ ALG . 1 ] [ 1 ] : 0.0000089999999772772
[ ALG . 2 ] [ 1 ] : 0.0000090000000000000
[ REL . ERR . ] [ ALG . 1 ] : 2.52 e−009
[ REL . ERR . ] [ ALG . 2 ] : 0.00 e+000

>>

Calcolo π

Eseguiamo un codice Matlab che valuti le successioni {u_{n}}, {z_{n}}, definite rispettivamente come

s_{1} = 1, s_{2} = 1 + \frac{1}{4}

u_{1} = 1, u_{2} = 1 + \frac{1}{4}
s_{n+1} = s_{n} + \frac{1}{(n+1)^{2}}
u_{n+1} = \sqrt{6 s_{n+1}}

e

\{ \begin{matrix}z_{1} = 1, z_{2} = 2 \\ z_{n+1} = 2^{n− \frac{1}{2}} \sqrt{1 − \sqrt{1 − 4^{1−n} · z_{n}^{2}}} \end{matrix}  (3)

che teoricamente convergono a π.

Alvise Sommariva  Matlab: complessità e stabilità degli algoritmi.  7/ 24

Calcolo π

Implementiamo poi la successione, diciamo {y_{n}}, che si ottiene razionalizzando (3), cioè moltiplicando numeratore e denominatore di

z_{n+1} = 2^{n− \frac{1}{2}} \sqrt{1 − \sqrt{1 − 4^{1−n} · z_{n}^{2}}}

per

\sqrt{1 + \sqrt{1 − 4^{1−n} · z_{n}^{2}}}

e calcoliamo u_{m}, z_{m} e y_{m} per m = 2, 3, …, 40 (che teoricamente dovrebbero approssimare π). Infine disegniamo in un unico grafico l’andamento dell’errore relativo di u_{n}, z_{n} e y_{n} rispetto a π aiutandoci con l’help di Matlab relativo al comando semilogy.

Calcolo π: metodo 1

In seguito scriviamo un’implementazione di quanto richiesto commentando i risultati. Si salvi in un file pigreco.m il codice

% SEQUENZE CONVERGENTI ” PI GRECO” .
% METODO 1 .
s ( 1 ) =1; u ( 1 ) =1;
s ( 2 ) =1.25; u ( 2 )=s ( 2 ) ;
f o r n=2:40
s ( n+1)=s ( n )+(n+1)ˆ ( −2) ;
u ( n+1)=s q r t (6∗s ( n+1) ) ;

end

rel_err_u=abs ( u−p i ) / p i ;
f p r i n t f ( ’ \n ’ ) ;

Calcolo π: metodo 2

% METODO 2 .
format long
z ( 1 ) =1;
z ( 2 ) =2;
f o r n=2:40
c=(4ˆ (1−n ) ) ∗ ( z ( n ) ) ˆ 2 ; inner_sqrt=s q r t (1−c ) ;
z ( n+1)=(2ˆ (n −0.5) ) ∗ s q r t ( 1−inner_sqrt ) ;

end

rel_err_z=abs ( z−p i ) / p i ;
f p r i n t f ( ’ \n ’ ) ;

Calcolo π: metodo 3

% METODO 3 .
y ( 1 ) =1;
y ( 2 ) =2;
f o r n=2:40
num =(2ˆ (1/2) ) ∗ abs ( y ( n ) ) ;
c=(4ˆ (1−n ) ) ∗ ( z ( n ) ) ˆ 2 ;

     inner_sqrt=s q r t (1−c ) ;
     den=s q r t ( 1+inner_sqrt ) ;
     y ( n+1)=num/den ;

end

rel_err_y=abs ( y−p i ) / p i ;

Calcolo π: plots

% SEMILOGY PLOT.
s e m i l o g y ( 1 : l e n g t h ( u ) , rel_err_u , ’ k . ’ , …
1 : l e n g t h ( z ) , rel_err_z , ’m+’ , …
1 : l e n g t h ( y ) , rel_err_y , ’ ro ’ ) ;

Di seguito digitiamo sulla shell di Matlab/Octave

>> pigreco

Plot risultati

Figura: Grafico che illustra le 3 successioni, rappresentate rispettivamente da ., + e o.

Alvise Sommariva  Matlab: complessità e stabilità degli algoritmi.  13/ 24

Discussione risultati.

▶ La prima successione converge molto lentamente a π, la seconda diverge mentre la terza converge velocemente a π.

▶ Per alcuni valori {z_{n}} e {y_{n}} coincidono per alcune iterazioni per poi rispettivamente divergere e convergere a π. Tutto ciò è naturale poichè le due sequenze sono analiticamente (ma non numericamente) equivalenti.

▶ Dal grafico dell’errore relativo, la terza successione, dopo aver raggiunto errori relativi prossimi alla precisione di macchina, si assesta ad un errore relativo di circa 10^{−15} (probabilmente per questioni di arrotondamento).

Una successione ricorrente.

Consideriamo la successione {I_{n}} definita da

I_{n} = e^{−1} ∫_{0}^{1} x^{n} e^{x} dx  (4)

▶ n = 0: I_{0} = e^{−1} ∫_{0}^{1} e^{x} dx = e^{−1}(e^{1} − 1). ▶ integrando per parti

I_{n+1} = e^{−1} (x^{n+1} e^{x} |_{0}^{1} −(n + 1) ∫_{0}^{1} x^{n} e^{x} dx)
= 1 − (n + 1) I_{n}.

▶ I_{n} > 0, decrescente e si prova che I_{n} → 0 come 1/n.

Alvise Sommariva  Matlab: complessità e stabilità degli algoritmi.  15/ 24

Problema.

Calcoliamo I_{n} per n = 1, …, 99:

▶ mediante la successione in avanti

\{ \begin{matrix}I_{0} = e^{−1}(e^{1} − 1) \\ I_{n+1} = 1 − (n + 1) I_{n}. \end{matrix}  (5)

▶ mediante la successione all’indietro

\{ \begin{matrix}t_{1000} = 0 \\ t_{n−1} = (1 − t_{n})/n. \end{matrix}

Si noti che se I_{n+1} = 1 − (n + 1) I_{n} allora I_{n} = (1 − I_{n+1})/(n + 1) e quindi I_{n−1} = (1 − I_{n})/n.

Alvise Sommariva  Matlab: complessità e stabilità degli algoritmi.  16/ 24

Successioni ricorrente in Matlab

Scriviamo il codice in un file succricorrente.m.

% SUCCESSIONE RICORRENTE .
c l e a r a l l ;
% SUCCESSIONE ” s n ” .
s ( 1 )=exp (−1) ;
f o r n=1:99
s ( n+1)=1−(n+1)∗s ( n ) ;

end

% SUCCESSIONE ” t n ” .
M=1000;
t=z e r o s ( M , 1 ) ; % INIZIALIZZAZIONE ” t ” .
f o r n=M : −1:2
j=n−1;
t ( j )=(1−t ( n ) ) /n ;

end

Successioni ricorrente in Matlab

% PLOT SEMI−LOGARITMICO .
c l f ;
s e m i l o g y ( 1 : l e n g t h ( s ) , abs ( s ) , ’ k−’ , …
1 : l e n g t h ( s ) , abs ( t ( 1 : l e n g t h ( s ) ) ) ) ;

Di seguito digitiamo sulla shell di Matlab/Octave

>> succricorrente

Plot risultati

Figura: Grafico che illustra i valori assoluti assunti dalla successione in avanti (in nero) e all’indietro (in rosa magenta).

Alvise Sommariva  Matlab: complessità e stabilità degli algoritmi.  19/ 24

Commento successione ricorrente.

La successione in avanti amplifica gli errori. Infatti se

\overline{I_{0}} = I_{0} + ϵ

\overline{I}_{n+1} = 1 − (n + 1)\overline{I}_{n}

allora

\overline{I}_{1} = 1 − \overline{I}_{0} = 1 − (I_{0} + ϵ) = I_{1} − ϵ
\overline{I}_{2} = 1 − 2\overline{I}_{1} = 1 − 2(I_{1} − ϵ) = I_{2} + 2 · 1 · ϵ
\overline{I}_{3} = 1 − 3\overline{I}_{2} = 1 − 3(I_{2} + 2ϵ) = I_{3} − 3 · 2 · 1 · ϵ
\overline{I}_{4} = 1 − 4\overline{I}_{3} = 1 − 4(I_{3} + 3ϵ) = I_{4} − 4 · 3 · 2 · 1 · ϵ

e in generale

\overline{I}_{n} = I_{n} + (−1)^{n}n! · ϵ.

con il termine n! · ϵ che tende velocemente a +∞ al crescere di n.

Alvise Sommariva  Matlab: complessità e stabilità degli algoritmi.  20/ 24

Commento successione ricorrente.

La successione all’indietro smorza gli errori. Infatti, se

\overline{I}_{m} = I_{m} + ϵ

e

\overline{I}_{n−1} = (1 − \overline{I}_{n})/n

allora

\overline{I}_{m−1} = (1 − \overline{I}_{m})/m = (1 − (I_{m} + ϵ))/m = I_{m} − ϵ/m
\overline{I}_{m−2} = (1 − (I_{m−1} + ϵ/m))/(m − 1) = I_{m−2} − ϵ/((m − 1) · m)
…
\overline{I}_{m−k} = … = I_{m−k} − ϵ/ ∏_{s=0}^{k}(m − s)

con il termine ϵ/ ∏_{s=0}^{k}(m − s) che tende velocemente a 0 al crescere di k.

Alvise Sommariva  Matlab: complessità e stabilità degli algoritmi.  21/ 24

Complessità : algoritmo di Horner.

Ci poniamo il problema di valutare il polinomio

p(x) = a_{0} + a_{1} · x + … + a_{n} · x^{n}  (6)

in un punto x. Osserviamo che

p(x) = a_{0} + x · (a_{1} + x · (a_{2} + … + x · (a_{n−1} + x · a_{n}))) (7)

Supponiamo sia a = (a_{0}, …, a_{n}) il vettore di dimensione n + 1 delle componenti del polinomio. Possiamo valutare il polinomio tramite i seguenti due algoritmi, il primo che valuta direttamente il polinomio secondo quanto descritto in (6), il secondo che effettua la stessa operazione come descritto in (7) calcolando dapprima s_{1} = a_{n−1} + x · a_{n}, poi s_{2} = a_{n−2} + x · s_{1} e così via.

Complessità : algoritmo di Horner.

In Matlab avremo allora

f u n c t i o n s=algoritmo1 ( a , x )

xk=1; s=a ( 1 ) ;
f o r i=2: l e n g t h ( a )

   xk=xk∗x ;
   s=s+a ( i ) ∗xk ;

end

e

f u n c t i o n s=algoritmo2 ( a , x )

L=l e n g t h ( a ) ;
s=a ( L ) ; % COMPONENTE a n IMMAGAZZINATA IN a ( n+1) .
f o r i=L−1: −1:1
s=a ( i )+x∗s ;

end

Complessità : algoritmo di Horner.

Se lanciamo il codice demo horner per la valutazione di p(x) = 1 + 2 · x + 3 · x^{2} + 4 · x^{3} in x = π

c l e a r a l l ;
a=[1 2 3 4 ] ;
x=p i ;

y1=algoritmo1 ( a , x ) ;
y2=algoritmo2 ( a , x ) ;

format long ;

y1
y2

otteniamo

>> demo_horner

ans = 1.609171052316469 e+02
y1 = 1.609171052316469 e+02
y2 =1.609171052316469 e+02

>>

Alvise Sommariva  Matlab: complessità e stabilità degli algoritmi.  24/ 24

Complessità : algoritmo di Horner.

La differenza sta nella complessità computazionale e non nel risultato numerico. Il primo codice richiede 2n moltiplicazioni e n somme, mentre il secondo algoritmo n moltiplicazioni e n somme.