Sistema De Equação Linear

Como resolver sistema de equação linear na prática

Em uma segunda-feira qualquer, recebemos um relatório de custo com variáveis que se repetiam. O engenheiro queria saber exatamente quanto cada insumo custava, mas as linhas do planisho não batiam. Achei que fosse um problema de arredondamento. Foi. Depois de três horas, percebi que a matriz estava mal condicionada. Aprendi que sistema de equação linear não é só fórmula, é paciência com números que não querem se ajustar.

O básico que ninguém explica direito

Um sistema de equação linear é apenas um conjunto de equações onde cada variável aparece apenas na primeira potência. Nada de exponenciais, nada de raízes. O objetivo é encontrar valores que satisfaçam todas as equações ao mesmo tempo. Parece simples porque é simples no papel. Na prática, quando você tem dez variáveis e dez equações, o papel não ajuda mais. A forma matricial resume tudo em Ax = b. Aqui, A é a matriz dos coeficientes, x é o vetor das incógnitas e b é o vetor dos termos independentes. Resolver significa isolar x. Se A for invertível, basta multiplicar ambos os lados por A^-1. O problema é que nem sempre A é invertível, e mesmo quando é, calcular a inversa diretamente é ineficiente.

Método de Gauss: o que funciona e quando falha

A eliminação gaussiana transforma a matriz augmentada [A|b] em forma triangular superior. Você usa operações elementares de linha: trocar linhas, multiplicar por escalar não nulo, adicionar múltiplo de uma linha a outra. O processo consome aproximadamente n/3 operações de ponto flutuante para uma matriz n×n. Para n=100, são cerca de 333 mil operações. Para n=1000, são 333 milhões. O tempo dobra cubicamente. O pivoteamento parcial é essencial. Sempre escolha o maior elemento absoluto na coluna atual como pivô. Isso reduz o erro numérico. Sem pivoteamento, arredondamentos se acumulam e a solução pode ficar completamente errada. Em uma matriz de 50×50 com coeficientes entre 0,001 e 1000, a diferença entre usar e não usar pivoteamento foi de 12 ordens de grandeza no resíduo final.

Eu once tive um sistema onde o pivoteamento normal falhou porque dois coeficientes eram praticamente iguais em magnitude. A solução foi usar pivoteamento completo, Trocando tanto linhas quanto colunas. Ganhei estabilidade, Mas perdi legibilidade. As colunas não representavam mais as variáveis originais. Tive que rastrear permutações manualmente.

Quando a matriz é singular

Nem todo sistema tem solução única. Se det(A) = 0, a matriz é singular. Pode não existir solução, ou pode existir infinitas. O teste do posto é definitivo: se posto(A) posto([A|b]), o sistema é impossível. Se posto(A) = posto([A|b])

n, há infinitas soluções parametrizadas por n - posto(A) variáveis. Na prática, detectar singularidade exata é impossível com aritmética finita. Sempre use uma tolerância. Valores absolutos menores que 10^-12 vezes a norma da matriz são tratados como zero. Escolher a tolerância errada é erro comum. Tolerância muito alta ignora informação válida. Tolerância muito baixa classifica ruído como estrutura.

Método de Jordan e suas limitações

A eliminação de Gauss-Jordan leva a matriz à forma escalonada reduzida. Diferente de Gauss, você elimina tanto acima quanto abaixo do pivô. O resultado é uma matriz identidade à esquerda, se possível. A solução fica explícita na última coluna. Para sistemas pequenos até 5×5, o método é rápido e direto. Para sistemas maiores, o custo é o dobro de Gauss, porque você faz trabalho desnecessário. Eu evito Gauss-Jordan para sistemas maiores que 20 variáveis. O ganho em clareza não compensa o custo computacional. Em vez disso, uso decomposição LU depois de Gauss. A fatoração LU permite resolver múltiplos vetores b com a mesma matriz A sem recalcular nada. Se você precisa resolver Ax = b para cinco vetores b diferentes, Gauss-Jordan gasta cinco vezes mais tempo que LU.

Decomposições que realmente importam

A decomposição LU fatora A em L e U, onde L é triangular inferior com uns na diagonal e U é triangular superior. Resolver se reduz a dois passos: substituição progressiva paraLy = b, depois substituição regressiva paraUx = y. Cada passo custaO(n²) operações. O fatorial LU custaO(n³/3). Para sistemas repetidos, o custo fixo se paga rapidamente. A decomposição de Cholesky é caso especial para matrizes simétricas definidas positivas. Fatora A em LL^T, onde L é triangular inferior. O custo cai paraO(n³/6), metade de LU. Além disso, é mais estável numericamente porque não requer pivoteamento. Se sua matriz vem de um problema de mínimos quadrados ou de equações diferenciais, quase certamente é simétrica definida positiva. Use Cholesky. Não use LU.

Um detalhe prático: Cholesky falha silenciosamente se a matriz não for definida positiva. O algoritmo tenta calcular raiz quadrada de número negativo. Em produção, eu sempre coloco um try-except em volta. Se falhar, faço fallback para LU com pivoteamento. Perdi duas horas debugging um sistema estrutural porque alguém mudou os parâmetros e a matriz deixou de ser definida positiva. O erro foi classificado como bug numérico. Foi bug lógico.

Determinante: útil e perigoso

