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)