Circuito RLC: RKF45 com e sem Adaptação de Passo

Comparação numérica entre RKF45 de passo fixo e adaptativo, usando um circuito RLC série livre como caso-teste. A EDO tem solução fechada nos três regimes de amortecimento, então cada integrador pode ser comparado direto contra o resultado exato.

Formulação física

Circuito RLC série sem fonte, pela lei das tensões de Kirchhoff:

$$L\ddot{q} + R\dot{q} + \frac{q}{C} = 0$$

Como sistema de primeira ordem, $x_1=q$, $x_2=\dot q = i$:

$$ \begin{cases} \dot{x}_1 = x_2 \\ \dot{x}_2 = -\dfrac{R}{L}x_2 - \dfrac{1}{LC}x_1 \end{cases} $$

O discriminante de $Ls^2+Rs+\dfrac{1}{C}=0$ define três regimes, em função de

$$\zeta = \frac{R}{2}\sqrt{\frac{C}{L}}, \qquad \omega_0=\frac{1}{\sqrt{LC}}, \qquad \alpha=\frac{R}{2L}=\zeta\omega_0$$

Subamortecido ($\zeta<1$): $$q(t) = e^{-\alpha t}\left(A\cos\omega_d t + B\sin\omega_d t\right), \qquad \omega_d=\sqrt{\omega_0^2-\alpha^2}$$

$$A = q_0, \qquad B = \frac{i_0+\alpha q_0}{\omega_d}$$

Crítico ($\zeta=1$): $$q(t) = (A+Bt)e^{-\alpha t}$$

$$A = q_0, \qquad B = i_0+\alpha q_0$$

Superamortecido ($\zeta>1$): $$q(t) = A e^{s_1 t} + B e^{s_2 t}, \qquad s_{1,2}=-\alpha\pm\sqrt{\alpha^2-\omega_0^2}$$

$$A = \frac{i_0-s_2 q_0}{s_1-s_2}, \qquad B = q_0-A$$

$A,B$ vêm de $q(0)=q_0, \dot q(0)=i_0$. As três formas foram conferidas simbolicamente (resíduo nulo na EDO original) e contra uma integração de referência do scipy.integrate.solve_ivp (rtol=1e-13):

Regime	R (Ω)	$\max\lvert q_{\text{exata}}-q_{\text{scipy}}\rvert$
Subamortecido	12,649	$1{,}47\times10^{-19}$
Crítico	63,246	$6{,}90\times10^{-21}$
Superamortecido	316,228	$9{,}10\times10^{-20}$

Métodos numéricos

RKF45

Runge-Kutta encaixado de 6 estágios, tableau de Fehlberg (1969):

$$k_i = hf\left(t+c_i h, y+\textstyle\sum_{j < i}a_{ij}k_j\right), \qquad i=1,\dots,6$$

$$y_4 = y+\sum_{i=1}^{6}b_i^{(4)}k_i, \qquad y_5 = y+\sum_{i=1}^{6}b_i^{(5)}k_i$$

$c$	$k_1$	$k_2$	$k_3$	$k_4$	$k_5$	$k_6$
0
1/4	1/4
3/8	3/32	9/32
12/13	1932/2197	-7200/2197	7296/2197
1	439/216	-8	3680/513	-845/4104
1/2	-8/27	2	-3544/2565	1859/4104	-11/40
$b^{(4)}$	25/216	0	1408/2565	2197/4104	-1/5	0
$b^{(5)}$	16/135	0	6656/12825	28561/56430	-9/50	2/55

$y_5-y_4$ dá o erro local de graça, sem avaliações extras de $f$. As duas versões implementadas avançam com $y_5$; a diferença entre elas é só como $h$ é escolhido.

Conferi os coeficientes contra $y'=-y$, a razão de redução do erro a cada $h\to h/2$ converge pra 32 em $y_4$ e 64 em $y_5$, exatamente $2^5$ e $2^6$, como esperado pra ordem 4 e 5.

Passo fixo

$h$ constante.

Passo adaptativo

Erro por norma RMS escalada:

$$e_{\text{norm}} = \sqrt{\frac{1}{n}\sum_{k=1}^{n}\left(\frac{y_{5,k}-y_{4,k}}{\text{atol}+\text{rtol}\cdot\max(\lvert y_k\rvert,\lvert y_{5,k}\rvert)}\right)^2}$$

