Resolvendo equações lineares na prática
A equacao linear é uma das ferramentas mais usadas em computação numérica, mas também uma das mais mal compreendidas por quem está começando. A forma padrão ax + b = 0 parece simples demais, e é exatamente essa aparência que prega peças. O problema não está em resolver um exemplo do livro. O problema aparece quando você tenta implementar isso de verdade. Vamos começar pelo básico mesmo assim, porque muita gente pula essa parte e depois se perde. Uma equação linear é qualquer expressão que pode ser escrita na forma ax + b = 0, onde a e b são constantes conhecidas e x é a incógnita. A solução é simplesmente x = -b/a, desde que a seja diferente de zero. Quando a é zero, você tem dois casos: se b também for zero, qualquer valor de x satisfaz a equação; se b for diferente de zero, não existe solução alguma.
Em sistemas de equações lineares, a coisa muda de figura. M matrizes e n variáveis. Você pode ter uma solução única, infinitas soluções ou nenhuma solução. Isso depende da posto da matriz dos coeficientes e da posto da matriz ampliada. Se os postos forem iguais ao número de variáveis, solução única. Se o posto for menor, depende do sistema ser compatível ou não.
Por que sua implementação de equacao linear pode estar errada
Eu passei duas semanas caçando um bug em um sistema de balanceamento estrutural em 2019. O modelo tinha 48 equações lineares com 48 incógnitas. Os resultados não batiam com a planilha de referência em menos de 0,1% para a maioria dos casos, mas havia três nós onde o erro saltava para 15%. Descobri que dois dos elementos da matriz de rigidez tinham magnitudes extremamente diferentes — uma na ordem de 10^8 e outra de 10^-2 — porque um nó conectava uma viga de aço a uma mola de calibração muito flexível. O solver direto (eliminação gaussiana pura) estava acumulando erro de arredondamento que destruía a solução naqueles graus de liberdade específicos. A solução foi pré-condicionar a matriz antes de resolver. Multipliquei as linhas correspondentes aos elementos de mola por um fator de 10^10 para colocar todas as grandezas na mesma escala. O tempo de resolução aumentou cerca de 30% porque o pré-condicionador adicionou uma etapa extra, mas a precisão voltou ao esperado. Isso é algo que nenhum tutorial básico ensina.
👉 Clique no botão abaixo para saber mais sobre o assunto!
O ponto é que resolver uma equação linear no papel e resolver uma equação linear no computador são coisas fundamentalmente diferentes. No papel você faz contas exatas. No computador você lida com números de ponto flutuante, e o erro se propaga de maneiras que você nem sempre prevê. Um insight que vejo quase ninguém mencionar: o método dos fatores de Gauss (eliminação gaussiana com pivotação parcial) é o padrão da indústria para sistemas pequenos e médios, digamos até 1000 equações. Para sistemas maiores, métodos iterativos como Gauss-Seidel ou, melhor ainda, o gradiente conjugado (quando a matriz é simétrica e definida positiva) são mais eficientes em memória e tempo. O gradiente conjugado converge em no máximo n iterações para uma matriz n x n em aritmética exata, mas na prática com ponto flutuante você para quando o resíduo fica abaixo de uma tolerância, geralmente 10^-6 ou mais restritivo dependendo da aplicação.
Outro detalhe que passa despercebido: pivotaçao. Sempre faça pivotação parcial. Resolver sem pivotação é pedir para o erro numérico explodir. Um exemplo clássico é a matriz [[1e-20, 1], [1, 1]]. Se você eliminar sem trocar linhas, divide por 1e-20 e o segundo passo introduz um erro catastrófico. Com pivotação, troca as linhas e resolve normalmente. Se você está implementando isso do zero, há bibliotecas prontas que já fazem tudo certo. A mais citada é o LAPACK, especificamente os rotinas dgesv para sistemas gerais e dpotrf/dpotrs para matrizes simétricas definidas positivas. Numérica e robustez estão décadas à frente do que você vai escrever num fine de tarde. A menos que você esteja desenvolvendo um solver novo para pesquisa, usar uma biblioteca estabelecida é a escolha certa. Leva cerca de 10 linhas de código dgesv contra centenas para fazer elimination gaussiana correta com pivotação e back-substitution do jeito certo.
Uma limitação importante que precisa ser dita claramente: equações lineares não modelam tudo. Se o seu problema tem não-linearidades — e a maioria dos problemas reais tem — linearizar em torno de um ponto de operação e resolver o sistema linear resultante é uma aproximação. O método de Newton-Raphson faz exatamente isso, resolvendo um sistema linear a cada iteração. Mas se a linearização estiver longe da solução real, o método pode divergir. Não há garantia de convergência sem condições específicas no domínio do problema. Outra limitação prática: custo computacional. Fatores de Gauss exigem O(n^3) operações e O(n^2) de memória. Para n = 10.000, isso é cerca de 10^12 operações e 800 GB de memória em ponto duplo. Inviável na maioria dos hardwares comuns. Para essas escalas, você precisa de métodos iterativos ou aproveitar a esparsidade da matriz com formatos como CSR (Compressed Sparse Row) e solvers especializados como o SuiteSparse.
Se o seu sistema tem milhões de equações e a matriz é esparsa, considere o PETSc ou o Trilinos. São frameworks que abstraem a complexidade de pré-condicionadores, métodos iterativos e paralelização. Configurar tudo do zero leva meses. Usar PETSc com um pré-condicionador ILU(0) e o solver GMRES resolve sistemas de 500 mil equações em minutos em uma máquina com 16 núcleos, dependendo da condicionamento da matriz. Na prática, o fluxo de trabalho mais comum é: montar a matriz dos coeficientes, verificar se ela é esparsa ou densa, escolher o solver adequado (direto para sistemas pequenos e densos, iterativo para grandes e esparsos), aplicar pré-condicionamento se necessário, e validar a solução substituindo-a de volta na equação original para calcular o resíduo. Um resíduo relativo menor que 10^-10 geralmente indica que a solução é confiável para a maioria das aplicações de engenharia.
