DESCOMPOSICIÓN LU

DESCOMPOSICIÓN LU

Implementar el algoritmo de descomposición LU visto en teoría.
Descomposición LU (I)
función lu_.m:
function [L,U]=lu_(a)
% descomposicion LU sin pivotacion
% a: matriz cuadrada
[n,m]=size(a);
if n~=m error('Matriz no cuadrada.'); end
for j=1:n-1
if a(j,j)==0 error('Pivote nulo.'); end
for i=(j+1):n
factor=a(i,j)/a(j,j);
a(i,j)=factor;
fork=j+1:n
a(i,k)=a(i,k)-factor*a(j,k);
end
end
end
L=tril(a,-1)+eye(n);
U=triu(a);



Calcula U triangular superior, L triangular inferior con 1s en la diagona de forma que A=L*U


function [l,u]=LUdoolitle(a)

n=length(a);
l=a*0; u=a*0;
for k=1:n
u(k,k:n)=a(k,k:n)-l(k,1:k-1)*u(1:k-1,k:n);
l(k,k)=1;
l(k+1:n,k)=(a(k+1:n,k)-l(k+1:n,1:k-1)*u(1:k-1,k))/u(k,k);
endEscribe un programa que permita obtener las matrices correspondientes a la factorización LU, a partir de la
matriz de coeficientes de un sistema de ecuaciones lineales.
La función en Matlab debe admitir como entradas la matriz de coeficientes de un sistema y un vector columna
con los términos independientes, y debe devolver, las matrices L,U y P y un vector columna con las
soluciones delsistema.
function [P, L, U,x] = plu(A,b)
[m, n] = size(A);
if m ~= n
error('La matriz debe ser cuadrada')
end
P = eye(n, n);
L = eye(n, n);
U = zeros(n, n);
tol = sqrt(eps);
sign = 1;

for k = 1:n
if abs(A(k, k)) < tol
for r = k:n
if abs(A(r, k)) >= tol
break
end
if r == n
if nargout == 4
sign = 0;
return
else
disp('A es singular con tolerancia')
error(['No hay pivotaje enlas columnas' int2str(k) '.'])
end
end
end
A([r k], 1:n) = A([k r], 1:n);
if k > 1, L([r k], 1:k-1) = L([k r], 1:k-1); 
end
P([r k], 1:n) = P([k r], 1:n);
sign = -sign;
end
for i = k+1:n
L(i, k) = A(i, k) / A(k, k);
for j = k+1:n
A(i, j) = A(i, j) - L(i, k)*A(k, j);
end
end
for j = k:n
U(k, j) = A(k, j) * (abs(A(k, j)) >= tol);
end
end

Z=P*b;
Y=inv(L)*Z;
x=inv(U)*Y;
L
U P
disp('La solución del sistema es:')
disp(x)
end

4 Factorizacion LU (sin intercambio de filas)
%
% [L, U] = slu(A) usa la eliminacion gausiana para computar una matriz
% triangular inferior L y una matriz triangular superior U tal forma que
% L*U = A. El algoritmo se detendra si la entrada de pivote es muy pequena.
[n, n] = size(A);
for k = 1:n
if abs(A(k, k)) < sqrt(eps)disp([’Pivote muy pequeno encontrado en columna’ int2str(k) ’.’])
end
L(k, k) = 1;
for i = k+1:n
L(i,k) = A(i, k) / A(k, k);
for j = k+1:n
A(i, j) = A(i, j) - L(i, k)*A(k, j);
end
end
for j = k:n
U(k, j) = A(k, j);
end
end






Calcula la descomposición LU con pivoteo
% parcial de una matriz A de orden N x N
%
% Uso
% [B, pinfo, error] = lud(A)
%
% entrada:
% A: La matriz A de ordenN x N a factorizar
%
% salida:
% B: La descomposicion LU de A
% pinfo: informacion del pivote
% error: Bandera de error
% error = 0 Indica que la reduccion se pudo
% realizar, B es triangular superior
% error = 1 Indica que la dimension de A o B
% no es correcta y hubo error
8%
iflag = 0;
% Obtener el tamanio de A y comprobar las dimensiones
[m,n] = size(A);
if ( m ~= n )
iflag =1;
return
end
% Empieza la factorizacion LU con descomposicion
for k = 1:n-1
% Buscar en pivote indice
[amax, l ] = max(abs(A(k:n,k)));
l = l + k - 1;
ipivt(k) = l;
% Intercambiar columnas si es necesario
if( k ~= l )
tmp = A(k,k:n);
A(k,k:n) = A(l,k:n);
A(l,k:n) = tmp;
end
if( amax > 0 ) % Si amax == 0 subcolumna A(k:n,k) es cero
for i = k+1:n
% Computa multiplicadores
A(i,k) =-A(i,k)/A(k,k);
% Realiza la eliminacion por filas
A(i,k+1:n) = A(i,k+1:n) + A(i,k) * A(k,k+1:n);
end
end
end




ECUACIONES TRIDIAGONALES

Si se tiene una matriz tridiagonal de la forma





Entonces en una matriz de 10000 elementos, n=10. Hay E=9702 elementos iguales a cero. Y 208 elementos distintos de cero. Por lo que resulta conveniente almacenar la matriz de una forma...