Te sugerimos que pruebes esta respuesta en un ambiente controlado antes de enviarlo a producción, saludos.
Solución:
de MATLAB lu siempre realiza el pivoteo por defecto. Si, por ejemplo, tenía un coeficiente diagonal que era igual a 0 cuando intentó hacer el algoritmo de descomposición LU convencional, no funcionará, ya que los coeficientes diagonales son necesarios cuando se realiza la eliminación gaussiana para crear la matriz triangular superior. U por lo que obtendría un error de división por cero. Se requiere pivotar para asegurar que la descomposición sea estable.
Sin embargo, si puede garantizar que los coeficientes diagonales de su matriz no sean cero, es muy simple, pero tendrá que escribir esto por su cuenta. Todo lo que tiene que hacer es realizar la eliminación gaussiana en la matriz y reducir la matriz a la forma escalonada reducida. La matriz de forma escalonada reducida resultante es U mientras que los coeficientes requeridos para eliminar la parte triangular inferior de L en la eliminación gaussiana se colocaría en la mitad triangular inferior para hacer U.
Algo como esto podría funcionar, suponiendo que su matriz esté almacenada en A. Recuerde que estoy asumiendo una matriz cuadrada aquí. La implementación del algoritmo de descomposición LU no pivotante se coloca en un archivo de función de MATLAB llamado lu_nopivot:
function [L, U] = lu_nopivot(A)
n = size(A, 1); % Obtain number of rows (should equal number of columns)
L = eye(n); % Start L off as identity and populate the lower triangular half slowly
for k = 1 : n
% For each row k, access columns from k+1 to the end and divide by
% the diagonal coefficient at A(k ,k)
L(k + 1 : n, k) = A(k + 1 : n, k) / A(k, k);
% For each row k+1 to the end, perform Gaussian elimination
% In the end, A will contain U
for l = k + 1 : n
A(l, :) = A(l, :) - L(l, k) * A(k, :);
end
end
U = A;
end
Como ejemplo, supongamos que tenemos la siguiente matriz de 3 x 3:
>> rng(123)
>> A = randi(10, 3, 3)
A =
7 6 10
3 8 7
3 5 5
Ejecutar el algoritmo nos da:
>> [L,U] = lu_nopivot(A)
L =
1.0000 0 0
0.4286 1.0000 0
0.4286 0.4474 1.0000
U =
7.0000 6.0000 10.0000
0 5.4286 2.7143
0 0 -0.5000
multiplicando L y U juntos da:
>> L*U
ans =
7 6 10
3 8 7
3 5 5
… que es la matriz original A.
Podrías usar este truco (aunque como ya se mencionó, podrías perder la estabilidad numérica):
[L, U] = lu(sparse(A), 0)
Es posible que desee considerar realizar una descomposición LDU en lugar de LU no pivotada. Vea, LU sin pivotar es numéricamente inestable, incluso para matrices que son de rango completo e invertibles. El algoritmo simple proporcionado anteriormente muestra por qué: hay una división por cada elemento diagonal de la matriz involucrada. Por lo tanto, si hay un cero en cualquier parte de la diagonal, la descomposición falla, aunque la matriz aún no sea singular.
Wikipedia habla un poco sobre la descomposición de LDU aquí:
https://en.wikipedia.org/wiki/LU_decomposition#LDU_decomposition
sin citar un algoritmo. Cita el siguiente libro de texto como prueba de existencia:
Cuerno, Roger A.; Johnson, Charles R. (1985), Matrix Analysis, Cambridge University Press, ISBN 978-0-521-38632-6. Consulte la Sección 3.5.
Se garantiza que LDU existe (al menos para una matriz invertible), es numéricamente estable y también es único (siempre que tanto L como U estén restringidos a tener elementos unitarios en la diagonal).
Entonces, si por alguna razón “D” se interpone en su camino, puede absorber la matriz diagonal D en L (L:=LD) o U (U:=DU), o divídalo simétricamente entre L y U (como L:=L*sqrt(D) y U:=sqrt(D)*U), o como quiera hacerlo. Hay un número infinito de formas de dividir LDU en LU, y es por eso que la descomposición de LU no es única.
Sección de Reseñas y Valoraciones
Puedes proteger nuestra tarea fijando un comentario o puntuándolo te damos las gracias.