Sistema Linear

Como resolver sistemas lineares na prática

Quando você começa a mexer com álgebra linear, a primeira coisa que encontra é um sistema de equações como 2x + 3y = 8 e x - y = 1. Parece fácil até você ter 15 variáveis e 12 equações e precisar de um método que não envolva substituição manual. Um sistema linear é basicamente um conjunto de equações onde todas as incógnitas aparecem apenas com grau 1. Nada de x², nada de xy. A forma padrão é Ax = b, onde A é uma matriz de coeficientes, x é o vetor de variáveis e b é o vetor dos termos independentes. Pronto. É isso. A representação matricial não é apenas notação elegante. Ela permite usar algoritmos estruturados em vez de manipular equações uma a uma. O problema é que a maioria dos tutoriais para quando chega na parte de como você realmente resolve na prática.

Métodos de resolução

O método mais básico é a eliminação gaussiana. Você transforma a matriz ampliada [A|b] em forma escalonada usando operações elementares por linhas e depois faz retrosubstição. Funciona para qualquer sistema quadrado com determinante diferente de zero. Na mão, é viável para 3x3 ou no máximo 4x4. Acima disso, o tempo de cálculo cresce rápido e o risco de erro aritmético também. Para matrizes maiores, a decomposição LU é o padrão da indústria. Você decompõe A em L (matriz triangular inferior) e U (triangular superior) e resolve dois sistemas triangulares em vez de um geral. Isso economiza tempo quando você precisa resolver o mesmo sistema A com vários vetores b diferentes, porque a decomposição é feita uma única vez. Já se o sistema é overdetermined — mais equações do que incógnitas — a eliminação direta não funciona da mesma forma. Nesse caso, recorre-se ao método dos mínimos quadrados, resolvendo A^T A x = A^T b. O problema aqui é que A^T A dobra o número de condição, então se A já era mal condicionada, o resultado fica instável.

Um caso real que quase me custou horas

Há alguns anosWorking com um sistema linear de 47 equações e 47 variáveis oriundo de uma malha estrutural. A matriz era esparsa, mas quando converti para formato denso e mandei resolver com uma rotina de eliminação gaussiana comum, o solver começou a gerar zeros na diagonal que não deveriam existir. O sistema tinha solução, mas a condição numérica estava tão ruim que o pivoteamento parcial padrão não bastava. A solução foi ativar o pivoteamento completo, que troca tanto linhas quanto colunas para garantir o maior pivô possível em cada etapa. Em vez de usar o solver genérico, parti para uma fatoração LU com pivoteamento parcial revisado e verifiquei a norma da matriz invertida para estimar o número de condição. Descobri que estava em torno de 10^8, o que explicava a perda de precisão. Rodar com precisão dupla e usar um pré-condicionador diagonal antes da eliminação estabilizou o resultado.

Pegadinhas que ninguém avisa

O número de condição é o dado mais importante e o mais ignorado. Um sistema pode ser perfeitamente determinável e ainda assim dar soluções completamente erradas se o número de condição for alto. A regra prática: se o número de condição estiver acima de 10^8 em precisão dupla, desconfie do resultado. Não há mágica que resolva isso só com mais iterações. Outro erro comum é achar que todo sistema tem solução única. Sistema homogêneo Ax = 0 sempre tem a solução trivial x = 0, mas se o determinante for zero, existem infinitas soluções formando o núcleo da matriz. Sistema incompatível, com mais restrições do que graus de liberdade, simplesmente não tem solução e o solver vai retornar algo sem sentido se você não verificar o rank antes.

Quando usar cada abordagem

Sistemas pequenos e densos até 20x20: eliminação gaussiana direta é suficiente e rápida. Sistemas médios 20x20 a 500x500: decomposição LU com bibliotecas como LAPACK ou scipy.linalg.solve. Sistemas grandes e esparsos acima de 500x500: use métodos iterativos como GMRES ou conjugate gradient, preferencialmente com um pré-condicionador. Para overdetermined: mínimos quadrados com QR decomposition em vez de normal equations, porque o QR preserva melhor a precisão numérica. Se você precisa implementar isso do zero, a referência padrão é o código de Golub e Van Loan. As rotinas dgetrf e dgesv do LAPACK resolvem sistemas densos com pivoteamento parcial. Para esparsos, o SuiteSparse ou o PARDISO são opções maduras.

Sobre sistema linear

O conceito em si é simples. A dificuldade está nas escolhas numéricas que você faz ao implementá-lo. Matrizes mal condicionadas, esparsidade mal explorada, e a tentação de usar a fórmula mais óbvia em vez da mais estável são os três problemas que vejo repetidamente. A maioria dos erros em simulações estruturais, ajuste de curvas e modelos de fluxo vêm de um desses três pontos, não de um erro conceitual no sistema linear em si.