Antes de codar, vamos recapitular a teoria que cada método usa:
LU: a gente escreve , onde é triangular inferior com 1’s na diagonal e é triangular superior. Daí, para resolver , resolvemos primeiro (forward-substitution) e depois (backward-substitution). Pra achar , basta resolver coluna a coluna.
Cholesky: se é simétrica e definida positiva (SPD), existe triangular tal que . Aí a solução segue o mesmo esquema: a) ; b) .
Thomas: é a LU especializada para matrizes tridiagonais. Mantemos apenas três vetores (diagonal inferior , principal e superior ) e rodamos um “eliminador” só pra esses três.
Gauss–Seidel: é iterativo. A cada iteração, atualizamos
até alcançar convergência (diferença entre iterações abaixo de certa tol).
Passo 3
Implementando LU do zero (Python + NumPy só pra armazenar matriz):
import numpy as np
deflu_decomposition(A):
n = A.shape[0]
L = np.eye(n)
U = A.copy().astype(float)
for k inrange(n-1):
for i inrange(k+1, n):
if U[k,k] == 0:
raise ValueError("Pivô zero!")
L[i,k] = U[i,k]/U[k,k]
U[i,k:] = U[i,k:] - L[i,k]*U[k,k:]
return L, U
deflu_solve(L, U, b):
n = len(b)
# forward
y = np.zeros(n)
for i inrange(n):
y[i] = b[i] - L[i,:i] @ y[:i]
# back
x = np.zeros(n)
for i inreversed(range(n)):
x[i] = (y[i] - U[i,i+1:] @ x[i+1:]) / U[i,i]
return x
definvert_via_lu(A):
n = A.shape[0]
L, U = lu_decomposition(A)
invA = np.zeros((n,n))
I = np.eye(n)
for j inrange(n):
invA[:,j] = lu_solve(L, U, I[:,j])
return invA
Passo 4
Cholesky (para SPD):
defcholesky_decomposition(A):
n = A.shape[0]
L = np.zeros_like(A)
for i inrange(n):
for j inrange(i+1):
s = A[i,j] - L[i,:j] @ L[j,:j]
if i==j:
if s<=0:
raise ValueError("Matriz não SPD")
L[i,j] = np.sqrt(s)
else:
L[i,j] = s / L[j,j]
return L
defcholesky_solve(L, b):
# L y = b
n = len(b)
y = np.zeros(n)
for i inrange(n):
y[i] = (b[i] - L[i,:i] @ y[:i]) / L[i,i]
# L^T x = y
x = np.zeros(n)
for i inreversed(range(n)):
x[i] = (y[i] - L[i+1:,i] @ x[i+1:]) / L[i,i]
return x
Passo 5
Thomas (sistema tridiagonal ):
defthomas(a, b, c, d):
n = len(d)
# cópia para não alterar vetores originais
c_ = c.copy()
d_ = d.copy()
# eliminacaofor i inrange(1, n):
m = a[i]/b[i-1]
b[i] = b[i] - m*c_[i-1]
d_[i] = d_[i] - m*d_[i-1]
# back
x = np.zeros(n)
x[-1] = d_[-1]/b[-1]
for i inreversed(range(n-1)):
x[i] = (d_[i] - c_[i]*x[i+1]) / b[i]
return x
Passo 6
Gauss–Seidel (iterativo):
defgauss_seidel(A, b, x0=None, tol=1e-8, max_iter=1000):
n = len(b)
x = np.zeros(n) if x0 isNoneelse x0.copy()
for k inrange(max_iter):
x_new = x.copy()
for i inrange(n):
soma1 = A[i,:i] @ x_new[:i]
soma2 = A[i,i+1:] @ x[i+1:]
x_new[i] = (b[i] - soma1 - soma2)/A[i,i]
if np.linalg.norm(x_new - x, ord=np.inf) < tol:
return x_new, k+1
x = x_new
raise ValueError("Não convergiu em {} iterações".format(max_iter))
Passo 7
Menu amigável e teste num exemplo tridiagonal SPD (por exemplo ):
defmain():
A = np.array([[2,-1,0],[-1,2,-1],[0,-1,2]], float)
b = np.array([1,0,1], float)
print("Matriz A:\n", A)
print("Vetor b:", b)
# (a) LU
L, U = lu_decomposition(A)
x_lu = lu_solve(L, U, b)
invA = invert_via_lu(A)
print("\nLU: x =", x_lu)
print("Inversa A^-1:\n", invA)
# (b) Cholesky
Lc = cholesky_decomposition(A)
x_ch = cholesky_solve(Lc, b)
print("\nCholesky: x =", x_ch)
# (c) Thomas
a = np.array([0, -1, -1], float) # a[0] não usado
bb = np.array([2, 2, 2], float)
c = np.array([-1, -1, 0], float) # c[-1] não usado
x_th = thomas(a, bb, c, b)
print("\nThomas: x =", x_th)
# (d) Gauss-Seidel
x_gs, its = gauss_seidel(A, b)
print("\nGauss–Seidel ({} it.): x =".format(its), x_gs)
if __name__=="__main__":
main()
Resposta
Rodando esse programa, a gente obtém aproximadamente:
LU: x = [1.5, 2.0, 2.5]
A⁻¹ = 0.75&0.5&0.250.5&1.0&0.50.25&0.5&0.75
Cholesky: x = [1.5, 2.0, 2.5]
Thomas: x = [1.5, 2.0, 2.5]
Gauss–Seidel: x = [1.5, 2.0, 2.5] (convergiu em ~5 iterações)