2.1.2 Sistemas Lineares — Métodos Diretos — Decomposição LU

1. Introdução

Um sistema linear do tipo
(eq1)


\( a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n = b_1 \)

(eq2)

\( a_{21}x_1 + a_{22}x_2 + \cdots + a_{2n}x_n = b_2 \)

\(\vdots \)

(eq3)

\( a_{n1}x_1 + a_{n2}x_2 + \cdots + a_{nn}x_n = b_n \)

pode ser representado na seguinte forma matricial:

\(AX = B \)

onde \(A\) é a matriz formada pelos coeficientes \(a_{11}, a_{12}, \ldots, a_{nn} \), \(X \) é o vetor formado pelas variáveis a serem determinadas \(x_1, x_2, \ldots, x_n \) e \(B\) é o vetor dos termos independentes.

Para determinação das incógnitas, o método da eliminação de Gauss desenvolve duas fases: a primeira é a eliminação progressiva, onde reduz o número de variáveis ao longo da execução para, então, aplicar a segunda fase, chamada de substituição regressiva, onde utiliza o resultado da primeira para determinar a solução geral do sistema.

Dois dois passos descritos, o primeiro é o que consome mais tempo de cálculo, uma vez que é nesta fase que consiste o maior número de operações aritméticas e de trocas de dados. Por isso, encontrar um método que minimize esta fase crítica, implica em aumentar o desempenho para realizar a tarefa de resolução de sistemas lineares.

Os métodos de decomposição LU consistem em separar a fase de eliminação da matriz dos coeficientes \(A \) , que consomem maior tempo, das manipulações envolvidas com o vetor dos termos independentes, \(B \) . Portanto, devemos deixar claro que, ao contrário da eliminação de Gauss, uma decomposição de LU é uma estratégia de melhoria na resolução de sistemas lineares. Sendo assim, não existe “o método” de decomposição LU, mas sim algumas abordagens a serem tomadas que permitem decompor o sistema. Uma implicação interessante disso é que a própria eliminação de Gauss pode ser descrita como uma decomposição LU. Ralson e Rabinowitz utilizam essa abordagem ao tratar os métodos diretos em seu livro A First Course in Numerical Analysis[rabinowitz], desenvolvendo os principais algoritmos de decomposição — Crout, Doolittle e Cholesky — como técnicas de resolução direta de sistemas lineares.

2. Desenvolvimento do método da decomposição LU

A equação
(eq4)


\( AX = B \)

pode ser reescrita como
(eq5)

\( AX – B = 0 \)

Aplicando a eliminação de Gauss, sabemos que este sistema pode ser reescrito como uma matriz triangular superior na forma:


\(\left[ \begin{array}{ccc} u_{11} & u_{12} & u_{13} \\ 0 & u_{22} & u_{23} \\ 0 & 0 & u_{33} \end{array} \right]\) \(\left[ \begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array} \right] = \) \(\left[ \begin{array}{c} d_1 \\ d_2 \\ d_3 \\ \end{array} \right]\)

Esta notação matricial também pode ser usada da seguinte maneira:
(eq6)

\( UX – D = 0 \)

Agora, vamos supor que existe uma matriz composta apenas pela formula triangular inferior, tal que


\(L = \left[ \begin{array}{ccc} 1 & 0 & 0 \\ l_{21} & 1 & 0 \\ l_{31} & l_{32} & 1 \end{array} \right]\)

O processo de decomposição LU consiste justamente de em decompor a matriz dos coeficientes \(A \) em duas matrizes, onde a primeira está na forma triangular inferior (Low), enquanto a outra está na forma triangula superior (Upper). Sendo assim, para \(L \) e \(U \) , comparadascom a Equação 5, temos que
(eq7)


\( L[UX – D] = AX – B = 0 \)

Agora, isolamos os termos dependentes de \(X \) , temos que
(eq8)


\( LU = A \)

