Algoritmo de Crout revisitado

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*}

Fórmula matemática

\begin{align*} j &= 2, \dots, n \\ & \quad a_{1j} = \frac{a_{1j}}{a_{11}} \\ j &= 2, \dots, n-1 \\ & \quad a_{jj} = a_{jj} - \sum_{k=1}^{j-1} a_{jk} a_{kj} \\ & \quad i = j+1, m \\ & \quad \quad a_{ij} = a_{ij} - \sum_{k=1}^{j-1} a_{ik} a_{kj} \\ & \quad \quad a_{ji} = \frac{a_{ji} - \sum_{k=1}^{j-1} a_{jk} a_{ki}}{a_{jj}} \\ a_{nn} &= a_{nn} - \sum_{k=1}^{n-1} a_{nk} a_{kn} \end{align*}

Seudocódigo

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

Implementación



In [1]:

    
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



In [2]:

    
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.]]



In [3]:

    
#resultado
B = lu_crout(A)
print(B)









    



[[  1.   1.   2.   3.]
 [  2.  -1.   5.   5.]
 [  3.  -4.  13.   1.]
 [ -1.   3. -10.  -3.]]



In [4]:

    
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.]]



In [5]:

    
#revisando el resultado
np.dot(L,U)









    Out[5]:





array([[ 1.,  1.,  2.,  3.],
       [ 2.,  1., -1.,  1.],
       [ 3., -1., -1.,  2.],
       [-1.,  2.,  3., -1.]])



In [6]:

    
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. ]]



In [7]:

    
#resultado
B = lu_crout(A)
print(B)









    



[[  3.          -0.03333333  -0.06666667]
 [  0.1          7.00333333  -0.04188482]
 [  0.3         -0.19        10.01204188]]



In [8]:

    
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.        ]]



In [9]:

    
#revisando el resultado
np.dot(L,U)









    Out[9]:





array([[  3. ,  -0.1,  -0.2],
       [  0.1,   7. ,  -0.3],
       [  0.3,  -0.2,  10. ]])