Algoritmo de Doolittle revisitado

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

Fórmula matemática

\begin{align*} k &= 1, \dots, m - 1 \\ & \quad i = 1 + k, \dots, m \\ & \quad \quad j = 1 + k, \dots, n \\ & \quad \quad \quad a_{ij} = a_{ij} - \frac{a_{ik}}{a_{kk}} a_{kj} \\ & \quad \quad a_{ik} = \frac{a_{ik}}{a_{kk}} \end{align*}

Seudocódigo

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

Implementación descomposición de Doolittle



In [1]:

    
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



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_doolittle(A)
print(B)









    



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



In [4]:

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



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_doolittle(A)
print(B)









    



[[  3.          -0.1         -0.2       ]
 [  0.03333333   7.00333333  -0.29333333]
 [  0.1         -0.02712994  10.01204188]]



In [8]:

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



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