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.  1/ 20

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.  1/ 20

Codice complessita’: algoritmo 1

Salviamo il seguente codice in radicesecgrado.m.

p=1000; q =0.018000000081; s o l =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 ;
Alvise Sommariva  Matlab: complessità e stabilità degli algoritmi.  2/ 20

Codice complessita’: 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 complessita’: 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 ( s o l ) > 0 & ( s o l ˜ = 0)
r e r r 1 =abs ( s1−s o l )/ abs ( s o l ) ;
r e r r 2=abs ( t1−s o l )/ abs ( s o l ) ;
f p r i n t f ( ’\ n \ t [ REL . ERR . ALG . 1 ] : %2.2e ’ , r e r r 1 ) ;
f p r i n t f ( ’\ n \ t [ REL . ERR . ALG . 2 ] : %2.2e ’ , r e r r 2 ) ;

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:

>> r a d i c e s e c g r a d o
[ 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.  6/ 20

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));
f p r i n t f ( ’\ n \ t [ SEQ . 1 ] [ INDEX ] : %3.0 f ’ , n ) ;
f p r i n t f ( ’ [ REL . ERR ] : %2.2e ’ , abs ( u ( n+1)− p i )/ p i ) ;

end

r e l e r r u=abs ( u−p i )/ p i ;
f p r i n t f ( ’\ n ’ ) ;
Alvise Sommariva  Matlab: complessità e stabilità degli algoritmi.  8/ 20

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 ; i n n e r s q r t=s q r t (1−c ) ;
z ( n+1)=(2ˆ (n −0.5))∗ s q r t ( 1− i n n e r s q r t ) ;
f p r i n t f ( ’\ n \ t [ SEQ . 2 ] [ N ] : %3.0 f ’ , n ) ;
f p r i n t f ( ’ [ REL . ERR ] : %2.2e ’ , abs ( z ( n+1)− p i )/ p i ) ;

end

r e l e r r 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 ;
i n n e r s q r t=s q r t (1−c ) ;
den=s q r t ( 1+ i n n e r s q r t ) ;
y ( n+1)=num/ den ;
r e l e r r 3=
f p r i n t f ( ’\ n \ t [ SEQ . 3 ] [ N ] : %3.0 f ’ , n ) ;
f p r i n t f ( ’ [ REL . ERR] : % 2 . 2 e ’ , abs ( y ( n+1)− p i )/ p i ) ;

end

r e l e r r y=abs ( y−p i )/ p i ;

Calcolo π: plots

% SEMILOGY PLOT.
hold on ;
s e m i l o g y ( 1 : l e n g t h ( u ) , r e l e r r u , ’ k . ’ ) ;
s e m i l o g y ( 1 : l e n g t h ( z ) , r e l e r r z , ’m+ ’);
s e m i l o g y ( 1 : l e n g t h ( y ) , r e l e r r y , ’ ro ’ ) ;
hold o f f ;

Di seguito digitiamo sulla shell di Matlab/Octave

>> p i g r e c o

Plot risultati

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

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

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

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

▶ 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.  14/ 20

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.  14/ 20

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.  14/ 20

Problema.

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

▶ mediante la successione in avanti

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

▶ mediante la successione all’indietro

\{ \begin{matrix}t_{1000} = 0 \\ t_{n−1} = \frac{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.  15/ 20

Problema.

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

▶ mediante la successione in avanti

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

▶ mediante la successione all’indietro

\{ \begin{matrix}t_{1000} = 0 \\ t_{n−1} = \frac{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.  15/ 20

Successioni ricorrente in Matlab

Scriviamo il codice in un file succricorrente.m.

% SUCCESSIONE RICORRENTE .
% 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=500; M=2∗m;
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

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

Successioni ricorrente in Matlab

% PLOT SEMI−LOGARITMICO .
s u b p l o t ( 2 , 1 , 1 )
s e m i l o g y ( 1 : l e n g t h ( s ) , abs ( s ) , ’ k − ’); hold on ;
s e m i l o g y ( 1 : l e n g t h ( s ) , abs ( t ( 1 : l e n g t h ( s ) ) ) , ’m− ’);
hold o f f ;

Di seguito digitiamo sulla shell di Matlab/Octave

>> s u c c r i c o r r e n t e

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.  18/ 20

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.  19/ 20

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