e
(eq9)

\( LD = B \)

desta forma, isolamos a dependência dos termos independentes \(B \) da matriz dos coeficientes \(A \) . Desta forma, também livramos as operações efetuadas sobre \(A \) (agora \(LU \) ) de serem feitas em \(B \) , diminuindo a demanda de recursos para resolução deste sistema.

2.1. Eliminação de Gauss como uma Decomposição LU

Como dito na seção anterior, a matriz \(U \) possui uma forma muito semelhante àquela gerada pelo processo de eliminação progressiva, estágio da eliminação de Gauss. Neste estágio, a matriz dos coeficientes \(A \) é reduzida para a forma


\(\left[ \begin{array}{ccc} a_{11} & a_{12} & a_{13} \\ 0 & a’_{22} & a’_{23} \\ 0 & 0 & a”_{33} \end{array} \right]\)

que é a forma triangular superior procurada. Desta forma, os termos retirados são aqueles que formam a matriz \(L \) . Podemos ver isso ao agruparmos em uma mesma matriz os elementos obtidos da eliminação progressiva, e os retirados pelo processo. Isto é, para o sistema original

\(\left[ \begin{array}{ccc} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{array} \right] \left[ \begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array} \right] = \left[ \begin{array}{c} b_1 \\ b_2 \\ b_3 \\ \end{array} \right] \)

a primeira etapa na eliminação de Gauss é multiplicar a primeira linha pelo fator de normalização

\(\phi_{21} = \dfrac{a_{21}}{a_{11}} \)

e subtrair este resultado da próxima linha, a fim de eliminar o termo \(a_{21} \) . Da mesma forma, o próximo passo se dá em multiplicar a primeira linha por

\(\phi_{31} = \dfrac{a_{31}}{a_{11}} \)

e o resultado é subtraído da terceira linha para eliminar \(a_{31} \) . Por fim, multiplicamos a segunda linha por

\(\phi_{32} = \dfrac{a’_{31}}{a’_{22}} \)

e subtraímos o resultado da terceira linha para eliminarmos \(a’_{32} \). Como a ideia neste processo é gerar zeros em \(a_{21} \) , \(a_{31} \) e \(a_{32} \) . Assim, após a eliminação, \(A \) pode ser escrita na forma:

\(\left[ \begin{array}{ccc} a_{11} & a_{12} & a_{13} \\ \phi_{21} & a’_{22} & a’_{23} \\ \phi_{31} & \phi_{32} & a”_{33} \end{array} \right]\)

Essa matriz, de fato, representa uma decomposição \(LU \) de \(A \) :
(eq10)

\( A = LU \)

pois

\(U = \left[ \begin{array}{ccc} a_{11} & a_{12} & a_{13} \\ 0 & a’_{22} & a’_{23} \\ 0 & 0 & a”_{33} \end{array} \right]\)

e

\(L = \left[ \begin{array}{ccc} 1 & 0 & 0 \\ \phi_{21} & 1 & 0 \\ \phi_{31} & \phi_{32} & 1 \end{array} \right]\)

2.2. O algoritmo de Crout

O algoritmo de Crout decompõe as matrizes \(U \) e \(L \) de forma não-singular da matriz quadrada \(A \) , sendo a fatoração de \(A \) nos produtos da matriz \(L \) (Lower) e \(U \) (Upper), cuja diagonal principal é formada por valores unitários. O processo de decomposição é feito a partir das seguintes fórmulas (Obtidas de Ralson e Rabinowitz A First Course in Numerical Analysis[rabinowitz], páginas 422 e 423):
(eq10)


\( l_{i1} = a_{i1} \)

para \(i = 1,2,\ldots, n \) .
(eq11)

\( u_{1k} = \dfrac{a_{1k}}{l_{11}} \)

para \(k = 2, 3, \ldots, n \). Para \(j = 2, 3, \ldots, n-1 \) , temos:
(eq12)

