2.1.2 Sistemas Lineares — Métodos Diretos — Decomposição LU
1. Introdução </p>
Um sistema linear do tipo
$$ a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n = b_1 $$
$$ a_{21}x_1 + a_{22}x_2 + \cdots + a_{2n}x_n = b_2 $$
$$ 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 matriz formada pelos coeficientes , é o vetor formado pelas variáveis a serem determinadas e é 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 , que consomem maior tempo, das manipulações envolvidas com o vetor dos termos independentes, . 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 </p>
A equação
$$ AX = B $$
pode ser reescrita como
$$ 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:
$$ 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 em duas matrizes, onde a primeira está na forma triangular inferior (Low), enquanto a outra está na forma triangula superior (Upper). Sendo assim, para e , comparadascom a Equação 5, temos que
$$ L[UX – D] = AX – B = 0 $$
Agora, isolamos os termos dependentes de , temos que
$$ LU = A $$
e
$$ LD = B $$
desta forma, isolamos a dependência dos termos independentes da matriz dos coeficientes . Desta forma, também livramos as operações efetuadas sobre (agora ) de serem feitas em , diminuindo a demanda de recursos para resolução deste sistema.
2.1. Eliminação de Gauss como uma Decomposição LU </p>
Como dito na seção anterior, a matriz 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 é 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 . 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 . 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 . Por fim, multiplicamos a segunda linha por
$$\phi_{32} = \dfrac{a’_{31}{a’_{22} $$
e subtraímos o resultado da terceira linha para eliminarmos . Como a ideia neste processo é gerar zeros em , e . Assim, após a eliminação, 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 de :
$$ 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 e de forma não-singular da matriz quadrada , sendo a fatoração de nos produtos da matriz (Lower) e (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):
$$ l_{i1} = a_{i1} $$
para .
$$ u_{1k} = \dfrac{a_{1k}{l_{11} $$
para . Para , temos:
$$ l_{ij} = a_{ij} – \sum^{j-1}_{k=1} l_{ik}u_{kj} $$
onde e também
$$ u_{ji} = \dfrac{a_{ji} – \sum^{j-1}_{k=1} l_{jk}u_{ki}{l_{ii} $$
onde , e
$$ l_{nn} = a_{nn} – \sum^{n-1}_{k=1} l_{nk}u_{kn} $$
3. Implementação </p>
Uma implementação eficiente para decomposição 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. </p>
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/.