Aceita se $e_{\text{norm}}\le 1$ (mas o custo dos 6 rhs já foi pago, aceito ou não). Ajuste de passo:

$$h_{\text{novo}} = \text{clip}\left(0.9h\cdot e_{\text{norm}}^{-1/5},\quad 0.1h,\quad 5h\right)$$

Validação e ordem de convergência

Regime	$\zeta$	Ordem estimada	Erro mínimo observado
Subamortecido	0,20	5,15	$3{,}39\times10^{-18}$
Crítico	1,00	5,20	$2{,}34\times10^{-18}$
Superamortecido	5,00	5,05	$2{,}38\times10^{-23}$

O superamortecido satura em ruído de arredondamento antes de esgotar a faixa de $h$ testada, por isso o erro mínimo ali é bem menor que nos outros.

Casos-base nos três regimes

Regime	Método	Passos	Avaliações	Erro final
Subamortecido	Fixo	784	4704	$3{,}39\times10^{-18}$
Subamortecido	Adaptativo	581 (22 rej.)	3618	$1{,}72\times10^{-16}$
Crítico	Fixo	160	960	$2{,}34\times10^{-18}$
Crítico	Adaptativo	151 (2 rej.)	918	$4{,}41\times10^{-17}$
Superamortecido	Fixo	15678	94068	$3{,}95\times10^{-22}$
Superamortecido	Adaptativo	319 (4 rej.)	1938	$5{,}61\times10^{-19}$

No superamortecido o fixo gasta 48,5× mais avaliações que o adaptativo (94068 vs. 1938), o caso onde a adaptação compensa mais.

Variação paramétrica

Variação de R

$h$ e tolerância fixos no valor tunado no ponto crítico; $\zeta$ de 0,05 a 30 variando só $R$:

$\zeta$	R (Ω)	Erro fixo	Erro adapt.	Aval. fixo	Aval. adapt.	Razão custo
0,050	3,16	$1{,}49\times10^{-17}$	$8{,}86\times10^{-16}$	19200	14448	1,33
0,180	11,37	$3{,}44\times10^{-18}$	$1{,}92\times10^{-16}$	5340	4032	1,32
0,834	52,77	$1{,}21\times10^{-18}$	$2{,}09\times10^{-17}$	1152	978	1,18 (mín.)
1,798	113,70	$2{,}16\times10^{-21}$	$4{,}87\times10^{-17}$	3162	1230	2,57
3,874	244,99	$7{,}75\times10^{-25}$	$4{,}17\times10^{-18}$	7314	1434	5,10 (máx.)
8,346	527,86	$3{,}79\times10^{-23}$	$2{,}68\times10^{-19}$	15966	4194	3,81
30,000	1897,37	$5{,}61\times10^{-22}$	$7{,}86\times10^{-20}$	57582	47430	1,21

Razão de custo entre 1,18 e 5,10: o adaptativo nunca perde por muito, e passa de 5× mais barato perto de $\zeta\approx 3{,}9$.

Variação de L e de C

$\alpha=R/2L$ depende de L mas não de C; aumentar R ou diminuir L empurra $\alpha$ sem limite, o que pode instabilizar o fixo; $\omega_0$ depende dos dois igual.

Variando L, o fixo instabilizou em 2 dos 24 pontos (L = $1{,}00\times10^{-5}$ H e $2{,}23\times10^{-5}$ H; nesse segundo o "erro" calculado chega a $2{,}44\times10^{131}$, finito em ponto flutuante mas sem sentido físico). Variando C, nenhum dos 24 pontos instabilizou, consistente com $\alpha$ não depender de C.

Parâmetro	Erro fixo	Erro adapt.	Aval. fixo	Aval. adapt.
L = $1{,}00\times10^{-5}$ H	instável	$3{,}02\times10^{-20}$	1914	5988
L = $2{,}23\times10^{-5}$ H	instável	$7{,}72\times10^{-21}$	1908	3108
L = $1{,}22\times10^{-3}$ H (≈ζ=1)	$1{,}26\times10^{-18}$	$3{,}69\times10^{-17}$	1170	936
L = $1{,}00\times10^{-1}$ H	$8{,}58\times10^{-21}$	$6{,}64\times10^{-16}$	96000	4872
C = $1{,}00\times10^{-8}$ F	$3{,}10\times10^{-13}$	$4{,}25\times10^{-18}$	960	9780
C = $8{,}19\times10^{-7}$ F (≈ζ=1)	$3{,}39\times10^{-18}$	$2{,}97\times10^{-17}$	960	966
C = $1{,}00\times10^{-4}$ F	$1{,}04\times10^{-21}$	$3{,}53\times10^{-18}$	191520	5568

