Superponiendo las matrices \( \mathbf{L} \) y \( \mathbf{U} \) $$ \begin{equation*} \begin{array}{llll} u_{11} = a_{11} & u_{12} = a_{12} & u_{13} = a_{13} & u_{14} = a_{14} \\ l_{21} = \frac{a_{21}}{a_{11}} & u_{22} = a_{22}' & u_{23} = a_{23}' & u_{24} = a_{24}' \\ l_{31} = \frac{a_{31}}{a_{11}} & l_{32} = \frac{a_{32}'}{a_{22}'} & u_{33} = a_{33}'' & u_{34} = a_{34}'' \\ l_{41} = \frac{a_{41}}{a_{11}} & l_{42} = \frac{a_{42}'}{a_{22}'} & l_{43} = \frac{a_{43}''}{a_{33}''} & u_{44} = a_{44}''' \end{array} \end{equation*} $$
Se usara una matriz para ahorrar memoria $$ \begin{equation*} \begin{array}{llll} a_{11} = a_{11} & a_{12} = a_{12} & a_{13} = a_{13} & a_{14} = a_{14} \\ a_{21} = \frac{a_{21}}{a_{11}} & a_{22} = a_{22}' & a_{23} = a_{23}' & a_{24} = a_{24}' \\ a_{31} = \frac{a_{31}}{a_{11}} & a_{32} = \frac{a_{32}'}{a_{22}'} & a_{33} = a_{33}'' & a_{34} = a_{34}'' \\ a_{41} = \frac{a_{41}}{a_{11}} & a_{42} = \frac{a_{42}'}{a_{22}'} & a_{43} = \frac{a_{43}''}{a_{33}''} & a_{44} = a_{44}''' \end{array} \end{equation*} $$
function lu_doolittle(a)
m, n = tamaño(a)
for k=1 to m-1 do
for i=1+k to m do
for j=1+k to n do
a(i,j) = a(i,j) - a(i,k)*a(k,j)/a(k,k)
end for
a(i,k) = a(i,k)/a(k,k)
end for
end for
end function
otra alternativa para reducir tiempo de cálculo
function lu_doolittle(a)
m, n = tamaño(a)
for k=1 to m-1 do
for i=1+k to m do
factor = a(i,k)/a(k,k)
for j=1+k to n do
a(i,j) = a(i,j) - factor*a(k,j)
end for
a(i,k) = factor
end for
end for
end function
import numpy as np
def lu_doolittle(A):
a = np.copy(A)
m, n = a.shape
for k in range(0,m-1):
for i in range(1+k,m):
factor = a[i,k]/a[k,k]
for j in range(1+k,n):
a[i,j] = a[i,j] - factor*a[k,j]
a[i,k] = factor
return a
A = np.array([[1,1,2,3],
[2,1,-1,1],
[3,-1,-1,2],
[-1,2,3,-1]],float)
print(A)
[[ 1. 1. 2. 3.]
[ 2. 1. -1. 1.]
[ 3. -1. -1. 2.]
[-1. 2. 3. -1.]]
#resultado
B = lu_doolittle(A)
print(B)
[[ 1. 1. 2. 3. ]
[ 2. -1. -5. -5. ]
[ 3. 4. 13. 13. ]
[ -1. -3. -0.76923077 -3. ]]
L = np.tril(B)
np.fill_diagonal(L,1)
U = np.triu(B)
print(L) #matriz triangular inferior
print(U) #matriz triangular superior
[[ 1. 0. 0. 0. ]
[ 2. 1. 0. 0. ]
[ 3. 4. 1. 0. ]
[-1. -3. -0.76923077 1. ]]
[[ 1. 1. 2. 3.]
[ 0. -1. -5. -5.]
[ 0. 0. 13. 13.]
[ 0. 0. 0. -3.]]
#revisando el resultado
np.dot(L,U)
array([[ 1., 1., 2., 3.],
[ 2., 1., -1., 1.],
[ 3., -1., -1., 2.],
[-1., 2., 3., -1.]])
A = np.array([[3,-0.1,-0.2],
[0.1,7,-0.3],
[0.3,-0.2,10]],float)
print(A)
[[ 3. -0.1 -0.2]
[ 0.1 7. -0.3]
[ 0.3 -0.2 10. ]]
#resultado
B = lu_doolittle(A)
print(B)
[[ 3. -0.1 -0.2 ]
[ 0.03333333 7.00333333 -0.29333333]
[ 0.1 -0.02712994 10.01204188]]
L = np.tril(B)
np.fill_diagonal(L,1)
U = np.triu(B)
print(L) #matriz triangular inferior
print(U) #matriz triangular superior
[[ 1. 0. 0. ]
[ 0.03333333 1. 0. ]
[ 0.1 -0.02712994 1. ]]
[[ 3. -0.1 -0.2 ]
[ 0. 7.00333333 -0.29333333]
[ 0. 0. 10.01204188]]
#revisando el resultado
np.dot(L,U)
array([[ 3. , -0.1, -0.2],
[ 0.1, 7. , -0.3],
[ 0.3, -0.2, 10. ]])