Python para Sistemas Lineares

Antes de programar os métodos de Gauss, Gauss-Jordan, Jacobi e Gauss-Seidel, é importante conectar a escrita matemática de um sistema linear com as operações matriciais feitas em Python. Como listas, funções, laços e NumPy já foram apresentados na introdução ao Python, aqui o foco será representar sistemas, manipular linhas, controlar pivôs, calcular resíduos e organizar iterações.

Representando um sistema linear

Um sistema linear pode ser escrito na forma matricial:

Nos algoritmos diretos, trabalharemos bastante com a matriz aumentada \([A|b]\). Em Python, uma forma simples de armazenar essa matriz é usar uma lista de listas:

M = [
    [1, -2, -1, 2, 1],
    [-2, 5, 3, -3, 4],
    [-2, 6, 5, -2, 13],
    [2, -4, -1, 5, 7]
]

Cada lista interna representa uma linha da matriz. A última coluna representa o vetor \(b\).

Usando NumPy com tipo real

Para calcular com matrizes de forma mais direta, usaremos a biblioteca numpy. Ela permite multiplicar linhas por constantes, subtrair linhas, calcular produtos matriciais e medir erros com menos código. O uso de dtype=float é importante porque os métodos envolvem divisões e produzem valores reais.

import numpy as np

M = np.array([
    [1, -2, -1, 2, 1],
    [-2, 5, 3, -3, 4],
    [-2, 6, 5, -2, 13],
    [2, -4, -1, 5, 7]
], dtype=float)

Quando for necessário preservar a matriz original para conferência, use uma cópia antes de aplicar eliminações:

M_original = M.copy()
M_trabalho = M.copy()

Índices de linhas e colunas

Na matemática, é comum numerar linhas e colunas começando em 1. Em Python, os índices começam em 0. Assim, a primeira linha é M[0], a segunda linha é M[1] e assim por diante.

Operações com linhas

Nos métodos de Gauss e Gauss-Jordan, a operação mais repetida é substituir uma linha por uma combinação dela com outra linha:

Em Python com NumPy, essa operação pode ser escrita de forma parecida com a expressão matemática:

fator = M[j, i] / M[i, i]
M[j] = M[j] - fator * M[i]

Como essa operação será repetida várias vezes, podemos criar uma função:

def eliminar(M, i, j):
    fator = M[j, i] / M[i, i]
    M[j] = M[j] - fator * M[i]

Nessa função, i é a linha pivô usada como referência para zerar um elemento, enquanto j é a linha alvo que será alterada. A expressão M[j] seleciona a linha inteira, permitindo atualizar todos os elementos da linha de uma vez.

Percorrendo pivôs e linhas alvo

Na eliminação de Gauss, cada coluna define uma linha pivô. Depois de escolher o pivô, percorremos apenas as linhas abaixo dele para zerar os elementos daquela coluna:

dimensao = 4

for i in range(dimensao - 1):
    for j in range(i + 1, dimensao):
        eliminar(M, i, j)

O primeiro laço escolhe a linha pivô. O segundo laço percorre as linhas abaixo dela, que receberão a operação de eliminação.

Pivô e troca de linhas

Antes de dividir pelo pivô, é preciso verificar se ele é nulo ou muito próximo de zero. Uma estratégia comum é trocar a linha atual pela linha abaixo que possui o maior valor absoluto naquela coluna:

coluna = i
linha_maior = i + np.argmax(abs(M[i:, coluna]))

if linha_maior != i:
    M[[i, linha_maior]] = M[[linha_maior, i]]

A expressão M[linha_pivo:, coluna] seleciona a parte da coluna abaixo da linha pivô. Essa técnica será útil para tornar a eliminação mais estável e evitar divisões inválidas.

Separando \(A\), \(b\) e \(x\)

Nos métodos iterativos, como Jacobi e Gauss-Seidel, trabalharemos separando a matriz de coeficientes \(A\), o vetor \(b\) e o vetor de aproximações \(x\):

A = np.array([
    [8, -2, 5],
    [1, -5, 3],
    [-1, 3, -10]
], dtype=float)

b = np.array([43, 3, -28], dtype=float)
x = np.array([0, 0, 0], dtype=float)

Essa separação permite escrever as contas de forma próxima da notação matricial. Por exemplo, o produto \(Ax\) é calculado com A @ x, e o resíduo do sistema pode ser obtido por b - A @ x.

Para calcular uma nova aproximação de uma incógnita, isolamos uma variável por vez. Em Python, o somatório dos outros termos pode ser feito com sum:

n = len(b)
i = 0

soma = sum(A[i, j] * x[j] for j in range(n) if j != i)
x_novo = (b[i] - soma) / A[i, i]

A condição if j != i remove o termo da diagonal do somatório, pois ele fica isolado no denominador. Esse padrão aparece tanto no método de Jacobi quanto no método de Gauss-Seidel.

Cálculo do erro

Nos métodos iterativos, precisamos verificar se a aproximação já está boa o suficiente. Uma forma de fazer isso é calcular o erro usando a diferença entre \(b\) e \(Ax\):

Em Python com NumPy, essa conta pode ser escrita assim:

residuo = b - A @ x
erro_por_equacao = abs(residuo)
maior_erro = max(erro_por_equacao)

Outra forma comum é usar uma norma vetorial. Com NumPy, a norma infinito mede o maior erro absoluto entre as equações:

erro = np.linalg.norm(residuo, ord=np.inf)

Registrando aproximações iterativas

Nos métodos iterativos, é importante guardar a aproximação antiga e a nova para comparar a evolução. Use copy para evitar que as duas variáveis apontem para o mesmo vetor:

x_anterior = x.copy()
x_novo = x.copy()

# depois de calcular a nova iteração
erro_iteracao = np.linalg.norm(x_novo - x_anterior, ord=np.inf)

Exibindo resultados

Para acompanhar os cálculos, usaremos print. Quando for necessário misturar texto e valores, usaremos f-strings:

iteracao = 1
maior_erro = 0.025

print(f"Iteração {iteracao}: erro = {maior_erro}")

Com esses comandos, já é possível acompanhar as implementações das próximas aulas de sistemas lineares.

Atividade

Observe o código abaixo e calcule o resíduo b - A @ x para o chute inicial. Depois, use np.linalg.norm para medir o maior erro absoluto entre as equações:

import numpy as np

A = np.array([
    [4, 1, -1],
    [2, 7, 1],
    [1, -3, 12]
], dtype=float)

b = np.array([3, 19, 31], dtype=float)
x = np.array([0, 0, 0], dtype=float)