3.2 Métodos diretos: decomposição LU#
Os métodos diretos para resolução de sistemas lineares teoricamente dão a solução exata do sistema em um número finito de passos. Na prática, as soluções podem sofrer erros de arredondamento devido à aritmética que está sendo usada.
Os métodos diretos de eliminação consistem no processo de transformar um sistema de equações lineares \(Ax=b\), em que \(A= [a_{ij}]_i,j=1,..,n\), \(x=[x_1,x_2,...,x_n]^t\) e \(b=[b_1, b_2,...,b_n]^t\), em um sistema equivalente, aplicando operações elementares sobre as linhas da matriz aumentada \([A \mid b]\). O sistema equivalente obtido ao final do processo, deve ser de resolução mais simples ou imediata, por substituição direta.
O método de eliminação de Gauss é o método direto usado para escalonar a matriz que representa um sistema linear com o objetivo de obter um ou mais sistemas triangulares, os quais podem ser resolvidos por subtituições subsequentes. Outros métodos diretos podem envolver fatorar a matriz dos coeficientes do sistema em um produto de duas ou mais matrizes, como é o caso do método por fatoração LU, que será abordado nesse texto.
Antes disso, vamos relembrar as função já desenvolvidas nas aulas anteriores.
Funções para escalonar uma matriz e para resolver sistemas triangulares#
import numpy as np
# função para escalonar (sem pivoteamento)
def escalona(M):
n = len(M)
for j in range(n-1):
for i in range(j,n-1):
m = M[i+1,j]/M[j,j]
M[i+1] = M[i+1]-m*M[j]
return (M)
# função para sistema triangular inferior
def solve_tril(L,b):
y = np.empty(len(b))
for i in range(len(b)):
y[i] = (b[i] - np.dot(L[i,0:i],y[0:i]))/L[i][i]
return y
# função para sistema triangular superior
def solve_triu(U,y):
n=len(y)
x = np.empty(n)
for i in range(n-1, -1, -1):
x[i] = (y[i] - np.sum(U[i,i+1:n+1]*x[i+1:n+1]))/U[i,i]
return x
Exemplo 1: Vamos neste exemplo revisar o sistema linear abaixo, já estudado anteriormente (fonte). $\( \begin{aligned} x_1-x_2+2 x_3-x_4 & =-8, \\ 2 x_1-2 x_2+3 x_3-3 x_4 & =-20, \\ x_1+x_2+x_3 & =-2, \\ x_1-x_2+4 x_3+3 x_4 & =4 . \end{aligned} \)$
Para isso, vamos analisar o escalonamento, coluna por coluna.
# escrevendo a matriz ampliada
M = np.array([[1,-1,2,-1,-8],
[2,-2,3,-3,-20],
[1,1,1,0,-2],
[1,-1,4,3,4]], dtype=float)
# zerando os elementos da primeira coluna
for i in range(3):
m = M[i+1,0]/M[0,0]
M[i+1] = M[i+1] - m*M[0]
print (M)
[[ 1. -1. 2. -1. -8.]
[ 0. 0. -1. -1. -4.]
[ 0. 2. -1. 1. 6.]
[ 0. 0. 2. 4. 12.]]
Observemos que o elemento pivô na segunda linha é zero, portanto, para que o processo de eliminação possa ser continuado, é necessário realizar uma troca de linhas.
M[[1,2]] = M[[2,1]]
print(M)
[[ 1. -1. 2. -1. -8.]
[ 0. 2. -1. 1. 6.]
[ 0. 0. -1. -1. -4.]
[ 0. 0. 2. 4. 12.]]
Agora podemos seguir o escalonamento na terceira coluna.
# zerando os elementos da terceira coluna
m = M[3,2]/M[2,2]
M[3] = M[3] - m*M[2]
print (M)
[[ 1. -1. 2. -1. -8.]
[ 0. 2. -1. 1. 6.]
[ 0. 0. -1. -1. -4.]
[ 0. 0. 0. 2. 4.]]
# solução do sistema
x = solve_triu(M[:,:-1], M[:,-1])
print(x)
[-7. 3. 2. 2.]
Uma função bastante útil para tratar esse tipo de situação é a função numpy.roll(), que “rola” os elementos do array ao longo do eixo especificado. Se o elemento a ser desclocado estiver na última posição, ele será levado de volta para a primeira posição. Veja um exemplo a seguir.
# rolando as linhas de M para baixo uma vez
np.roll(M, 1, axis=0)
array([[ 0., 0., 0., 2., 4.],
[ 1., -1., 2., -1., -8.],
[ 0., 2., -1., 1., 6.],
[ 0., 0., -1., -1., -4.]])
# rolando as linhas de M para cima duas vezes
np.roll(M, -2, axis=0)
array([[ 0., 0., -1., -1., -4.],
[ 0., 0., 0., 2., 4.],
[ 1., -1., 2., -1., -8.],
[ 0., 2., -1., 1., 6.]])
Exercício 1:
Use a função roll para modificar a função escalona mostrada abaixo, de modo que elementos nulos na posição dos pivôs sejam substituídos pela troca de linhas da matriz.
Use a matriz do Exemplo 1 para testar.
M = np.array([[1,-1,2,-1,-8],
[2,-2,3,-3,-20],
[1,1,1,0,-2],
[1,-1,4,3,4]], dtype=float)
def escalona2(M):
n = len(M)
for j in range(n-1):
for i in range(j,n-1):
# rolando as linhas
while M[j,j]==0:
M[j:,:] = np.roll(M[j:,:], -1, axis=0)
m = M[i+1,j]/M[j,j]
M[i+1] = M[i+1]-m*M[j]
return (M)
M = escalona2(M)
print(M)
[[ 1. -1. 2. -1. -8.]
[ 0. 2. -1. 1. 6.]
[ 0. 0. 2. 4. 12.]
[ 0. 0. 0. 1. 2.]]
# solução do sistema
x = solve_triu(M[:,:-1], M[:,-1])
print(x)
[-7. 3. 2. 2.]
Observe que a matriz escalonada é diferente da matriz obtida no Exemplo 1 pois a troca de linha ocorreu de forma diferente durante o processo. Mesmo assim, o resultado final é o mesmo, como já era esperado.
Exercício 2: Resolva o seguinte sistema pelo método de eliminação de Gauss.
M = np.array([[2,1,-1,1,-3,7],
[1,0,2,-1,1,2],
[0,-2,-1,1,-1,-5],
[3,1,-4,0,5,6],
[1,-1,-1,-1,1,3]], dtype=float)
Mt = escalona(M.copy())
np.round(Mt,3)
array([[ 2. , 1. , -1. , 1. , -3. , 7. ],
[ 0. , -0.5 , 2.5 , -1.5 , 2.5 , -1.5 ],
[ 0. , 0. , -11. , 7. , -11. , 1. ],
[ 0. , 0. , 0. , -3.182, 12. , -3.455],
[ 0. , 0. , 0. , 0. , -4.886, 5.543]])
# solução do sistema
x = solve_triu(Mt[:,:-1], Mt[:,-1])
print(x)
[ 1.91812865 1.96491228 -0.98830409 -3.19298246 -1.13450292]
M
array([[ 2., 1., -1., 1., -3., 7.],
[ 1., 0., 2., -1., 1., 2.],
[ 0., -2., -1., 1., -1., -5.],
[ 3., 1., -4., 0., 5., 6.],
[ 1., -1., -1., -1., 1., 3.]])
M[:,:-1]@x
array([ 7., 2., -5., 6., 3.])
Exercício 3: Use o algoritmo de eliminação de Gauss para resolver os seguintes sistemas lineares, se possível, e determine se são necessárias trocas de linhas (fonte)
a)\( \begin{aligned} x_1+x_2+x_4 & =2, \\ 2 x_1+x_2-x_3+x_4 & =1, \\ 4 x_1-x_2-2 x_3+2 x_4 & =0, \\ 3 x_1-x_2-x_3+2 x_4 & =-3 . \end{aligned} \quad \quad\) b) \( \begin{aligned} x_1-\frac{1}{2} x_2+x_3 & =4 \\ 2 x_1-x_2-x_3+x_4 & =5 \\ x_1+x_2+\frac{1}{2} x_3 & =2 \\ x_1-\frac{1}{2} x_2+x_3+x_4 & =5 \end{aligned} \quad \quad\) c) \(\begin{aligned} x_1+x_2+x_4 & =2, \\ 2 x_1+x_2-x_3+x_4 & =1, \\ -x_1+2 x_2+3 x_3-x_4 & =4, \\ 3 x_1-x_2-x_3+2 x_4 & =-3 . \end{aligned}\)
M1 = np.array([[1,1,0,1,2],[2,1,-1,1,1],[4,-1,-2,2,0],[3,-1,-1,2,-3]],dtype=float)
M2 = np.array([[1,-0.5,1,0,4],[2,-1,-1,1,5],[1,1,0.5,0,2],[1,-0.5,1,1,5]],dtype=float)
M3 = np.array([[1,1,0,1,2],[2,1,-1,1,1],[-1,2,3,-1,4],[3,-1,-1,2,-3]],dtype=float)
Solução letra a):
M1e = escalona(M1)
M1e
array([[ 1., 1., 0., 1., 2.],
[ 0., -1., -1., -1., -3.],
[ 0., 0., 3., 3., 7.],
[ 0., 0., 0., 0., -4.]])
Sistema impossível.
Solução letra b):
M2e = escalona(M2.copy())
M2e
/tmp/ipykernel_4617/899190976.py:8: RuntimeWarning: divide by zero encountered in scalar divide
m = M[i+1,j]/M[j,j]
/tmp/ipykernel_4617/899190976.py:9: RuntimeWarning: invalid value encountered in multiply
M[i+1] = M[i+1]-m*M[j]
/tmp/ipykernel_4617/899190976.py:8: RuntimeWarning: invalid value encountered in scalar divide
m = M[i+1,j]/M[j,j]
array([[ 1. , -0.5, 1. , 0. , 4. ],
[ 0. , 0. , -3. , 1. , -3. ],
[ nan, nan, inf, -inf, inf],
[ nan, nan, nan, nan, nan]])
M2e = escalona2(M2.copy())
M2e
array([[ 1. , -0.5, 1. , 0. , 4. ],
[ 0. , 1.5, -0.5, 0. , -2. ],
[ 0. , 0. , -3. , 1. , -3. ],
[ 0. , 0. , 0. , 1. , 1. ]])
# solução do sistema
x = solve_triu(M2e[:,:-1], M2e[:,-1])
print(x)
[ 2.22222222 -0.88888889 1.33333333 1. ]
M2[:,:-1]@x
array([4., 5., 2., 5.])
Solução letra c):
M3
array([[ 1., 1., 0., 1., 2.],
[ 2., 1., -1., 1., 1.],
[-1., 2., 3., -1., 4.],
[ 3., -1., -1., 2., -3.]])
M3e = escalona2(M3.copy())
M3e
array([[ 1., 1., 0., 1., 2.],
[ 0., -1., -1., -1., -3.],
[ 0., 0., 3., 3., 3.],
[ 0., 0., 0., -3., -3.]])
# solução do sistema
x = solve_triu(M3e[:,:-1], M3e[:,-1])
print(x)
[-1. 2. 0. 1.]
M3[:,:-1]@x
array([ 2., 1., 4., -3.])
Decomposição LU#
A fatoração \(LU\) é útil para resolver sistemas lineares, calcular determinantes e encontrar inversas de matrizes de forma mais eficiente. Uma vez que a matriz \(A\) é fatorada, resolver sistemas do tipo \(Ax=b\) fica mais rápido, pois basta resolver dois sistemas triangulares mais simples, primeiramente \(Ly=b\) e, então, \(Ux=y\).
Além disso, é mais fácil obter o determinante de \(A\) calculando o produto dos elementos da diagonal da matriz triangular. Por fim, evita fazer a eliminação gaussiana para obter a solução do sistema \(Ax=b\) toda vez que muda o vetor \(b\).
As mesmas operações com linhas utilizadas para resolver um sistema da forma \(Ax = b\) por eliminação de Gauss podem ser utilizados para fatorar uma matriz \(A\) em um produto de duas matrizes triangulares \(L\) (triangular inferior) e \(U\) (triangular superior). Vejamos um exemplo:
Exemplo 3: Seja o sistema
Vamos encontrar a solução usando o processo de fatoração \(LU\).
# criamos as matrizes
A = np.array([[2,0,1],[0,2,1],[1,1,3]], dtype=float)
b = np.array([3,3,5], dtype=float)
Então, vamos elimintar os elementos da 1ª coluna, abaixo do pivô:
# segunda linha
m21 = A[1,0]/A[0,0]
A[1] = A[1] - m21*A[0]
#terceira linha
m31 = A[2,0]/A[0,0]
A[2] = A[2] - m31*A[0]
print (A)
[[2. 0. 1. ]
[0. 2. 1. ]
[0. 1. 2.5]]
Observe que o resultado é o mesmo que se multiplicar a matriz
$\(M_1= \left[\begin{array}{lll} 1 & 0 & 0 \\ -m21 & 1 & 0 \\ -m31 & 0 & 1 \end{array}\right]= \left[\begin{array}{lll} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -0.5 & 0 & 1 \end{array}\right] \)$.
pela matriz \(A\) como é mostrado abaixo
M1 = np.array([[1,0,0],[0,1,0],[-0.5,0,1]])
A =np.array([[2,0,1],[0,2,1],[1,1,3]])
A = M1@A
print(A)
[[2. 0. 1. ]
[0. 2. 1. ]
[0. 1. 2.5]]
Por fim, concluímos a eliminação de Gauss fazendo:
#terceira linha segunda coluna
m32 = A[2,1]/A[1,1]
A[2] = A[2] - m32*A[1]
print (A)
[[2. 0. 1.]
[0. 2. 1.]
[0. 0. 2.]]
O resultado é o mesmo que obtemos ao multiplicar a matriz $\(M_2 = \left[\begin{array}{lll} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -m_{32} & 1 \end{array}\right]= \left[\begin{array}{lll} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -0.5 & 1 \end{array}\right] \)$
por \(M_1 A\).
Ou seja, o escalonamento da matriz \(A\) pode ser realizado pelo produto \(M_2M_1 A\), como podemos verificar a seguir:
A =np.array([[2,0,1],[0,2,1],[1,1,3]])
M1 = np.array([[1,0,0],[0,1,0],[-0.5,0,1]])
M2 = np.array([[1,0,0],[0,1,0],[0,-0.5,1]])
U = M2@M1@A
print(U)
[[2. 0. 1.]
[0. 2. 1.]
[0. 0. 2.]]
Então, em resumo, tem-se:
Observe que podemos obter a matriz \(A\) a partir da matriz triangular superior \(U\) fazendo \(M_1^{-1}M_2^{-1}U=A\), pois
from numpy.linalg import inv
inv(M1)@inv(M2)@U
array([[2., 0., 1.],
[0., 2., 1.],
[1., 1., 3.]])
Assim, fazendo
podemos escrever a fatoração \(LU\) como
Como é verificado abaixo
L = np.array([[1,0,0],[0,1,0],[0.5,0.5,1]])
U = np.array([[2,0,1],[0,2,1],[0,0,2]])
L@U
array([[2., 0., 1.],
[0., 2., 1.],
[1., 1., 3.]])
Por fim, para obter a solução do sistema linear \(Ax=b\) por decomposição \(LU\) fazemos
escrevendo \(Ux=y\) obtemos a solução resolvendo dois sitemas triangulares
Assim, a resolução do sistema \(Ax=b\) é feita resolvendo-se o sistema triangular inferior e, em seguida, o sistema triagular superior como é mostrado a seguir:
y = solve_tril(L,b)
print(y)
[3. 3. 2.]
x = solve_triu(U,y)
print(x)
[1. 1. 1.]
Agora, podemos sintetizar os resultados escrevendo uma função para fatorar uma matriz por eliminação de Gauss , fazendo:
def fatoraLU(A):
L = np.eye(len(A))
U = A.copy()
n = len(A)
for j in range(n-1):
for i in range(j,n-1):
m = U[i+1,j]/U[j,j]
U[i+1] = U[i+1]-m*U[j]
L[i+1,j] = m
return (L,U)
# testando
A = np.array([[2,0,1],[0,2,1],[1,1,3]], dtype=float)
L,U = fatoraLU(A)
print(L)
print(U)
[[1. 0. 0. ]
[0. 1. 0. ]
[0.5 0.5 1. ]]
[[2. 0. 1.]
[0. 2. 1.]
[0. 0. 2.]]
y = solve_tril(L,b)
print("y=",y)
y= [3. 3. 2.]
x = solve_triu(U,y)
print("x=",x)
x= [1. 1. 1.]
Exemplo 3: Agora vamos considerar novamente o sistema do exemplo 1 e resolver usando fatoração LU.
M = np.array([[2,1,-1,1,-3,7],
[1,0,2,-1,1,2],
[0,-2,-1,1,-1,-5],
[3,1,-4,0,5,6],
[1,-1,-1,-1,1,3]], dtype=float)
A = M[:,0:5]
b = M[:,-1]
L,U = fatoraLU(A)
y = solve_tril(L,b)
x = solve_triu(U,y)
print(x)
[ 1.91812865 1.96491228 -0.98830409 -3.19298246 -1.13450292]
A@x
array([ 7., 2., -5., 6., 3.])
Fatoração LU por redução de Doolittle#
Uma forma de obter a fatoração \(LU\) da matriz \(A\), conhecido como processo Redução de Doolittle será apresentado a seguir (fonte). Considere a fatoração da matriz \(A = (a_{ij})_{i,j=1,...,n}\) nas matrizes \(L = (l_{ij})_{i,j=1,...,n}\) e \(U = (u_{ij})_{i,j=1,...,n}\) como é mostrado a seguir
o processo de múltiplicação de matrizes nos leva às seguintes equações, que fornecem a primeira linha da matriz \(U\): \(u_{11}=a_{11}\), \(u_{12}=a_{12}\), \(u_{12}=a_{12}\),…, \(u_{1n}=a_{1n}\). Seguindo o procedimento de multiplicar matrizes, podemos encontrar a 1ª coluna da matriz \(L\): \(l_{21}=a_{21}/u_{11}\), \(l_{31}=a_{31}/u_{11}\),…, \(u_{n1}=a_{n1}/u_{11}\). Continuando o processo para a 2ª linha de \(U\), 2ª coluna de \(L\), 3ª linha de \(U\), 3ª coluna de \(L\), e assim por diante, chegamos nas seguintes fórmulas: $\(u_{ij}=a_{ij}-\sum\limits_{k=1}^{i-1}l_{ik}u_{kj} \qquad i,j=1,...,n\)\( e \)\(l_{ij}=\frac{a_{ij}-\sum\limits_{k=1}^{j-1}l_{ik}u_{kj}}{ujj} \qquad i,j=1,...,n\)$
Uma implementação desse algoritmo é mostrado a seguir.
def fatoracao_LU(A):
n = len(A)
U = np.zeros((n,n))
L = np.identity(n)
for m in range(n):
for j in range(m, n):
U[m,j] = A[m,j] - np.sum(L[m,0:m] * U[0:m,j])
for i in range(m+1, n):
L[i,m] = (A[i,m] - np.sum(L[i,0:m] * U[0:m,m]))/U[m,m]
return L,U
Testando para a mesma matriz do exemplo 3.
A = np.array([[2,0,1],
[0,2,1],
[1,1,3]])
L,U = fatoracao_LU(A)
print ("Matriz L:")
print (np.array(L))
print ("Matriz U:")
print(np.array(U))
Matriz L:
[[1. 0. 0. ]
[0. 1. 0. ]
[0.5 0.5 1. ]]
Matriz U:
[[2. 0. 1.]
[0. 2. 1.]
[0. 0. 2.]]
Exercícios (fonte):
1. Determine a fatoração \(LU\) da matriz \(A\) e resolva o sistema linear \(Ax = b\), em que
2. Use a fatoração do exercício anterior para resolver o sistema $\( \begin{aligned} x_1+x_2+3 x_4 & =8, \\ 2 x_1+x_2-x_3+x_4 & =7, \\ 3 x_1-x_2-x_3+2 x_4 & =14, \\ -x_1+2 x_2+3 x_3-x_4 & =-7 . \end{aligned} \)$
Usando o módulo scipy.linalg para fatorar uma matriz#
https://docs.scipy.org/doc/scipy/reference/linalg.html
import scipy
import scipy.linalg
A = np.array([[1,2,4],[2,8,10],[4,10,26]])
Fatoração LU
scipy.linalg.lu(A)
(array([[0., 0., 1.],
[0., 1., 0.],
[1., 0., 0.]]),
array([[ 1. , 0. , 0. ],
[ 0.5 , 1. , 0. ],
[ 0.25 , -0.16666667, 1. ]]),
array([[ 4., 10., 26.],
[ 0., 3., -3.],
[ 0., 0., -3.]]))
Fatoração Cholesky
A fatoração de Cholesky é uma forma especial de decompor uma matriz simétrica e definida positiva. Ela escreve a matriz \(A\) como:
onde \(L\) é uma matriz triangular inferior (com números reais positivos na diagonal), \(L^T\) é a transposta de \(L\). Nesse caso a solução de um sistema linear \(Ax=b\) por fatoração de Cholesky pode ser ainda mais rápida do que por fatoração \(LU\).
A decomposição Cholesky é especialmente útil em problemas de otimização, sistemas de equações diferenciais, simulações físicas e finanças.
L = np.linalg.cholesky(A)
L
array([[1., 0., 0.],
[2., 2., 0.],
[4., 1., 3.]])
#verificando
L@L.T
array([[ 1., 2., 4.],
[ 2., 8., 10.],
[ 4., 10., 26.]])
Fatoração QR
A fatoração \(QR\) é uma decomposição que escreve uma matriz \(A\) como o produto de duas matrizes, uma matriz ortogonal \(Q\) e uma matriz triangular superior \(R\). É útil na resolução de sistemas lineares superdeterminados, no cálculo de autovalores e na ortonormalização de vetores. Além disso, a fatoração \(QR\) é mais estável numericamente que a \(LU\), especialmente para sistemas mal condicionados.
Q,R = np.linalg.qr(A)
print(Q)
print(R)
[[-0.21821789 0.27263927 -0.93704257]
[-0.43643578 -0.88607763 -0.15617376]
[-0.87287156 0.374879 0.31234752]]
[[ -4.58257569 -12.65663763 -27.93188995]
[ 0. -2.79455252 1.97663471]
[ 0. 0. 2.81112771]]
Q@R
array([[ 1., 2., 4.],
[ 2., 8., 10.],
[ 4., 10., 26.]])
Exercícios#
1. Obtenha, caso exista, a solução para os seguintes sistemas lineares utilizando os métodos estudados (eliminação de Gauss e decomposição LU).
a) \(\begin{cases} 2x_1 + 3x_2 + x_3 +5x_4= 11\\ x_1 + 3.5x_2 + x_3 +7.5x_4= 13\\ 1.4x_1 + 2.7x_2 + 5.5x_3 + 12x_4 = 21.6\\ -2x_1 + 1x_2 + 3x_3 +28x_4 = 30 \end{cases}\)
b) \(\begin{cases} 6.1x_1 + 0.32x_2 + 1.3x_3 +2.1x_4 + 0.11x_5 = 19.52\\ 0.82x_1 + 8.81x_2 + 1.01x_3 +3x_4 + 3.12x_5= 15.83\\ 0.5x_1 + 1.78x_2 + 15.2x_3 + 4.2x_4 +8.1x_5= -22.14\\ 4.2x_1 + 5.3x_2 + 1.8x_3 +20.9x_4 +7.51x_5 = 27.28\\ 0.2x_1 + 9.1x_2 + 4.68x_3 +4.3x_4 +20.1x_5 = -21.78 \end{cases}\)
c)\( \begin{cases} 12.1756 x_1 + 4.0231 x_2 - 2.1732 x_3 + 5.1967 x_4 = 17.1020\\ -4.0231 x_1 + 6.0030 x_2 + 1.1973 x_4 = -6.1593\\ -1.0000 x_1 - 5.2107 x_2 + 11.1111 x_3 = 3.0004\\ 6.0235 x_1 + 7.0000 x_2 + - 14.1561 x_4 = 0.0000 \end{cases} \)
2. (Burden et al, 2016, p. 420) Resolva o seguinte sistema linear por eliminação de Gauss e decomposição LU (se possível).
3. (Burden et al, 2016, p.455) Use fatoração LU para resolver o sistema linear a seguir.