10ma Clase Métodos Numéricos

estamos comparando:
1) Tiempos de ejecución
2) Error

de distintos métodos de resolución de problemas matriciales tipo
AX=b

probamos con matrices random de 2000x2000 y 2000x1

Link a carpeta con funciones usadas 

tabla:

	SVD	LU	INV	QR
Tiempo de ejecución	76.888625	0.350403	1.573359	3.20302
Error	4.06E-12	4.49E-11	1.24E-10	2.73E-12

>> A=randn(2000,2000);
>> b=randn(2000,1);

SVD:

Función:
_____________________________

  function X=resuelve_svd(A,b)
[U,S,V]=svd(A);
Z=U'*b;
W=diag(1./diag(S))*Z;
X=V*W;
_________________________

>> tic
X=resuelve_svd(A,b);
toc
norm(A*X-b)
Elapsed time is 76.888625 seconds.

ans =

  4.0640e-012

>>

LU:

Función:
_______________________

function X=resuelve_LU(A,b)
[L,U,P]=lu(A);
Z=resuelve_ti(L,P*b);
X=resuelve_ts(U,Z);
_______________________
aux1:

function Z=resuelve_ti(L,C)
Z=zeros(size(C));
Z(1)=C(1)/L(1,1);
for k=2:length(Z)
    Z(k)=(C(k)-L(k,1:k-1)*Z(1:k-1))/L(k,k);
end
____________________
aux2:

function X=resuelve_ts(U,z)
n=length(z);
X=zeros(size(z));
X(n)=z(n)/U(n,n);
for k=n-1:-1:1
    X(k)=(z(k)-U(k,k+1:n)*X(k+1:n))/U(k,k);
end





Resuelto:

>> tic
X=resuelve_LU(A,b);
toc
norm(A*X-b)
Elapsed time is 0.350403 seconds.

ans =

  4.4870e-011

>>




INV:

Función:

________________

function X=resuelve_inv(A,b)
X=(A^(-1))*b;

_______________


resuelto:
>> tic
X=resuelve_inv(A,b);
toc
norm(A*X-b)
Elapsed time is 1.573359 seconds.

ans =

  1.2354e-010

>>



La conclusión es que INV no sirve, SVD da poco error pero tarda mucho, y LU tarda poco pero da error.

Por eso presentaremos una nueva descomposición que tiene menos error que LU y es mas rapida que SVD.

Es la descomposición QR

 Función:
_________________

function X=resuelve_qr(A,b)
[Q,R,P]=qr(A);
X=P*resuelve_ts(R,Q'*b);

________________

resuelto:

>> tic
X=resuelve_qr(A,b);
toc
norm(A*X-b)
Elapsed time is 3.203020 seconds.

ans =

  2.7314e-012

>> 



En el caso de matrices definidas positivas, simétricas (Mij=Mji) y reales.


Descomposición de Cholesky


Ejercicio:

A)
Programar una función que resuelva usando cholesky


______________

function X=resuelve_chol(A,b)
R=chol(A);
Z=resuelve_ti(R',b);
X=resuelve_ts(R,Z);

_____________
B)
Probamos.

Para eso necesitamos matris simetrica positiva:
M=A'*A;

resolver con cholesky y con LU para comparar.



tic
X=resuelve_chol(M,b);
toc
norm(M*X-b)
Elapsed time is 0.266135 seconds.

ans =

  1.0242e-008

>> tic
X=resuelve_LU(M,b);
toc
norm(M*X-b)
Elapsed time is 0.281710 seconds.

ans =

  9.5414e-009


_____________
TABLA:

	LU	Chol
Tiempo de ejecución	0.28171	0.266135
Error	9.54E-09	1.02E-08

Conclusión de la clase:

Usaré siempre LU, salvo que necesite mucha presición en cuyo caso usaré QR, o Cholesky en el caso de que la matriz sea simétrica definida positiva.

Hay una manera de salvar y usar colesky a pesar de que no sea simetrica positiva. Sin embargo, se pierde tiempo:


>> tic
X=resuelve_chol(A'*A,A'*b);
toc
norm(A*X-b)
Elapsed time is 0.446775 seconds.

ans =

  2.0051e-011

>>



Ejercicio:

Sea :

2X1+3X2-X3=7
3X1+2X2+X3=7.5
-X1+3X2+X3=0.8
X1-X2+X3=1.1

(tiene mas ecuaciones que incognitas).

a)
Hallar X=
X1
X2
X3

solución de cuadrados mínimos.

b)
Comparar [A*X0 y b]

c)
Hallar la norm(AX-b)

d)
Elegir a gusto o al azar otra Xbis, y comparar AX0-b y AXbis-b



a)
>> A=[2 3 -1; 3 2 1; -1 3 1;1 -1 1];
>> b=[7;7.5;0.8;1.1];
>> X=resuelve_chol(A'*A,A'*b);
>> X

X =

    1.9569
    0.9466
   -0.1259

>>

b)
>> A*X

ans =

    6.8793
    7.6379
    0.7569
    0.8845

>> b

b =

    7.0000
    7.5000
    0.8000
    1.1000

>>

c)
>> norm(A*X-b)

ans =

    0.2862


d)
>> Xbis=[2;1;0.125];
>> norm(A*Xbis-b)

ans =

    0.7159

>>