Sistemas Lineares na Prática Computacional: Métodos Diretos e Fatoração LU
Contents
import matplotlib.pyplot as plt
plt.style.use('../styles/gcpeixoto-book.mplstyle')
8. Sistemas Lineares na Prática Computacional: Métodos Diretos e Fatoração LU#
Objetivos de aprendizagem
- Compreender os fundamentos dos métodos diretos para a resolução de sistemas lineares e sua relação com a fatoração de matrizes.;
- Aplicar a fatoração LU como estratégia eficiente para resolver sistemas lineares;
- Implementar algoritmos de fatoração LU e resolver sistemas computacionalmente;
- Analisar o condicionamento de matrizes e compreender seu impacto na estabilidade numérica das soluções computacionais;
""
8.1. Métodos diretos#
A resolução de sistemas lineares é uma das tarefas mais fundamentais da engenharia computacional, estando no cerne de aplicações que vão desde a análise estrutural até a simulação de escoamentos e processos térmicos. Um sistema linear é tipicamente representado na forma matricial \(\mathbf{A} \mathbf{x} = \mathbf{b}\), onde \(\mathbf{A}\) é uma matriz de coeficientes conhecida, \(\mathbf{b}\) é um vetor de termos independentes conhecido, e \(\mathbf{x}\) é o vetor de incógnitas a ser determinado. Para resolver esse tipo de problema, pode-se recorrer a métodos iterativos ou diretos. Os métodos diretos visam encontrar a solução exata (dentro dos limites da precisão numérica da máquina) por meio de um número finito de operações. Esses métodos são preferidos quando a matriz \(\mathbf{A}\) é de tamanho moderado e bem condicionada, pois permitem resolver o sistema com robustez e sem necessidade de aproximações sucessivas.
8.1.1. Fatoração de matrizes#
O conceito de fatoração em álgebra linear tem forte analogia com a fatoração algébrica aprendida ainda no ensino médio, como quando escrevemos \(x^2 - 9\) como \((x - 3)(x + 3)\), ou \(x^2 + 5x + 6\) como \((x + 2)(x + 3)\). Nessas expressões, decompor um polinômio em fatores mais simples permite resolver equações com mais facilidade. Da mesma forma, na álgebra matricial, a fatoração de uma matriz busca expressar uma matriz \(\mathbf{A}\) como o produto de duas ou mais matrizes com estrutura conhecida e propriedades computacionalmente vantajosas. Por exemplo, dada uma matriz
podemos reescrevê-la como o produto
onde \(\mathbf{L}\) é uma matriz triangular inferior com 1s na diagonal, e \(\mathbf{U}\) é uma matriz triangular superior. Essa fatoração, em particular, é chamada de fatoração LU, que generaliza a ideia da eliminação Gaussiana, de transformar a matriz densa associada ao sistema em uma matriz triangular superior, para facilitar a resolução. Essa técnica é especialmente poderosa e serve como base para variações mais sofisticadas, como a fatoração com pivotamento e a decomposição de Cholesky para matrizes simétricas positivas definidas. O quadro abaixo resume os principais esquemas de fatoração conhecidos e as condições de funcionamento.
Fatoração |
Forma da Decomposição |
Condições / Observações |
|---|---|---|
LU (genérica) |
\(\mathbf{A}\) = \(\mathbf{LU}\) |
A quadrada e não singular. Pode requerer pivotamento para estabilidade. |
LU de Doolittle |
\(\mathbf{A}\) = \(\mathbf{LU}\), com \(\textrm{diag}(L) = 1\) |
A quadrada. \(\mathbf{L}\) com 1s na diagonal; facilita implementação. |
LU de Crout |
\(\mathbf{A}\) = \(\mathbf{LU}\), com \(\textrm{diag}(U) = 1\) |
A quadrada. \(\mathbf{U}\) com 1s na diagonal; alternativa ao Doolittle. |
LU com pivotamento |
\(\mathbf{P}\mathbf{A} = \mathbf{LU}\) ou \(\mathbf{PAQ}\) = \(\mathbf{LU}\) |
Usa permutação \(\mathbf{P}\) e \(\mathbf{Q}\) para melhorar estabilidade numérica. |
Cholesky |
\(\mathbf{A}\) = \(\mathbf{LL}^T\) |
A simétrica, definida positiva. Muito eficiente; reduz pela metade o custo. |
LDL\(^T\) |
\(\mathbf{A}\) = \(\mathbf{LDL}^T\) |
A simétrica. \(\mathbf{D}\) é diagonal; útil em substituição ao Cholesky quando \(\mathbf{A}\) não é positiva definida. |
QR |
\(\mathbf{A} = QR\) |
\(\mathbf{Q}\) ortogonal, \(\mathbf{R}\) triangular superior. Muito usada em problemas de mínimos quadrados. |
Thomas |
\(\mathbf{A}\) = \(\mathbf{LU}\) otimizada |
Para matrizes tridiagonais. Complexidade linear \(\mathcal{O}(n)\). |
LU blocada |
\(\mathbf{A}\) = \(\mathbf{LU}\), por blocos |
Usado para matrizes muito grandes ou estruturadas em blocos. |
SVD |
\(\mathbf{A}\) = \(\mathbf{U\Sigma V}^T\) |
Decomposição de valores singulares. Mais geral, usada em problemas mal-condicionados. |
Neste capítulo, exploraremos o processo de fatoração LU de Doolittle cujo algoritmo busca resolver o sistema \(\mathbf{A}\mathbf{x} = \mathbf{b}\), em duas etapas principais. Ambas as etapas envolvem apenas substituições diretas — progressiva e regressiva — que são muito mais eficientes do ponto de vista computacional do que operar diretamente com \(\mathbf{A}\). Assim, a fatoração transforma o problema original em uma sequência de tarefas mais simples e rápidas de resolver, além de permitir reaproveitamento em casos com múltiplos vetores \(\mathbf{b}\).
8.2. Fatoração LU por substituição de recorrência#
Consideraremos que as matrizes triangulares inferiores \(\mathbf{L}\) sempre terão a sua diagonal principal formada por entradas iguais a 1.
onde
Multiplicando \(\mathbf{L}\mathbf{U}\) e definindo a resposta igual a \(\mathbf{A}\) temos:
Agora, através de substituição de recorrência, facilmente encontramos \(\mathbf{L}\) e \(\mathbf{U}\).
8.3. Algoritmo da fatoração LU#
Uma vez decomposta a matriz \(\mathbf{A}\) no produto \(\mathbf{L}\mathbf{U}\), podemos obter a solução \(\mathbf{x}\) de forma direta. O procedimento, equivalente à substituição de recorrência mencionada anteriormente, pode ser resumido como segue: dada \(\mathbf{A}\), encontre \(\mathbf{L}\) e \(\mathbf{U}\) tal que \(\mathbf{A} = \mathbf{L}\mathbf{U}\), ou seja, \((\mathbf{L}\mathbf{U})\mathbf{x} = \mathbf{b}\). Em outras palavras:
Defina \(\mathbf{y} = \mathbf{U}\mathbf{x}\).
Resolva o sistema triangular \(\mathbf{L}\mathbf{y} = \mathbf{b}\) para \(\mathbf{y}\).
Finalmente, resolva o sistema triangular \(\mathbf{U}\mathbf{x} = \mathbf{y}\) para \(\mathbf{x}\).
O benefício desta abordagem é a resolução de somente sistemas triangulares. Por outro lado, o custo empregado é termos de resolver dois sistemas.
8.3.1. Exemplo#
Encontre a solução \(\mathbf{x} = \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix}\) do sistema
As matrizes \(\mathbf{L}\) e \(\mathbf{U}\) já foram obtidas anteriormente.
A próxima etapa é resolver \(\mathbf{L}\mathbf{y} = \mathbf{b}\), para o vetor \(\mathbf{y} = \begin{bmatrix} y_1 \\ y_2 \\ y_3 \end{bmatrix}\).
Este sistema pode ser resolvido por substituição direta, obtendo \(\mathbf{y} = \begin{bmatrix} 3 \\ 4 \\ -6 \end{bmatrix}\).
Agora que encontramos \(\mathbf{y}\), concluímos o procedimento resolvendo \(\mathbf{U}\mathbf{x} = \mathbf{y}\) para \(\mathbf{x}\). Ou seja, resolvemos:
Realizando a substituição regressiva (baixo para cima; direita para esquerda), obtemos a solução do problema.
8.3.2. Implementação computacional#
Em implementações computacionais, é comum armazenar as matrizes \(\mathbf{L}\) e \(\mathbf{U}\) sobre a própria estrutura da matriz original, economizando memória: os elementos abaixo da diagonal representam \(\mathbf{L}\) (sem armazenar os 1s da diagonal), e os elementos da diagonal para cima representam \(\mathbf{U}\). Além disso, o custo computacional da fatoração LU é da ordem de \(\mathcal{O}(n^3)\) para matrizes densas, com \(n\) sendo o número de variáveis, mas isso pode ser reduzido consideravelmente em matrizes esparsas ou com estrutura especial (como tridiagonais).
Abaixo, temos uma implementação de uma fatoração LU sem pivoteamento.
import numpy as np
def lu_nopivot(A):
'''
Realiza fatoração LU para a matriz A
entrada:
A - matriz: array (nxn)
saida:
L,U - matriz triangular inferior, superior : array (nxn)
'''
n = np.shape(A)[0] # dimensao
L = np.eye(n) # auxiliar
for k in np.arange(n):
L[k+1:n,k] = A[k+1:n,k]/A[k,k]
for l in np.arange(k+1,n):
A[l,:] = A[l,:] - np.dot(L[l,k],A[k,:])
U = A
return (L,U)
Vamos testar a implementação anterior para uma matriz de quarta ordem:
A = np.array([[ 4., -2., -3., 6.],[ 1., 4., 2., 3.],[ 2., -3., 3., -2.],[ 1., 5., 3., 4.]])
L,U = lu_nopivot(A)
L, U
(array([[ 1. , 0. , 0. , 0. ],
[ 0.25 , 1. , 0. , 0. ],
[ 0.5 , -0.44444444, 1. , 0. ],
[ 0.25 , 1.22222222, 0.06796117, 1. ]]),
array([[ 4. , -2. , -3. , 6. ],
[ 0. , 4.5 , 2.75 , 1.5 ],
[ 0. , 0. , 5.72222222, -4.33333333],
[ 0. , 0. , 0. , 0.96116505]]))
8.4. Resolução de sistemas com scipy.linalg#
Métodos adequados para a resolução de sistemas lineares e realizar operações no escopo da Álgebra Linear são encontrados no submódulo linalg do scipy. Importamos essas funcionalidades com:
from scipy import linalg
Vamos calcular a solução do sistema linear \({\bf A}{\bf x} = {\bf b}\) com
Vamos importar os módulos e escrever a matriz \({\bf A}\).
import numpy as np
from scipy import linalg
A = np.array([[4,-2,-3,6],[-6,7,6.5,-6],[1,7.5,6.25,5.5],[-12,22,15.5,-1]])
print(A)
[[ 4. -2. -3. 6. ]
[ -6. 7. 6.5 -6. ]
[ 1. 7.5 6.25 5.5 ]
[-12. 22. 15.5 -1. ]]
Agora, vamos escrever o vetor \({\bf b}\).
b = np.array([12,-6.5,16,17])
print(b)
[12. -6.5 16. 17. ]
Podemos checar as dimensões com
# dimensões de A
A.shape
(4, 4)
# dimensão de b
b.shape
(4,)
A solução do sistema pode ser obtida através do método linalg.solve.
x = linalg.solve(A,b)
print(x)
[ 2. 4. -3. 0.5]
Matematicamente, a solução do sistema anterior é dada por \({\bf x} = {\bf A}^{-1}{\bf b}\). Podemos até invocar a matriz inversa aqui com linalg.inv(A).dot(b) e a solução é a mesma que no caso anterior.
x2 = linalg.inv(A).dot(b)
print(x2)
[ 2. 4. -3. 0.5]
Por outro lado, este método é ineficiente. Computacionalmente, a inversão de matrizes não é aconselhável.
8.4.1. Verificação da solução#
Podemos usar o fato de que \({\bf A}{\bf A}^{-1}{\bf b} - {\bf b} = {\bf 0}\).
x3 = A.dot(linalg.inv(A).dot(b)) - b
print(x3)
[-2.48689958e-14 4.26325641e-14 3.55271368e-15 8.52651283e-14]
Note que o vetor é próximo do vetor nulo, mas não identicamente nulo.
Podemos também computar a norma do resíduo (erro): \(||{\bf r}|| = ||{\bf b} - {\bf A}{\bf x}|| = \langle {\bf b} - {\bf A}{\bf x}, {\bf b} - {\bf A}{\bf x} \rangle^{1/2}\)
r = b - A.dot(x)
np.sqrt(r.dot(r))
5.0242958677880805e-15
Como a norma do resíduo é próxima de zero, a solução do sistema linear é assumida como correta.
8.5. Eliminação Gaussiana#
A Eliminação Gaussiana (EG) é um algoritmo utilizado para resolver sistemas de equações lineares ao reduzir a matriz plena associada do sistema a uma matriz triangular. Este processo também é chamado de escalonamento. Abaixo, usaremos uma matriz genérica 3x3 para exemplificação.
Os passos para a EG são os seguintes:
Escrever o sistema linear na forma de matriz estendida usando os coeficientes das variáveis como elementos da matriz e o vetor independente como sendo a última coluna;
Realizar operações elementares de combinação linear e permutação entre linhas;
Multiplicação por escalar:
$$ \begin{array}{c} L_2 \leftarrow L_2 .w\ \Rightarrow\
\end{array} $$
Combinação linear:
$$ \begin{array}{c} L_2 \leftarrow L_2 - L_1.w\ \Rightarrow\
\end{array} $$
Permutação:
Anular todos os elementos na porção triangular inferior da matriz original, isto é, todas as entradas exatamente abaixo das entradas dispostas na diagonal principal;
A partir da forma triangular, realizar a substituição regressiva.
Vejamos um exemplo numérico de como funciona a Eliminação Gaussiana.
# matriz
M = np.array([[1.0,1.5,-2.0],[2.0,1.0,-1.0],[3.0,-1.0,2.0]])
print(M)
[[ 1. 1.5 -2. ]
[ 2. 1. -1. ]
[ 3. -1. 2. ]]
# zeramento da segunda linha
m1 = M[1,0]/M[0,0]
M[1,:] += - m1*M[0,:]
print(M)
[[ 1. 1.5 -2. ]
[ 0. -2. 3. ]
[ 3. -1. 2. ]]
# zeramento da terceira linha
m2 = M[2,0]/M[0,0]
M[2,:] += - m2*M[0,:]
print(M)
[[ 1. 1.5 -2. ]
[ 0. -2. 3. ]
[ 0. -5.5 8. ]]
# eliminação
M = np.array([[1.0,1.5,-2.0],[2.0,1.0,-1.0],[3.0,-1.0,2.0]])
b = np.array([-2,3,1])
m,n = M.shape
for i in range(m):
for j in range(i+1,n):
pivo = M[j,i]/M[i,i]
for k in range(m):
M[j,k] += -pivo*M[i,k]
print(M)
[[ 1. 1.5 -2. ]
[ 0. -2. 3. ]
[ 0. 0. -0.25]]
# função simples para eliminação
def eliminacao(M):
m,n = M.shape
for i in range(m):
for j in range(i+1,n):
pivo = M[j,i]/M[i,i]
for k in range(m):
M[j,k] += -pivo*M[i,k]
return M
# matriz randômica 5x5
M2 = np.random.rand(5,5)
print(eliminacao(M2))
[[ 2.17290028e-01 1.53718524e-01 3.21754779e-01 3.57012363e-01
6.51351085e-01]
[ 2.77555756e-17 1.05316462e-01 3.77045560e-01 -2.32541306e-01
-6.30255349e-01]
[-3.41761440e-17 0.00000000e+00 -3.75397868e-01 -3.83992779e-01
-7.23027923e-01]
[ 3.90760482e-17 0.00000000e+00 -2.77555756e-17 -7.46062713e-01
-1.94455006e+00]
[-1.00315604e-16 0.00000000e+00 5.47203352e-17 0.00000000e+00
-7.23561110e-01]]
8.6. Condicionamento de matrizes#
O condicionamento de uma matriz está relacionado à sensibilidade da solução de um sistema linear \(\mathbf{A}\mathbf{x} = \mathbf{b}\) a pequenas perturbações nos dados, sejam elas no vetor \(\mathbf{b}\) ou nos coeficientes da própria matriz \(\mathbf{A}\). Uma matriz é considerada bem condicionada quando pequenas variações em \(\mathbf{b}\) produzem pequenas variações na solução \(\mathbf{x}\); caso contrário, diz-se que a matriz é mal condicionada, o que compromete a estabilidade numérica e a confiabilidade da solução computacional.
O número de condição da matriz, denotado por \(\kappa(\mathbf{A})\), é uma medida quantitativa desse comportamento, sendo definido, por exemplo, como
para uma dada norma. Quanto maior esse número, pior o condicionamento. Em problemas de engenharia e ciências aplicadas, identificar o condicionamento da matriz é fundamental para avaliar a precisão da solução e escolher algoritmos adequados, como métodos com pivotamento ou regularizações.
Vejamos como pequenas “perturbações” em matrizes podem provocar mudanças drásticas nas soluções de sistemas. Isto ocorre quando temos um problema mal condicionado.
A1 = np.array([[1,2],[1.1,2]])
b1 = np.array([10,10.4])
print('matriz')
print(A1)
print('vetor')
print(b1)
matriz
[[1. 2. ]
[1.1 2. ]]
vetor
[10. 10.4]
# solução do sistema A1x1 = b1
x1 = linalg.solve(A1,b1)
print(x1)
[4. 3.]
d = 0.045
A2 = np.array([[1,2],[1.1-d,2]])
b2 = np.array([10,10.4])
print('matriz')
print(A2)
print('vetor')
print(b2)
matriz
[[1. 2. ]
[1.055 2. ]]
vetor
[10. 10.4]
# solução do sistema perturbado A2x1 = b2
x2 = linalg.solve(A2,b2)
print(x2)
[7.27272727 1.36363636]
A solução muda drasticamente aqui! Isto se deve à quase dependência linear em que a matriz se encontra. Ou seja, \(\det({\bf A}_2) \approx 0\).
print(linalg.det(A1),linalg.det(A2))
-0.20000000000000018 -0.11000000000000032
linalg.norm(A)*linalg.norm(linalg.inv(A))
169.2838804582743
8.7. Método de Gauss-Jordan#
O método de Gauss-Jordan é uma variação da EG, em que não apenas as entradas da porção inferior da matriz plena do sistema são anuladas, mas também as entradas da porção superior, resultando em uma matriz diagonal.
Além disso, todas as linhas são normalizadas através da sua divisão pelo respectivo elemento pivô. Por exemplo, \(a_{11}\) é o elemento pivô da primeira equação, \(a_{22}\) da segunda, e assim por diante). A partir daí, a obtenção dos valores das variáveis é imediata.
O método é melhor ilustrado no seguinte exemplo.
# Matriz aumentada
AB = np.array([[3. , -0.1, -0.2, 7.85], [0.1, 7., -0.3, -19.3], [0.3, -0.2, 10., 71.4]])
print(AB)
[[ 3. -0.1 -0.2 7.85]
[ 0.1 7. -0.3 -19.3 ]
[ 0.3 -0.2 10. 71.4 ]]
# Normalização da 1a. linha
AB[0,:] /= AB[0,0] # L1 <- L1/a11
print(AB)
[[ 1.00000000e+00 -3.33333333e-02 -6.66666667e-02 2.61666667e+00]
[ 1.00000000e-01 7.00000000e+00 -3.00000000e-01 -1.93000000e+01]
[ 3.00000000e-01 -2.00000000e-01 1.00000000e+01 7.14000000e+01]]
# Eliminação de x1 da 2a. e 3a. linhas
m1 = AB[1,0]
AB[1,:] -= m1*AB[0,:] # L2 <- L2 - m1*L1
m2 = AB[2,0]
AB[2,:] -= m2*AB[0,:] # L3 <- L3 - m2*L1
print(AB)
[[ 1.00000000e+00 -3.33333333e-02 -6.66666667e-02 2.61666667e+00]
[ 0.00000000e+00 7.00333333e+00 -2.93333333e-01 -1.95616667e+01]
[ 0.00000000e+00 -1.90000000e-01 1.00200000e+01 7.06150000e+01]]
# Normalização da 2a. linha
AB[1,:] /= AB[1,1] # L2 <- L2/a22
print(AB)
[[ 1.00000000e+00 -3.33333333e-02 -6.66666667e-02 2.61666667e+00]
[ 0.00000000e+00 1.00000000e+00 -4.18848168e-02 -2.79319372e+00]
[ 0.00000000e+00 -1.90000000e-01 1.00200000e+01 7.06150000e+01]]
# Eliminação de x2 da 1a. e 3a. linhas
m3 = AB[0,1]
AB[0,:] -= m3*AB[1,:] # L1 <- L1 - m3*L2
m4 = AB[2,1]
AB[2,:] -= m4*AB[1,:] # L3 <- L3 - m4*L2
print(AB)
[[ 1.00000000e+00 0.00000000e+00 -6.80628272e-02 2.52356021e+00]
[ 0.00000000e+00 1.00000000e+00 -4.18848168e-02 -2.79319372e+00]
[ 0.00000000e+00 0.00000000e+00 1.00120419e+01 7.00842932e+01]]
# Normalização da 3a. linha
AB[2,:] /= AB[2,2] # L3 <- L3/a33
print(AB)
[[ 1. 0. -0.06806283 2.52356021]
[ 0. 1. -0.04188482 -2.79319372]
[ 0. 0. 1. 7. ]]
# Eliminação de x3 da 1a. e 2a. linhas
m5 = AB[0,2]
AB[0,:] -= m5*AB[2,:] # L1 <- L1 - m5*L3
m6 = AB[1,2]
AB[1,:] -= m6*AB[2,:] # L2 <- L2 - m5*L3
print(AB)
[[ 1. 0. 0. 3. ]
[ 0. 1. 0. -2.5]
[ 0. 0. 1. 7. ]]
Do último resultado, vemos que a matriz identidade é obtida, apontando para o vetor solução \([3 \ \ -2.5 \ \ 7]^T\).