Cálculo Numérico
Universidade Federal do Pará
Uma EDO de ordem \(n\) converte-se em um sistema de \(n\) equações de 1ª ordem. Para \(m\ddot{x} + c\dot{x} + kx = F(t)\), define-se \(x_1 = x,\ x_2 = \dot{x}\):
\[\begin{cases} \dot{x}_1 = x_2 \\ \dot{x}_2 = \dfrac{1}{m}\bigl(F - c\,x_2 - k\,x_1\bigr) \end{cases}\]
Os mesmos métodos (Euler, RK4, solve_ivp) aplicam-se ao vetor de estado \(\mathbf{y} = [x_1, x_2]^\top\).
import numpy as np
from scipy.integrate import solve_ivp
m, c, k = 1.0, 0.5, 20.0 # massa-mola-amortecedor
def sistema(t, y):
x, v = y
return [v, (-c*v - k*x) / m]
sol = solve_ivp(sistema, [0, 10], [1.0, 0.0], method='RK45',
t_eval=np.linspace(0, 10, 500), rtol=1e-8)
x, v = sol.y # deslocamento e velocidadeA solução vetorial fornece todos os estados simultaneamente.
Tip
📓 Notebook complementar — mesmo sistema, com força de excitação, comparando Euler explícito e solve_ivp, e relacionando os polos (\(m s^2+cs+k\)) à estabilidade.
Aplica-se a equação-teste \(y' = \lambda y\), com \(\mathrm{Re}(\lambda) < 0\). Para o Euler explícito:
\[y_{k+1} = (1 + h\lambda)\,y_k \;\Rightarrow\; \text{estável se } |1 + h\lambda| < 1\]
Para \(\lambda\) real negativo, exige-se:
\[h < \frac{2}{|\lambda|}\]
Important
Passo grande demais torna a solução numérica divergente, ainda que a solução exata decaia. Métodos explícitos têm região de estabilidade limitada.
Sistemas com escalas de tempo muito distintas (autovalores de magnitudes díspares) forçam passo minúsculo nos métodos explícitos.
Sintoma: RK45 fica lentíssimo ou diverge, mesmo com solução suave.
Solução: métodos implícitos com região de estabilidade ampla — Radau, BDF.
| Situação | Método em solve_ivp |
|---|---|
| Não rígido, uso geral | RK45 (padrão) |
| Precisão alta, não rígido | DOP853 |
| Rígido | Radau, BDF |
| Solução oscilatória suave | LSODA (alterna automaticamente) |
Tip
Na dúvida sobre rigidez, LSODA detecta e troca entre estratégias explícita e implícita automaticamente.
A análise de estabilidade anterior é linear — vale para \(y'=\lambda y\) ou para a linearização de um sistema em torno de um ponto de equilíbrio. Sistemas não-lineares podem exibir comportamento que a teoria linear não prevê: o caos determinístico.
O exemplo clássico é o sistema de Lorenz, modelo simplificado de convecção atmosférica:
\[\begin{cases} \dot{x} = \sigma(y-x) \\ \dot{y} = x(\rho - z) - y \\ \dot{z} = xy - \beta z \end{cases}\]
Sem nenhum termo aleatório, o sistema é totalmente determinístico — mas duas trajetórias com condições iniciais que diferem por \(10^{-6}\) divergem exponencialmente após poucas dezenas de unidades de tempo, tornando a previsão de longo prazo impraticável.
Important
Sensibilidade às condições iniciais não é instabilidade numérica: mesmo a solução exata do sistema tem esse comportamento. Ainda assim, a trajetória permanece confinada a uma região limitada do espaço de estados — o atrator estranho.
📓 Notebook complementar — integra duas trajetórias com condições iniciais quase idênticas e traça o atrator de Lorenz (retrato de fase e vista 3D).
solve_ivp.RK45 e com BDF, comparando número de passos e tempo de execução.