Mapa de calor R×L

	Fixo	Adaptativo
Erro mínimo	$0{,}00\times10^{0}$	n/d
Erro mediano	$3{,}72\times10^{-19}$	n/d
Erro máximo (onde estável)	$9{,}37\times10^{-14}$	n/d
Avaliações mínimas	n/d	186
Avaliações medianas	n/d	933
Avaliações máximas	n/d	23976

28 dos 196 pontos (14,3%) divergiram no fixo, todos do lado superamortecido, como esperado.

Análise computacional

Diagrama de trabalho-precisão

Caso suave ($\zeta=0{,}2$) e rígido ($\zeta=5{,}0$). Com custo parecido (~3000 avaliações) no caso suave, fixo e adaptativo ficam próximos: $1{,}67\times10^{-17}$ (3006 avaliações) contra $2{,}92\times10^{-17}$ (3306). No rígido a diferença é dramática: com $h$ grande (192 avaliações) o fixo diverge, erro chegando a $3{,}4\times10^{141}$, só estabilizando a partir de ≈1896 avaliações; o adaptativo não diverge em nenhum ponto testado, e com 1296 avaliações já está em $1{,}15\times10^{-9}$.

Taxa de rejeição vs. tolerância

Entre 1,04% (tolerância $10^{-14}$) e 18,75% ($1{,}27\times10^{-4}$).

Tempo de CPU real

Tunado pra erro final $\approx 10^{-10}$, 20 repetições por caso:

Caso	Fixo (ms)	Adaptativo (ms)
Suave	$1{,}78\pm0{,}05$	$10{,}26\pm0{,}41$
Rígido	$12{,}61\pm0{,}22$	$17{,}59\pm0{,}30$

Em Python puro o adaptativo perde nos dois casos em tempo real, mesmo gastando menos avaliações de rhs no rígido; o overhead do laço de controle domina. Numa implementação vetorizada essa vantagem provavelmente se confirmaria.

Conservação de energia

$$E(t) = \frac{1}{2}Li(t)^2 + \frac{q(t)^2}{2C}$$

t (s)	E numérica (J)	E exata (J)	Erro relativo
0	$5{,}000000\times10^{-7}$	$5{,}000000\times10^{-7}$	0
$1{,}05\times10^{-4}$	$1{,}386294\times10^{-7}$	$1{,}386118\times10^{-7}$	$3{,}51\times10^{-5}$ (máx.)
$4{,}80\times10^{-4}$	$9{,}727738\times10^{-10}$	$9{,}723517\times10^{-10}$	$8{,}44\times10^{-7}$

Erro relativo máximo: $3{,}512\times10^{-5}$, na mesma ordem do erro de $q(t)$ isolado.

Estrutura do notebook (rlc_rkf45.ipynb)

Definições → validação → funções auxiliares → casos-base (um bloco por regime) → variação paramétrica (R, L, C, mapa de calor) → análise computacional (convergência, trabalho-precisão, rejeição, CPU, energia). Toda comparação de dados tem tabela impressa no notebook, além do gráfico.

Métrica de custo

Avaliações de rhs, contadas por um contador global incrementado a cada chamada; cada passo do RKF45 custa 6, aceito ou rejeitado. Independente de hardware, mas não necessariamente de tempo real: ver seção de CPU acima.

Referências

• Fehlberg, E. (1969). Low-order classical Runge-Kutta formulas with stepsize control. NASA TR R-315.
• Boyce, W. E., DiPrima, R. C., Meade, D. B. (2017). Elementary Differential Equations and Boundary Value Problems (11th ed.). Wiley.
• Azarani, N. A., Md Sudin, A. (2025). Fehlberg Runge-Kutta Method for Solving System of Ordinary Differential Equations. Enhanced Knowledge in Sciences and Technology, 5(2), 123–132. https://doi.org/10.30880/ekst.2025.05.02.015