Aviso: A notação usada nesta página não é complicada, entretanto você pode não compreender se não entender a notação, portanto recomento que vejas antes Equações Diferenciais Ordinárias onde explico a notação.
----------------------------------------
1.0 - O método de Runge-Kutta de 4ª ordem:
----------------------------------------
O Método de Runge-Kutta de 4ª ordem, ou simplesmente RK4, é um método de cálculo numérico que permite colher resultados com maior precisão do que se usássemos o método de Euler ou, se quisermos a mesma precisão, com menos tempo do que usando Euler.
Para o PVI: \[ \left \{ \begin{align*} \frac{d}{dx}f(x) &= F(x, \; f(x) ) \\ f(x_0) &= C \end{align*}\right . \; \; \; \; (1) \] o método de Euler de primeira ordem pode ser expresso por: \[ f(x + dx) = f(x) + I dx \] Onde $I$ é uma inclinação. No método de Euler essa inclinação é a derivada $\frac{d}{dx}f(x)$ no ponto de abscissa $x$, no RK4 essa inclinação é mais "detalhada", expressa por: \[ I=\frac{1}{6} \left ( k_1 + 2k_2 + 2k_3 + k_4 \right ) \; \; \; \; (2) \] onde: \[ \left \{ \begin{align*} k_1 &= F\left ( x, \; \; f(x) \right )\\ k_2 &= F\left ( x + \frac{dx}{2}, \; \; f(x) + k_1 \, \frac{dx}{2} \right ) \\ k_3 &= F\left ( x + \frac{dx}{2}, \; \; f(x) + k_2 \, \frac{dx}{2} \right ) \\ k_4 &= F\left ( x + dx, \; \; f(x) + k_3 \, dx \right ) \end{align*} \right . \] ou seja:
· $k_1$ é a inclinação de $f$ no início do intervalo $[x, x + dx]$, exatamente a mesma de Euler;
· $k_2$ é a inclinação de $f$ na metade do intervalo $[x, x + dx]$ estimada através do método de Euler usando $k_1$;
· $k_3$ é a inclinação de $f$ na metade do intervalo $[x, x + dx]$ estimada através do método de Euler usando $k_2$;
· $k_4$ é a inclinação de $f$ no fim do intervalo $[x, x + dx]$ estimada através do método de Euler usando $k_3$.
----------------------------------------
2.0 - Exemplo (1ª ordem):
----------------------------------------
Seja o PVI: \[ \left \{ \begin{align*} \frac{d}{dx}f(x) &= x \; f(x) \\ f(x_0) &= C \end{align*}\right . \] Para estimar o valor de $f(x_0 + dx)$ basta aplicar RK4: \[ f(x_0 + dx) = f(x_0) + \frac{1}{6} \left ( k_1 + 2k_2 + 2k_3 + k_4 \right ) dx \] com: \[ \left \{ \begin{align*} k_1 &= F(x_0, f(x_0)) = x_0 \; C \\ k_2 &= (x_0 + \frac{dx}{2}) \; (C + k_1 \frac{dx}{2}) \\ k_3 &= (x_0 + \frac{dx}{2}) \; (C + k_2 \frac{dx}{2}) \\ k_4 &= (x_0 + dx) \; (C + k_3 dx) \end{align*}\right . \] Finalmente, as equações do método RK4 para o problema, que seriam colocadas dentro de um loop em um programa, serão: \[ \left \{ \begin{align*} k_1 &= x \; f(x) \\ k_2 &= (x + \frac{dx}{2}) \; (f(x) + k_1 \frac{dx}{2}) \\ k_3 &= (x + \frac{dx}{2}) \; (f(x) + k_2 \frac{dx}{2}) \\ k_4 &= (x + dx) \; (f(x) + k_3 dx) \\ f(x + dx) &= f(x) + \frac{1}{6} \left ( k_1 + 2k_2 + 2k_3 + k_4 \right ) \; dx \\ x &= x + dx \end{align*}\right . \] lembrando de fazer $f(x) = C$, $x = x_0$ e $dx = 1.0 \cdot 10^{-3}$ antes do loop.
----------------------------------------
3.0 - Exemplo (2ª ordem):
----------------------------------------
Para ficar mais claro, vamos voltar ao problema do oscilador harmônico simples. A EDO era da forma: \[ \ddot x(t) = - \omega_0^2 \, x(t) \] Além disso, fazendo $q_1(t) = x(t)$ e $q_2(t) = \dot x(t)$ o sistema de equações de 1ª ordem era: \[ \left \{ \begin{align*} \dot q_1(t) &= q_2(t) & (3)\\ \dot q_2(t) &= -\omega_0^2 \;\; q_1(t) & (4) \end{align*}\right . \] Agora basta aplicar o RK4 para as duas equações: \[ \left \{ \begin{align*} q_1(t + dt) = q_1(t) + \frac{1}{6} \left ( kx_1 + 2kx_2 + 2kx_3 + kx_4 \right ) dt \\ q_2(t + dt) = q_2(t) + \frac{1}{6} \left ( kv_1 + 2kv_2 + 2kv_3 + kv_4 \right ) dt \end{align*}\right . \] onde:
Note que, para o PVI são conhecidos os valores de $q_2(t)$ , $q_1(t)$, mas o valor de $q_2(t + \frac{dt}{2})$, por exemplo, não é conhecido. Contudo, lembre-se da definição, $kx_2$ deve ser a velocidade na metade do intervalo, estimada usando algum $k_1$ pelo método de Euler, note que, pela eq. (4), $kv_1$ pode ser utilizado como a inclinação de Euler.
De forma similar, pela eq. (3) $kx_1$ pode ser usado para estimar $q_1(t + \frac{dt}{2})$ em $kv_2$ por Euler.
Dessa forma, devemos aplicar o método de Euler em $kx_2$, $kx_3$, $kx_4$, $kv_2$, $kv_3$ e $kv_4$, daí:
Finalmente, retornando à nossa notação em termos de posição e velocidade, as equações de RK4 para o problema, que seriam colocadas dentro de um loop em um programa, serão: \[ \left \{ \begin{align*} kx_1 &= \dot x(t) \\ kv_1 &= -\omega_0^2 \;\; x(t) \\ kx_2 &= \dot x(t) + kv_1 \; \; \frac{dt}{2} \\ kv_2 &= -\omega_0^2 \;\; ( x(t) + kx_1 \;\; \frac{dt}{2} ) \\ kx_3 &= \dot x(t) + kv_2 \; \; \frac{dt}{2} \\ kv_3 &= -\omega_0^2 \;\; ( x(t) + kx_2 \;\; \frac{dt}{2} ) \\ kx_4 &= \dot x(t) + kv_3 \; \; dt \\ kv_4 &= -\omega_0^2 \;\; ( x(t) + kx_3 \;\; dt ) \\ x(t + dt) &= x(t) + \frac{1}{6} \left ( kx_1 + 2kx_2 + 2kx_3 + kx_4 \right ) dt \\ \dot x(t + dt) &= \dot x(t) + \frac{1}{6} \left ( kv_1 + 2kv_2 + 2kv_3 + kv_4 \right ) dt \\ t &= t + dt \end{align*}\right . \] lembrando de fazer $x(t) = x_0$, $\dot x(t) = v_0$, $t = t_0$ e $dt = 1.0 \cdot 10^{-3}$ antes do loop.
Veja um exemplo de programa em c no GitHub.
Veja um exemplo de programa em fontran no GitHub.
----------------------------------------
1.0 - O método de Runge-Kutta de 4ª ordem:
----------------------------------------
O Método de Runge-Kutta de 4ª ordem, ou simplesmente RK4, é um método de cálculo numérico que permite colher resultados com maior precisão do que se usássemos o método de Euler ou, se quisermos a mesma precisão, com menos tempo do que usando Euler.
Para o PVI: \[ \left \{ \begin{align*} \frac{d}{dx}f(x) &= F(x, \; f(x) ) \\ f(x_0) &= C \end{align*}\right . \; \; \; \; (1) \] o método de Euler de primeira ordem pode ser expresso por: \[ f(x + dx) = f(x) + I dx \] Onde $I$ é uma inclinação. No método de Euler essa inclinação é a derivada $\frac{d}{dx}f(x)$ no ponto de abscissa $x$, no RK4 essa inclinação é mais "detalhada", expressa por: \[ I=\frac{1}{6} \left ( k_1 + 2k_2 + 2k_3 + k_4 \right ) \; \; \; \; (2) \] onde: \[ \left \{ \begin{align*} k_1 &= F\left ( x, \; \; f(x) \right )\\ k_2 &= F\left ( x + \frac{dx}{2}, \; \; f(x) + k_1 \, \frac{dx}{2} \right ) \\ k_3 &= F\left ( x + \frac{dx}{2}, \; \; f(x) + k_2 \, \frac{dx}{2} \right ) \\ k_4 &= F\left ( x + dx, \; \; f(x) + k_3 \, dx \right ) \end{align*} \right . \] ou seja:
· $k_1$ é a inclinação de $f$ no início do intervalo $[x, x + dx]$, exatamente a mesma de Euler;
· $k_2$ é a inclinação de $f$ na metade do intervalo $[x, x + dx]$ estimada através do método de Euler usando $k_1$;
· $k_3$ é a inclinação de $f$ na metade do intervalo $[x, x + dx]$ estimada através do método de Euler usando $k_2$;
· $k_4$ é a inclinação de $f$ no fim do intervalo $[x, x + dx]$ estimada através do método de Euler usando $k_3$.
----------------------------------------
2.0 - Exemplo (1ª ordem):
----------------------------------------
Seja o PVI: \[ \left \{ \begin{align*} \frac{d}{dx}f(x) &= x \; f(x) \\ f(x_0) &= C \end{align*}\right . \] Para estimar o valor de $f(x_0 + dx)$ basta aplicar RK4: \[ f(x_0 + dx) = f(x_0) + \frac{1}{6} \left ( k_1 + 2k_2 + 2k_3 + k_4 \right ) dx \] com: \[ \left \{ \begin{align*} k_1 &= F(x_0, f(x_0)) = x_0 \; C \\ k_2 &= (x_0 + \frac{dx}{2}) \; (C + k_1 \frac{dx}{2}) \\ k_3 &= (x_0 + \frac{dx}{2}) \; (C + k_2 \frac{dx}{2}) \\ k_4 &= (x_0 + dx) \; (C + k_3 dx) \end{align*}\right . \] Finalmente, as equações do método RK4 para o problema, que seriam colocadas dentro de um loop em um programa, serão: \[ \left \{ \begin{align*} k_1 &= x \; f(x) \\ k_2 &= (x + \frac{dx}{2}) \; (f(x) + k_1 \frac{dx}{2}) \\ k_3 &= (x + \frac{dx}{2}) \; (f(x) + k_2 \frac{dx}{2}) \\ k_4 &= (x + dx) \; (f(x) + k_3 dx) \\ f(x + dx) &= f(x) + \frac{1}{6} \left ( k_1 + 2k_2 + 2k_3 + k_4 \right ) \; dx \\ x &= x + dx \end{align*}\right . \] lembrando de fazer $f(x) = C$, $x = x_0$ e $dx = 1.0 \cdot 10^{-3}$ antes do loop.
----------------------------------------
3.0 - Exemplo (2ª ordem):
----------------------------------------
Para ficar mais claro, vamos voltar ao problema do oscilador harmônico simples. A EDO era da forma: \[ \ddot x(t) = - \omega_0^2 \, x(t) \] Além disso, fazendo $q_1(t) = x(t)$ e $q_2(t) = \dot x(t)$ o sistema de equações de 1ª ordem era: \[ \left \{ \begin{align*} \dot q_1(t) &= q_2(t) & (3)\\ \dot q_2(t) &= -\omega_0^2 \;\; q_1(t) & (4) \end{align*}\right . \] Agora basta aplicar o RK4 para as duas equações: \[ \left \{ \begin{align*} q_1(t + dt) = q_1(t) + \frac{1}{6} \left ( kx_1 + 2kx_2 + 2kx_3 + kx_4 \right ) dt \\ q_2(t + dt) = q_2(t) + \frac{1}{6} \left ( kv_1 + 2kv_2 + 2kv_3 + kv_4 \right ) dt \end{align*}\right . \] onde:
$
\left \{ \begin{align*}
kx_1 &= q_2(t) \\
kx_2 &= q_2(t + \frac{dt}{2}) \\
kx_3 &= q_2(t + \frac{dt}{2}) \\
kx_4 &= q_2(t + dt)
\end{align*}\right .
$$\; \; \; \; \; \; \; \;$$
\left \{ \begin{align*}
kv_1 &= -\omega_0^2 \;\; q_1(t) \\
kv_2 &= -\omega_0^2 \;\; q_1(t + \frac{dt}{2}) \\
kv_3 &= -\omega_0^2 \;\; q_1(t + \frac{dt}{2}) \\
kv_4 &= -\omega_0^2 \;\; q_1(t + dt)
\end{align*}\right .
$
Note que, para o PVI são conhecidos os valores de $q_2(t)$ , $q_1(t)$, mas o valor de $q_2(t + \frac{dt}{2})$, por exemplo, não é conhecido. Contudo, lembre-se da definição, $kx_2$ deve ser a velocidade na metade do intervalo, estimada usando algum $k_1$ pelo método de Euler, note que, pela eq. (4), $kv_1$ pode ser utilizado como a inclinação de Euler.
De forma similar, pela eq. (3) $kx_1$ pode ser usado para estimar $q_1(t + \frac{dt}{2})$ em $kv_2$ por Euler.
Dessa forma, devemos aplicar o método de Euler em $kx_2$, $kx_3$, $kx_4$, $kv_2$, $kv_3$ e $kv_4$, daí:
$
\left \{ \begin{align*}
kx_1 &= q_2(t) \\
kx_2 &= q_2(t) + kv_1 \; \; \frac{dt}{2} \\
kx_3 &= q_2(t) + kv_2 \; \; \frac{dt}{2} \\
kx_4 &= q_2(t) + kv_3 \; \; dt \\
\end{align*}\right .
$$\; \; \; \; \; \; \; \;$$
\left \{ \begin{align*}
kv_1 &= -\omega_0^2 \;\; q_1(t) \\
kv_2 &= -\omega_0^2 \;\; ( q_1(t) + kx_1 \;\; \frac{dt}{2} ) \\
kv_3 &= -\omega_0^2 \;\; ( q_1(t) + kx_2 \;\; \frac{dt}{2} ) \\
kv_4 &= -\omega_0^2 \;\; ( q_1(t) + kx_3 \;\; dt )
\end{align*}\right .
$
Finalmente, retornando à nossa notação em termos de posição e velocidade, as equações de RK4 para o problema, que seriam colocadas dentro de um loop em um programa, serão: \[ \left \{ \begin{align*} kx_1 &= \dot x(t) \\ kv_1 &= -\omega_0^2 \;\; x(t) \\ kx_2 &= \dot x(t) + kv_1 \; \; \frac{dt}{2} \\ kv_2 &= -\omega_0^2 \;\; ( x(t) + kx_1 \;\; \frac{dt}{2} ) \\ kx_3 &= \dot x(t) + kv_2 \; \; \frac{dt}{2} \\ kv_3 &= -\omega_0^2 \;\; ( x(t) + kx_2 \;\; \frac{dt}{2} ) \\ kx_4 &= \dot x(t) + kv_3 \; \; dt \\ kv_4 &= -\omega_0^2 \;\; ( x(t) + kx_3 \;\; dt ) \\ x(t + dt) &= x(t) + \frac{1}{6} \left ( kx_1 + 2kx_2 + 2kx_3 + kx_4 \right ) dt \\ \dot x(t + dt) &= \dot x(t) + \frac{1}{6} \left ( kv_1 + 2kv_2 + 2kv_3 + kv_4 \right ) dt \\ t &= t + dt \end{align*}\right . \] lembrando de fazer $x(t) = x_0$, $\dot x(t) = v_0$, $t = t_0$ e $dt = 1.0 \cdot 10^{-3}$ antes do loop.
Veja um exemplo de programa em c no GitHub.
Veja um exemplo de programa em fontran no GitHub.