\( l_{ij} = a_{ij} – \sum^{j-1}_{k=1} l_{ik}u_{kj} \)

onde \(i = j, j+1, j+2, \ldots, n \) e também
(eq13)

\( u_{ji} = \dfrac{a_{ji} – \sum^{j-1}_{k=1} l_{jk}u_{ki}}{l_{ii}} \)

onde \(i = j + 1, j + 2, \ldots, n \) , e
(eq14)

\( l_{nn} = a_{nn} – \sum^{n-1}_{k=1} l_{nk}u_{kn} \)

3. Implementação

Uma implementação eficiente para decomposição \(LU \) segue abaixo:

def LU(A, b):
    """
    Return a list with the solution of linear system.
 
    LU(A, b)
 
    INPUT:
    * A: a matrix of coefficients
    * b: a list, relative a vector of independent terms
 
    return: a list with all unknown variables of system.
 
    Author: Pedro Garcia [sawp @sawp.com.br]
    see: http://www.sawp.com.br
 
    License: Creative Commons
             http://creativecommons.org/licenses/by-nc-nd/2.5/br/
 
    Mar 2010
    """
 
    def pivot(A, k):
        n = len(b)
        p = k
        bigger = abs(A[k][k])
        for i in xrange(k+1, n):
            sentinel = abs(A[i][k])
            if (sentinel > bigger):
                bigger = sentinel
                p = i
        if p != k:
            for i in xrange(k, n):
                (A[p][i], A[k][i]) = (A[k][i], A[p][i])
        return A
 
    def crout_decomposition(M):
        """ Perform the Crout-based algorithm for decomposition the matriz
        A in two L and U using the Croud's Method.
        """
        n = len(M)
        L = [[0.0 for i in xrange(n)] for j in xrange(n)]
        U = [[0.0 for i in xrange(n)] for j in xrange(n)]
 
        for j in xrange(n):
            U[0][j] = M[0][j]
        for i in xrange(n):
            L[i][0] = M[i][0]/U[0][0]
 
        for i in xrange(n):
            # Calculate L
            for j in xrange(i+1):
                suma = 0.0
                for k in xrange(j):
                    suma += L[i][k] * U[k][j]
                L[i][j] = M[i][j] - suma
 
            # Calculate U
            for j in xrange(i, n):
                suma = 0.0
                for k in xrange(i):
                    suma += L[i][k] * U[k][j]
                U[i][j] = (M[i][j] - suma)/L[i][i]
        return (L, U)
 
    def solveLU(A, b):
        n = len(b)
        y = [0.0 for i in xrange(n)]
        x = [0.0 for i in xrange(n)]
        (L, U) = crout_decomposition(A)
 
        # Step 1: Solve the L system: L * y = b
        y[0] = b[0]/L[0][0]
        for i in xrange(1, n):
            suma = 0.0
            for j in xrange(i):
                suma += L[i][j] * y[j]
            y[i] = (b[i] - suma)/L[i][i]
 
        # Step 2: Solve the U system: U * x = y (regressive replacement)
        x[n-1] = y[n-1]/U[n-1][n-1]
        for i in xrange(n-1, -1, -1):
            suma = y[i]
            for j in xrange(i+1, n):
                suma = suma - U[i][j] * x[j]
            x[i] = suma/U[i][i]
        return x
 
    # Function body
    return solveLU(A, b)

Esta função está disponível em http://www.sawp.com.br/code/linear/LU.py, assim como um exemplo de como utilizá-la.

4. Copyright

Este documento é disponível sob a licença Creative Commons. As regras dos direitos de cópia deste conteúdo estão acessíveis em http://creativecommons.org/licenses/by-nc-nd/2.5/br/.

References

[rabinowitz] Ralston, Anthony and Rabinowitz, Philip, A First Course in Numerical Analysis (2nd ed.) McGraw-Hill and Dover, (2001), 425.