Superponiendo las matrices \( \mathbf{L} \) y \( \mathbf{U} \) $$ \begin{equation*} \begin{array}{llll} l_{11} = a_{11} & u_{12} = \cfrac{a_{12}}{l_{11}} & u_{13} = \cfrac{a_{13}}{l_{11}} & u_{14} = \cfrac{a_{14}}{l_{11}} \\ l_{21} = a_{21} & l_{22} = a_{22} - l_{21} u_{12} & u_{23} = \cfrac{a_{23} - l_{21} u_{13}}{l_{22}} & u_{24} = \cfrac{a_{24} - l_{21} u_{14}}{l_{22}} \\ l_{31} = a_{31} & l_{32} = a_{32} - l_{31} u_{12} & l_{33} = a_{33} - l_{31} u_{13} - l_{32} u_{23} & u_{34} = \cfrac{a_{34} - l_{31} u_{14} - l_{32} u_{24}}{l_{33}} \\ l_{41} = a_{41} & l_{42} = a_{42} - l_{41} u_{12} & l_{43} = a_{43} - l_{41} u_{13} - l_{42} u_{23} & l_{44} = a_{44} - l_{41} u_{14} - l_{42} u_{24} - l_{43} u_{34} \end{array} \end{equation*} $$
Se usara una matriz para ahorrar memoria $$ \begin{equation*} \begin{array}{llll} a_{11} = a_{11} & a_{12} = \cfrac{a_{12}}{a_{11}} & a_{13} = \cfrac{a_{13}}{a_{11}} & a_{14} = \cfrac{a_{14}}{a_{11}} \\ a_{21} = a_{21} & a_{22} = a_{22} - a_{21} a_{12} & a_{23} = \cfrac{a_{23} - a_{21} a_{13}}{a_{22}} & a_{24} = \cfrac{a_{24} - a_{21} a_{14}}{a_{22}} \\ a_{31} = a_{31} & a_{32} = a_{32} - a_{31} a_{12} & a_{33} = a_{33} - a_{31} a_{13} - a_{32} a_{23} & a_{34} = \cfrac{a_{34} - a_{31} a_{14} - a_{32} a_{24}}{a_{33}} \\ a_{41} = a_{41} & a_{42} = a_{42} - a_{41} a_{12} & a_{43} = a_{43} - a_{41} a_{13} - a_{42} a_{23} & a_{44} = a_{44} - a_{41} a_{14} - a_{42} a_{24} - a_{43} a_{34} \end{array} \end{equation*} $$
function lu_crout(A)
for j=2 to n do
a(1,j) = a(1,j)/a(1,1)
end for
for j=2 to n-1 do
suma = a(j,j)
for k=1 to j-1 do
suma = suma - a(j,k)*a(k,j)
end for
a(j,j) = suma
for i=j+1 to m do
sumav = a(i,j)
sumah = a(j,i)
for k=1 to j-1 do
sumav = sumav - a(i,k)*a(k,j)
sumah = sumah - a(j,k)*a(k,i)
end for
a(i,j) = sumav
a(j,i) = sumah/a(j,j)
end for
end for
suma = a(n,n)
for k=1 to n-1 do
suma = suma - a(n,k)*a(k,n)
end for
a(n,n) = suma
return a
end function
import numpy as np
def lu_crout(A):
a = np.copy(A)
m, n = a.shape
for j in range(1,n):
a[0,j] = a[0,j]/a[0,0]
for j in range(1,n-1):
suma = a[j,j]
for k in range(0,j):
suma = suma - a[j,k]*a[k,j]
a[j,j] = suma
for i in range(j+1,m):
sumav = a[i,j]
sumah = a[j,i]
for k in range(0,j):
sumav = sumav - a[i,k]*a[k,j]
sumah = sumah - a[j,k]*a[k,i]
a[i,j] = sumav
a[j,i] = sumah/a[j,j]
suma = a[n-1,n-1]
for k in range(0,n-1):
suma = suma - a[n-1,k]*a[k,n-1]
a[n-1,n-1] = suma
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_crout(A)
print(B)
[[ 1. 1. 2. 3.]
[ 2. -1. 5. 5.]
[ 3. -4. 13. 1.]
[ -1. 3. -10. -3.]]
L = np.tril(B)
U = np.triu(B)
np.fill_diagonal(U,1)
print(L) #matriz triangular inferior
print(U) #matriz triangular superior
[[ 1. 0. 0. 0.]
[ 2. -1. 0. 0.]
[ 3. -4. 13. 0.]
[ -1. 3. -10. -3.]]
[[ 1. 1. 2. 3.]
[ 0. 1. 5. 5.]
[ 0. 0. 1. 1.]
[ 0. 0. 0. 1.]]
#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_crout(A)
print(B)
[[ 3. -0.03333333 -0.06666667]
[ 0.1 7.00333333 -0.04188482]
[ 0.3 -0.19 10.01204188]]
L = np.tril(B)
U = np.triu(B)
np.fill_diagonal(U,1)
print(L) #matriz triangular inferior
print(U) #matriz triangular superior
[[ 3. 0. 0. ]
[ 0.1 7.00333333 0. ]
[ 0.3 -0.19 10.01204188]]
[[ 1. -0.03333333 -0.06666667]
[ 0. 1. -0.04188482]
[ 0. 0. 1. ]]
#revisando el resultado
np.dot(L,U)
array([[ 3. , -0.1, -0.2],
[ 0.1, 7. , -0.3],
[ 0.3, -0.2, 10. ]])