Por volta do ano 1900, dois matemáticos alemães, Carl David Tolmé Runge e Martin Wilhelm Kutta, aprimoraram os métodos iterativos para resoluções de Equações Diferenciais Ordinárias de uma forma extremamente precisa.

Como Heun, eles perceberam que quanto mais inclinações (derivadas) fossem calculadas, melhor seria a precisão, sendo cada inclinação associada a um determinado coeficiente de compensação.

Assim, analisando o trabalho de Heun, foi possível dar continuidade àquele raciocínio já desenvolvido pelo Alemão.

A versão de segunda ordem para o método de RK segue a seguinte equação:
$$ {y_{i+1} = y_i + (a_1k_1 + a_2k_2)h}, $$

em que
$$k_1 = f(x_i,y_i),$$ $$k_2 = f(x_i + p_1h, y_i + q_{11}k_1h).$$
Os valores $a_1, a_2, p_1$ e $q_{11}$ são calculados igualando-se a primeira equação à expansão em série de Taylor até os termos de 2º grau. Fazendo isso, deduzimos três equações para calcular as quatro constantes desconhecidas. As três equações são:

$${a_1 + a_2 = 1,}$$ $${a_2p_1 = \frac{1}{2}},$$ $${a_2q_{11} = \frac{1}{2}}.$$
Nós podemos adotar quaisquer valores para as constantes $a_2$ ou $a_1$, posteriormente fazendo o valor adotado obedecer o conjunto de equações, ou seja, nós teremos infinitas maneiras de escrever o método de segunda ordem de RK.

Se fizermos, por exemplo, $a_1 = 0$, $a_2$ será igual a $1. Então teremos

$$y_{i+1} = y_i + (a_1k_1 + a_2k_2)h.$$ Substituindo $a_1$ e $a_2$,
$$y_{i+1} = y_i + k_2h,$$ e nossas inclinações serão:

$$k_1 = f(x_i,y_i).$$
Para encontrarmos $k_2$ teremos que resolver as equações que restaram. Assim sendo, temos que

$$a_2p_1 = \frac{1}{2}.$$ Subtituindo o valor de $a_2 (a_2 = 1)$,
$$p_1 = \frac{1}{2}$$ e $$q_{11} = \frac{1}{2}.$$ Dessa forma, a segunda inclinação será: $$k_2 = f(x_i + \frac{1}{2}h, y_i + \frac{1}{2}k_1h).$$ Perceba que essa maneira que escrevemos o método de segunda ordem, na verdade, corresponde ao método do Ponto Médio, que já conhecemos.

Se fizermos agora $a_1 = \frac{1}{2}$, a partir disso (usando o conjunto de equações) podemos calcular que $a_2 = \frac{1}{2}$. Isso nos dará a seguinte equação para o calculo do próximo valor de $y$: $$ {y_{i+1} = y_i + (\frac{1}{2}k_1 + \frac{1}{2}k_2)h}, $$ $$k_1 = f(x_i, y_i),$$ e $a_1p_1 = \frac{1}{2}$. Então, $p_1 = 1$ e $a_2q_{11} = \frac{1}{2}$ implica que $q_{11} = 1$. Sendo assim, $$k_2 = f(x_i + h, y_i + q_{11}k_1h).$$ Essa expressão para o método de RK de segunda ordem também já vimos, ela é o método desenvolvido por $Heun$.

Dessa forma, podemos ver que realmente há infinitas formas de escrevermos a expressão para calcular numéricamente uma EDO usando o método de RK2, sendo Ponto Médio e Heun casos particulares dos métodos de Runge-Kutta.

A dedução do Método de Runge-Kutta de Terceira Ordem é bem parecida com a dos métodos de segunda ordem. O resultado é um conjunto de seis equações com oito incógnitas. Portanto, devem ser especificados valores, a priori, para as duas incógnitas (como fizemos em RK2) para determinar os parâmetros restantes. Uma, dentre as vesões mais comuns para RK3 é:

$$ {y_{i+1} = y_i + \frac{1}{6}(k_1 + 4k_2 + k_3)h}, $$ em que

Nesta versão, $k_1$ e $k3$ tem o peso $1$ e $k_2$ tem peso $4$. Devemos lembrar que esse método se trata de um caso particular dos métodos de RK de 3ª ordem.

Os métodos de RK4 são os de maior precisão que usaremos aqui no curso. Esses métodos já nos dão uma boa aproximação comparada à resolução analítica da EDO. A seguir, veremos a forma mais comumente usada, que ficou conhecida como método RK de quarta ordem clássico:

Podemos observar que a ordem do método, na verdade, corresponde ao número de inclinações que teremos de calcular, ou seja, fica claro que para o método de RK4, nós teremos quatro inclinações com novos pesos:

$$ {k_1 = f(x_i,y_i)}, $$ $$ \begin{aligned} k_2 &= f(x_{i+1}, \tilde{y}_{i+1}) \\ &= f(x_i + \frac{1}{2}h, y_i + \frac{1}{2}k_1 h), \end{aligned} $$ $$ \begin{aligned} k_3 &= f(x_{i+1}, \tilde{y}_{i+1}) \\ &= f(x_i + \frac{1}{2}h, y_i + \frac{1}{2}k_2 h), \end{aligned} $$ $$ \begin{aligned} k_4 &= f(x_{i+1}, \tilde{y}_{i+1}) \\ &= f(x_i + h, y_i + k_3h). \end{aligned} $$ A inclinação resultante será: $$ {k = \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4)}. $$ Sendo assim, já podemos calcular o próximo valor de $y$:
$$ {y_{i+1} = y_i + kh}. $$ Ilustrando o método graficamente:

Usando o método de Runge Kutta de quarta ordem, resolva a seguinte equação diferencial
de $a = 0$ até $b = 2.5$ com um passo $h = 0.5$ e condição inicial $y = 3$ em $x = 0$. Compare o resultado com a solução analítica: Solução:

O primeiro ponto da solução é $(0;3)$, que corresponde a condição inicial. Ou seja, no ponto ${i = 1}$ temos ${x_1 = 0}$ e ${y_1 = 3}$. Para o segundo ponto, calcularemos primeiro as inclinações: Logo, Comparando com o resultado da solução analítica:

O restante da solução podemos ver na tabela a seguir: