Desenvolva um programa amigável ao usuário para solução de sistemas algébricos lineares que utilizem:
(a) Decomposição LU, incluindo o recurso de calcular a matriz inversa
(b) Decomposição de Cholesky
© Algoritmo de Thomas
(d) Método de Gauss-Seidel
Teste o programa em um dos exemplos de aplicação.

Passo 1

Beleza, então o exercício pede pra gente criar um programinha “amigável” (em Python, por exemplo) que resolva sistemas lineares usando quatro jeitos diferentes:
(a) Decomposição LU (incluindo cálculo da inversa)
(b) Decomposição de Cholesky
© Algoritmo de Thomas (sistemas tridiagonais)
(d) Método de Gauss–Seidel (iterativo)
Pra isso a gente vai precisar de noções básicas de álgebra linear (fatoração de matrizes, solução de sistemas triangular, condição de SPD para Cholesky, etc.) e saber programar loops, funções e menus.

Passo 2

Antes de codar, vamos recapitular a teoria que cada método usa:
  1. LU: a gente escreve xi^(k+1) = 1aiil(bi - sumji aij,xj^(k)r) , 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.
  2. Cholesky: se é simétrica e definida positiva (SPD), existe triangular tal que . Aí a solução segue o mesmo esquema: a) ; b) .
  3. 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.
  4. 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

def lu_decomposition(A):
    n = A.shape[0]
    L = np.eye(n)
    U = A.copy().astype(float)
    for k in range(n-1):
        for i in range(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

def lu_solve(L, U, b):
    n = len(b)
    # forward
    y = np.zeros(n)
    for i in range(n):
        y[i] = b[i] - L[i,:i] @ y[:i]
    # back
    x = np.zeros(n)
    for i in reversed(range(n)):
        x[i] = (y[i] - U[i,i+1:] @ x[i+1:]) / U[i,i]
    return x

def invert_via_lu(A):
    n = A.shape[0]
    L, U = lu_decomposition(A)
    invA = np.zeros((n,n))
    I = np.eye(n)
    for j in range(n):
        invA[:,j] = lu_solve(L, U, I[:,j])
    return invA

Passo 4

Cholesky (para SPD):
def cholesky_decomposition(A):
    n = A.shape[0]
    L = np.zeros_like(A)
    for i in range(n):
        for j in range(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

def cholesky_solve(L, b):
    # L y = b
    n = len(b)
    y = np.zeros(n)
    for i in range(n):
        y[i] = (b[i] - L[i,:i] @ y[:i]) / L[i,i]
    # L^T x = y
    x = np.zeros(n)
    for i in reversed(range(n)):
        x[i] = (y[i] - L[i+1:,i] @ x[i+1:]) / L[i,i]
    return x

Passo 5

Thomas (sistema tridiagonal ):
def thomas(a, b, c, d):
    n = len(d)
    # cópia para não alterar vetores originais
    c_ = c.copy()
    d_ = d.copy()
    # eliminacao
    for i in range(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 in reversed(range(n-1)):
        x[i] = (d_[i] - c_[i]*x[i+1]) / b[i]
    return x

Passo 6

Gauss–Seidel (iterativo):
def gauss_seidel(A, b, x0=None, tol=1e-8, max_iter=1000):
    n = len(b)
    x = np.zeros(n) if x0 is None else x0.copy()
    for k in range(max_iter):
        x_new = x.copy()
        for i in range(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
2 & -1 & 0-1 & 2 & -10 & -1 & 2):
def main():
    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)

Exercícios relacionados