9na clase de Métodos numéricos

Ejercicio:

A=

2  -1  1

4  3  0

-5  -4  1

X=

X1

X2

X3

B=

1

-1.5

2

>> A=[-2 -1 -1;4 3 0 ; -5 -4 1]

A =

    -2    -1    -1

     4     3     0

    -5    -4     1

>> b

b =

    1.0000

   -1.5000

    2.0000

>> [U,S,V]=svd(A);

>> U

U =

   -0.2530   -0.8658    0.4316

    0.5920    0.2143    0.7769

   -0.7652    0.4521    0.4584

>> S

S =

    8.4296         0         0

         0    1.3909         0

         0         0    0.0853

>> V

V =

    0.7948    0.2361   -0.5590

    0.6038   -0.2154    0.7675

   -0.0608    0.9475    0.3138

>> Z=U'*b

Z =

   -2.6714

   -0.2831

    0.1831

>> S^(-1)

ans =

    0.1186         0         0

         0    0.7190         0

         0         0   11.7247

>> diag(1./diag(S))

ans =

    0.1186         0         0

         0    0.7190         0

         0         0   11.7247

>> W=diag(1./diag(S))*Z

W =

   -0.3169

   -0.2035

    2.1467

>> X=V*W

X =

   -1.5000

    1.5000

    0.5000

>> A*X

ans =

    1.0000

   -1.5000

    2.0000

>> 

Vamos a armarlo en una 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;

Prueba 1:

>> A

A =

    -2    -1    -1

     4     3     0

    -5    -4     1

>> b

b =

    1.0000

   -1.5000

    2.0000

>> X=resuelve_svd(A,b);

>> A*X

ans =

    1.0000

   -1.5000

    2.0000

>>checkea bien

Ejercicio:

a)
Tomar A=randn(10,10)

B=randn(10,1)

Hallar X sol de AX=b usando resuelve_svd y verificar que A*X=b

b)

Tomar A=randn(1500,1500)

>> A=randn(10,10);

>> b=randn(10,1);

>> X=resuelve_svd(A,b)

X =

    4.4404

    0.3800

   -0.6022

   -2.2681

   -0.0346

    1.6231

   -1.0064

   -3.5450

    4.0005

    3.0330

>> A*X

ans =

   -1.1878

   -2.2023

    0.9863

   -0.5186

    0.3274

    0.2341

    0.0215

   -1.0039

   -0.9471

   -0.3744

>> b

b =

   -1.1878

   -2.2023

    0.9863

   -0.5186

    0.3274

    0.2341

    0.0215

   -1.0039

   -0.9471

   -0.3744

>> norm(A*X-b)

ans =

  7.9930e-015

b)

>> A=randn(1500,1500);

>> b=randn(1500,1);

>> X=resuelve_svd(A,b);

>> norm(A*X-b)

ans =

  2.7654e-010

Nunca va a dar 0, por errores de redondeo, pero 10^-10 es suficientemente parecido a 0 para considerarlos iguales. Checkea bien :)


Cuando n es muy grande, 1500 p ejemplo, hay problemas con el tiempo que tarda el programa en calcular la descomposicion de matrices.
Para eso, se usan otras descomposiciones posibles, como por ejemplo la LU:
 

 Resuelve_LU requiere de dos funciones previas:

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


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

luego:

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

solcionado