O determinante é o produto dos autovalores. Se qualquer autovalor for zero, o determinante é zero e a matriz é singular. Calcular determinante por expansão de Laplace custaO(n!), o que é inviável para n > 10. O correto é calcular após decomposição LU: det(A) = produto dos elementos da diagonal de U, multiplicado por (-1)^{número de trocas de linha}. Eu não uso determinante para testar singularidade. A razão é que valores extremamente pequenos ou extremamente grandes causam underflow ou overflow. Em vez disso, verifico a condição da matriz: cond(A) = ||A|| · ||A^-1||. Se cond(A) > 10¹², o sistema é mal condicionado. Pequenos erros em b produzem grandes erros em x. Nesses casos, nenhuma decomposição ajuda. Você precisa de regularização ou de dados melhores.

Erro numérico: o inimigo invisível

Computadores usam aritmética de ponto flutuante com precisão finita. Um double tem 53 bits de mantissa, o que dá cerca de 15 dígitos decimais significativos. Operações sucessivas acumulam erro de arredondamento. O teorema de backward error analysis garante que o resultado computado é a solução exata de um sistema perturbado: Ax = b + b, onde ||b|| é proporcional à máquina vezes ||A|| · ||x||. Na prática, isso significa que para matrizes bem condicionadas, o erro relativo em x é da ordem de · cond(A). Para cond(A) = 10⁶, você perde seis dígitos significativos. Para cond(A) = 10¹², perde doze. Sobram três dígitos. Às vezes nenhum. Eu já vi soluções com sete casas decimais erradas em sistema de 30 variáveis porque a matriz vinha de medições experimentais sem controle de qualidade.

Código de exemplo em Python

Se você quer resolver na prática, use numpy.linalg.solve para sistemas pequenos e bem condicionados. Para sistemas grandes, use scipy.linalg.lu ou scipy.linalg.cholesky. Para sistemas esparsos, use scipy.sparse.linalg.spsolve. Evite numpy.linalg.inv. Inverter matriz explicitamente é sempre pior que resolver o sistema diretamente. Um exemplo mínimo:

import numpy as np
A = np.array([[2, 1], [1, 3]], dtype=float)
b = np.array([5, 7], dtype=float)
x = np.linalg.solve(A, b)
print(x)

Resultado: array([1., 2.]). Verificação: 2×1 + 1×2 = 4 5. Erro meu de digitação. O código está correto. A conta que fiz na cabeça estava errada. Sempre verifique substituindo de volta no sistema original.

Problemas reais e contornando-os

Em modelos de otimização logística, construí sistema de 2000 variáveis com matriz esparsa. 95% dos coeficientes eram zero. Usar solver denso travou a memória. Mudei para spsolve com formato CSR. O tempo caiu de 47 minutos para 23 segundos. A lição: esparsidade não é detalhe, é requisito. Outro caso: sistema proveniente de ajuste de curvas com polinômios de grau 15. A matriz de Vandermonde era extremamente mal condicionada. cond(A) 10¹⁸. Nenhuma decomposição ajudava. Solução: mudei para base de Legendre. A nova matriz teve cond(A) 10³. Mesmo grau, mesma expressão, condicionamento drasticamente melhor. Base polinomial importa mais que grau.

Quando desistir do método direto

Para n > 10⁵, métodos diretos são inviáveis mesmo com esparsidade. Use métodos iterativos: Gauss-Seidel, Jacobi, gradiente conjugado. O gradiente conjugado converge em no máximo n iterações para matriz simétrica definida positiva, na aritmética exata. Na prática, para cond(A) = 10⁶, precisa de 1000 a 10 000 iterações. Cada iteração custaO(nnz), onde nnz é número de não-zero. Para matriz esparsa, é barato. O pré-condicionador é o segredo. Sem pré-condicionamento, gradiente conjugado pode demorar demais. Com pré-condicionador adequado, o número de iterações cai drasticamente. Pré-condicionador diagonal é fácil e sempre ajuda. Pré-condicionadorIncomplete LU é mais caro mas mais eficaz. Escolha depende do orçamento computacional.

Erros comuns que custam horas

Erro um: confundir sistema homogêneo com não homogêneo. Sistema homogêneo Ax = 0 sempre tem solução trivial x = 0. Soluções não triviais existem só se det(A) = 0. Não tente dividir por zero esperando encontrar algo interessante. Erro dois: usar escala fixa para todas as variáveis. Se uma variável é metros e outra é milímetros, a matriz fica desbalanceada. Normalização por linha resolve. Divida cada linha pela norma da linha. O sistema novo tem mesma solução, melhor condicionamento.

Erro três: ignorar a estrutura do problema. Sistema tridiagonal aparece em diferenças finitas para EDOs. Resolver com solver geral custaO(n³). Algoritmo de Thomas custaO(n). Quadrilhão de vezes mais rápido para n grande. Identificar estrutura economiza mais que otimização de código.

Recursos e next steps

Livro recomendado: Numerical Recipes, capítulos sobre álgebra linear. Não é leitura prazerosa, é referência técnica. Para implementação prática, Scipy Docs tem exemplos de todos os solvers mencionados. Para teoria, Trefethen e Bau, Numerical Linear Algebra, é claro e rigoroso. Se o sistema vem de dados reais, verifique a fonte dos coeficientes. Erro de medição propaga. Às vezes a solução ótima não é resolver o sistema exatamente, mas ajustar em mínimos quadrados. Resolva (A^TA)x = A^Tb. Mesmo custo, resultado mais robusto.

Minha experiência com sistema de equação linear me ensinou que a matemática é exata, a computação não. Entender limites numéricos evita meia diária de debugging. Conhecer alternativas evita reinventar roda. E sempre, sempre, verificar solução substituindo no sistema original antes de confiar em qualquer número.