Esercitazione n. 3

 

Programmi

%--------------------------------------------------------

% lsolve.m - Forward substitution per Lx=b

%

function x=lsolve(L,b)

n=size(L,1);

x=zeros(n,1);

for i=1:n,

xi=b(i);

for j=1:i-1,

xi=xi-L(i,j)*x(j);

end

x(i)=xi/L(i,i);

end

%--------------------------------------------------------

% usolve.m - Backward substitution per Ux=b

%

function x=usolve(U,b)

n=size(U,1);

x=zeros(n,1);

for i=n:-1:1,

xi=b(i);

for j=i+1:n,

xi=xi-U(i,j)*x(j);

end

x(i)=xi/U(i,i);

end

%--------------------------------------------------------

% lusolve.m - Risoluzione LUx=Pb

%

function x=lusolve(L,U,P,b)

n=size(L,1);

b=P*b; % permuta il termine noto in accordo a P

c=lsolve(L,b);

x=usolve(U,c);

%--------------------------------------------------------

% mylu.m - Fattorizzazione PA=LU

%

function [L,U,P]=mylu(A,piv)

n=size(A,1);

p=zeros(1,n-1); % preallocazione di P

%

if strcmp(piv,'nopivot') == 1,

%------------------------------------------------------------------

% Fattorizzazione senza pivot

%------------------------------------------------------------------

for k=1:n-1,

maxa=abs(A(k,k));

if maxa == 0,

display('Matrice non fattorizzabile a passo k='),k

return

end

p(k)=k;

for i=k+1:n, % Calcolo dei moltiplicatori e eliminazione

Aik=A(i,k)/A(k,k);

A(i,k)=Aik; % memorizza il moltiplicatore

for j=k+1:n,

A(i,j)=A(i,j)-Aik*A(k,j);

end

end

end

elseif strcmp(piv,'parziale') == 1,

%------------------------------------------------------------------

% Fattorizzazione con pivot parziale

%------------------------------------------------------------------

for k=1:n-1,

maxa=abs(A(k,k));

riga_max=k;

p(k)=k;

for i=k+1:n, % scelta pivot (parziale)

abs_Aik=abs(A(i,k));

if abs_Aik > maxa,

maxa=abs_Aik;

riga_max=i;

p(k)=riga_max; % memorizzazione per costruzione matrice di permutazione

end

end

if maxa == 0,

display('Matrice non fattorizzabile a passo k='),k

return

end

%

if k ~= riga_max, % Eventuale scambio righe

for j=1:n,

tmp=A(k,j); A(k,j)=A(riga_max,j); A(riga_max,j)=tmp;

end

end

for i=k+1:n, % Calcolo dei moltiplicatori e eliminazione

Aik=A(i,k)/A(k,k);

A(i,k)=Aik; % memorizza il moltiplicatore

for j=k+1:n,

A(i,j)=A(i,j)-Aik*A(k,j);

end

end

end

elseif strcmp(piv,'antipivot')== 1,

%------------------------------------------------------------------

% Fattorizzazione con antipivot

%------------------------------------------------------------------

for k=1:n-1,

mina=abs(A(k,k));

riga_min=k;

p(k)=k;

%

ik=k;

while mina == 0 & ik<n,

ik=ik+1;

mina=abs(A(ik,k));

riga_min=ik;

p(k)=ik;

end

for i=ik+1:n, % scelta antipivot

abs_Aik=abs(A(i,k));

if abs_Aik < mina & abs_Aik~=0,

mina=abs_Aik;

riga_min=i;

p(k)=riga_min; % memorizzazione per costruzione matrice di permutazione

end

end

if mina == 0,

display('Matrice non fattorizzabile a passo k='),k

return

end

%

if k ~= riga_min,% Eventuale scambio righe

for j=1:n,

tmp=A(k,j); A(k,j)=A(riga_min,j); A(riga_min,j)=tmp;

end

end

for i=k+1:n,% Calcolo dei moltiplicatori e eliminazione

Aik=A(i,k)/A(k,k);

A(i,k)=Aik; % memorizza il moltiplicatore

for j=k+1:n,

A(i,j)=A(i,j)-Aik*A(k,j);

end

end

end

%

else

displ('Errore parametro piv')

return

end

%

% Copia della fattorizzazione in L e U

%

U=triu(A,0);

L=tril(A,-1)+eye(n);

%

% Costruzione della matrice di permutazione

%

P=eye(n);

for i=1:n-1,

for j=1:n,

tmp=P(i,j); P(i,j)=P(p(i),j); P(p(i),j)=tmp;

end

end

%--------------------------------------------------------

%--------------------------------------------------------

 

Esercizio n.1

 

for n=1:100

A=rand(n);

xesatta=(1:n)';

b=A*xesatta;

norm_xesatta=norm(xesatta);

[L,U,P]=mylu(A,'nopivot');

xnopivot=lusolve(L,U,P,b);

errore_relativo_nopivot(n)=norm(xesatta-xnopivot)/norm_xesatta;

[L,U,P]=mylu(A,'parziale');

xparziale=lusolve(L,U,P,b);

errore_relativo_parziale(n)=norm(xesatta-xparziale)/norm_xesatta;

[L,U,P]=mylu(A,'antipivot');

xantipivot=lusolve(L,U,P,b);

errore_relativo_antipivot(n)=norm(xesatta-xantipivot)/norm_xesatta;

xmatlab=A\b;

errore_relativo_matlab(n)=norm(xesatta-xmatlab)/norm_xesatta;

end

 

xn=(1:n);

semilogy(xn,errore_relativo_nopivot,'r',xn,errore_relativo_parziale,'g',xn,errore_relativo_antipivot,'b',xn,errore_relativo_matlab,'y')

xlabel('dimensione n')

ylabel('norma 2 errore relativo')

title('nopivot=rosso, parziale=verde, antipivot=blu, matlab=giallo')

 

Esercizio n.2

 

display('precisione di macchina'),eps

for n=1:100

n

A=rand(n);

[Lparziale,Uparziale,Pparziale]=mylu(A,'parziale');

[Lmatlab,Umatlab,Pmatlab]=lu(A);

Ea_P(n)=norm(Pparziale-Pmatlab,inf);

Ea_L(n)=norm(Lparziale-Lmatlab,inf);

Ea_U(n)=norm(Uparziale-Umatlab,inf);

%

Er_L(n)=norm(Lparziale-Lmatlab,inf)/norm(Lmatlab,inf);

Er_U(n)=norm(Uparziale-Umatlab,inf)/norm(Umatlab,inf);

end

 

eps= 2.220446049250313e-016

norm(Ea_P,inf)

0

norm(Ea_L,inf),norm(Ea_U,inf)

1.309421321371573e-013 5.463407504180395e-013

norm(Er_L,inf),norm(Er_U,inf)

4.340831379544653e-015 8.056932543180377e-015

 

xn=(1:n);

semilogy(xn,Er_L,'r',xn,Er_U,'b')

xlabel('dimensione n')

ylabel('norma infinito errore relativo')

title('L=rosso, U=blu')