Simulação Numérica: Runge-Kutta 4ª ordem

Os métodos de Runge-Kutta são um dos métodos de integração mais populares por aliar simplicidade e precisão. Esses métodos são usualmente classificados segundo a ordem de precisão associada a eles. Consistem em métodos para integração de equações diferenciais ordinárias (EDO) para solução de problemas de valor inicial (PVI) do tipo:

\begin{equation} \begin{aligned} \dot{\textbf{x}} =& \textbf{f}(\textbf{x},t)\\ \textbf{x}(t_{0}&) = \textbf{x}_{0} \end{aligned} \end{equation}

onde $\textbf{x}$ é uma função desconhecida do tempo $t$ cuja solução aproximada é desejada e $\textbf{x}_{0}$ são as condições iniciais em $t_{0}$. O método de 4ª ordem (RK4) considera erros associados à quarta potência do passo de integração. Desta forma, a solução $\textbf{x}_{n+1}$ é aproximada conforme segue:

\begin{equation} \textbf{x}_{n+1} = \textbf{x}_{n} + \frac{dt}{6} (\textbf{k}_{1} + 2 \textbf{k}_{2} + 2 \textbf{k}_{3} + \textbf{k}_{4}) \end{equation}

onde $dt$ é o passo de integração e os coeficiente $\textbf{k}_{i}$ são:

\begin{equation} \begin{aligned} \textbf{k}_{1} &= \textbf{f}(\textbf{x}_{n}, t_{n}) \\ \textbf{k}_{2} &= \textbf{f}(\textbf{x}_{n} + dt \, \textbf{k}_{1}/2, t_{n} + dt/2) \\ \textbf{k}_{3} &= \textbf{f}(\textbf{x}_{n} + dt \, \textbf{k}_{2}/2, t_{n} + dt/2) \\ \textbf{k}_{4} &= \textbf{f}(\textbf{x}_{n} + dt \, \textbf{k}_{3}, t_{n} + dt) \\ \end{aligned} \end{equation}

A Figura a seguir ilusta graficamente as inclinações utilizadas no método de Runge-Kutta de 4ª ordem.

def rk4(f, x, dt, t):
    x = np.asarray(x)
    k1 = f(x, t)
    k2 = f(x+k1*dt/2.0, t+dt/2.0)
    k3 = f(x+k2*dt/2.0, t+dt/2.0)
    k4 = f(x+k3*dt, t + dt)

    return x + dt*(k1 + 2.0*(k2 + k3) + k4)/6.0

Para considerar a solução numérica, o problema é discretizado, de forma que o instante de tempo $t_{n+1}$ é definido como:

\begin{equation} t_{n+1} = t_n + dt \end{equation}

A função integrate a seguir, integra a função $f$ de $t_0$ por $n$ passos de integração $dt$ a partir da condição inicial $x_0$ e retorna uma matriz contendo a solução.

def integrate(t0, dt, n, x0, f):
    data = []
    x = x0
    t = t0
    data.append([t]+list(x))
    for i in range(n):
        x = rk4(f, x, dt, t)
        t = t + dt
        data.append([t]+list(x))

    return np.array(data)