Calculo Numerico

O que realmente é isso quando a prática encontra a teoria

Você abre um livro de cálculo numérico e encontra definições bonitas sobre métodos iterativos, erros de truncamento e convergência assintótica. Depois vai para o código e descobre que o método que o livro dizia ser estável simplesmente explode em produção porque uma matriz mal condicionada não se importa com a sua teoria. É isso que o assunto é na vida real: a arte de encontrar respostas aproximadas quando a resposta exata não cabe no computador. O cálculo numérico existe porque a maioria dos problemas práticos não tem solução analítica fechada. Equações diferenciais parciais de engenharia, sistemas lineares com milhões de variáveis, integrais em domínios irregulares, otimização com restrições não-lineares. A matemática pura te diz que a solução existe. A numérica te diz como_approximá-la com um custo computacional que não quebre o orçamento da máquina.

Metodologias práticas de cálculo numérico

Vou começar pelo que todo mundo ignora: a análise de erro antes de qualquer método. A maioria das pessoas pula direto para implementar o algoritmo. Eu sempre começo perguntando qual tipo de erro domina o meu problema. Erro de truncamento vem da discretização. Erro de arredondamento vem da representação finita dos números. Erro de dados vem de medições imperfeitas. Se você não sabe qual deles é o gargalo, vai otimizar a coisa errada. Pegando um exemplo concreto que me custou duas semanas de debugging no passado: eu estava resolvendo um sistema linear Ax = b com o método de Gradientes Conjugados para uma matriz esparsa de elementos finitos. O residuo Caía lindamente nas primeiras cem iterações e depois travava num patamar de 10^-6 que não descia mais. O livro dizia que o GC converge em no máximo n iterações para matrizes simétricas definidas positivas. O problema era que a matriz tinha condição na casa de 10^14 devido a elementos com aspect ratio extremamente diverso na malha. Sem pré-condicionador, o método simplesmente não conseguia avançar. A solução foi implementar um pré-condicionador diagonal (Jacobi) que reduziu o número de condição efetivo para algo próximo de 10^3 e o residuo caiu para 10^-12 em quinze iterações. Lição: o método certo na presença errada é pior que o método errado.

Interpolacao polinomial é outro campo minado que os livros tratam de forma ingênua. Você tenta interpolar doze pontos com um único polinômio de grau onze e vê o comportamento de Runge acontecendo na sua cara: oscilações violentas nas bordas. A resposta pratica nao é usar Lagrange ou Newton com todos os pontos. É usar splines cúbicos naturais ou, se a funcao for suave mas o dominio grande, divisao em sub-intervalos com polinomios de baixo grau. Eu recentemente precisei interpolar dados experimentais de uma turbina a gas com sessenta pontos medidos em intervalos desiguais. Splines cúbicos com knots adaptativos baseados no criterio de Crooks deram um erro máximo de interpolacao da ordem de 10^-4, enquanto um polinômio global de Lagrange atingiu 10^2 nas extremidades. A diferença era pura estabilidade numerica. Quanto a integracao numerica, regra de Simpson parece suficiente ate voce tentar integrar uma funcao com descontinuidade ou pico agudo. Nesse caso, integracao adaptativa por quadratura de Gauss-Legendre com subdivisao recursiva é o padrão da indústria. Um pacote como o QUADPACK (disponível no SciPy como scipy.integrate.quad) aplica a regra de Gauss-Kronrod em sub-intervalos e refina onde a estimativa de erro local excede a tolerancia. Em testes comigo, integrar e^(-x^2) de 0 a infinito com tolerancia de 10^-10 levou menos de dois milissegundos. A mesma integral usando Simpson com passo fixo precisaria de aproximadamente dez milhões de avaliaçoes para atingir precisao similar.

Equações diferenciais ordinárias merecem uma discussão separada porque a escolha do método depende quase inteiramente da rigidez do sistema. Um sistema não-rígido responde bem a métodos de Runge-Kutta explícitos de quarta ordem (RK4). Um sistema rígido como os que aparecem em cinética química ou circuitos elétricos com constantes de tempo muito diferentes quebra métodos explícitos porque exige passos de tempo na casa de 10^-9 para estabilidade, mesmo que a solução variei suavemente em escala de segundos. Nesse caso, métodos implícitos como BDF (Backward Differentiation Formulas) ou backward Euler são obrigatórios. O custo por passo é maior porque exige resolução de sistemas lineares, mas o passo pode ser milhares de vezes maior. O scikits.odes e o solves do MATLAB resolvem isso automaticamente escolhendo entre métodos explícitos e implícitos baseado na estimativa da Jacobiana local. Aqui vai uma verdade contraintuitiva que pouca gente ensina: maior ordem do método não significa necessariamente melhor resultado final. Métodos de alta ordem como RK8 ou métodos de Adams-Moulton de ordem dez exigem muitas avaliações de função por passo e seu coeficiente de erro pode ser pior que um RK4 para problemas com ruído numérico significativo. Para a maioria dos problemas de engenharia, um RK4 com passo adaptativo e tolerância bem ajustada entrega precisão suficiente com menos sobrecarga. A ordem alta só compensa quando a função é extremamente suave e você precisa de precisão abaixo de 10^-12 consistentemente.

O outro ponto que as pessoas subestimam é a validade dos dados de entrada. Cálculo numérico opera sob o principio de que lixo entra e lixo sai. Eu já vi colegas tentarem refinar malhas e reduzir passos de tempo infinitamente sem perceber que os dados experimentais que alimentavam as condições de contorno tinham incerteza de cinco por cento. Nenhuma refinamento numerico vai eliminar esse erro. O trabalho correto nessa situacao é fazer análise de incerteza propagação usando métodos como Monte Carlo ou polinômios caóticos, nao simplesmente reduzir o passo numerico ate o computador travar. Para quem quer implementar isso de verdade, a hierarquia de ferramentas que eu recomendo baseada em experiencia pratica é a seguinte. Para prototipagem rapida e aprendizado, Python com NumPy, SciPy e Matplotlib cobre setenta por cento dos problemas do dia a dia. O SciPy implementa quase todos os métodos mencionados acima com interfaces consistentes. Para problemas maiores que exigem performance, substitua o NumPy por uma versao ligada a BLAS/LAPACK otimizada como a do Intel MKL ou OpenBLAS. Para elementos finidos, o deal.II ou o FEniCS Project sao mais adequados que programar do zero. Para EDOs rígidas em producao, o SUNDIALS (CVODE) é amplamente usado em simulacoes automotivas e aeroespaciais.

Se voce esta começando agora, o caminho é: dominar aritmética de ponto flutuante e seus limitacoes antes de tocar em qualquer metodo numerico. Entenda o que é machine epsilon, o que causa cancelamento catastrófico, porque a associatividade nao vale para operaçoes e como o formato IEEE 754 representa numeros no seu hardware. Esse conhecimento base evita que voce perca dias tentando debuggar um algoritmo que na verdade está matematicamente correto mas numericamente instavel para os dados que voce passou.