Solução Numérica de Equações Diferenciais Ordinárias usando o Método de Runge-Kutta

1.0 INTRODUÇÃO

A resolução de equações diferenciais é um problema importantíssimo porque possui aplicações a diversas áreas do conhecimento tais como Matemática Aplicada, Física, Engenharia, Biologia, Economia e Computação Gráfica.

Devido à impossibilidade de se determinar a solução exata na maioria dos casos, desenvolveram-se técnicas de determinação de solução numérica aproximadas da equação.

Aqui, iremos nos concentrar em problemas de valor inicial (P.V.I.) em sua forma mais simples que são as equações diferenciais ordinárias de primeira ordem:

\begin{equation}  
\left\{
\begin{array}{ll}
\frac{dy}{dx}=f(x,y)\\
y(x_0)=y_0
\end{array}
\right.
\end{equation} 

Uma solução exata de um P.V.I. do tipo (1) é uma função derivável cujo gráfico passa pelo ponto $(x_0, y_0)$. Uma solução aproximada é uma tabela de valores que inicia com $(x_0, y_0)$, próximos do gráfico da função que seria a solução da equação.

Nesta postagem, nosso intuito é o de apresentar as idéias básicas e as principais propriedades dos métodos de Runge-Kutta. Uma demonstração breve do método de Runge-Kutta de $2^a$ ordem é apresentada. E, em seguida, apresentamos as regras  de implementação dos métodos de Runge-Kutta de $3^a$ e $4^a$ ordem, concluindo com exemplos práticos de sua aplicação.

1.1 - OS MÉTODOS DE RUNGE-KUTTA

Os métodos de Runge-Kutta são uma família de métodos iterativos para aproximar numericamente a solução de uma equação diferencial ordinária. Os matemáticos alemães Carl David Runge (1856-1927) e Martin Wilhelm Kutta (1867-1944), construíram os métodos de tal modo que sejam equivalentes a aproximar a solução exata de uma equação diferencial pelos primeiros $n$ termos de uma expansão em série de Taylor.

As características de eficiência e simplicidade de implementação, conferem aos métodos de Runge-Kutta grande popularidade, sendo um dos mais utilizados por engenheiros e cientistas.

1.1.1 - Propriedades dos Métodos de Runge-Kutta

          A idéia básica destes métodos é aproveitar as qualidades dos métodos de série de Taylor (ordem elevada) e ao mesmo tempo eliminar sua maior dificuldade que é o cálculo de derivadas de $f(x,y)$ que torna os métodos de série de Taylor computacionalmente inaceitáveis.

Características dos métodos de Runge-Kutta de ordem $p$:

• são de passo simples (auto-iniciantes);
• não exigem o cálculo de derivadas parciais de $f(x,y)$;
• necessitam apenas do cálculo de $f(x,y)$ em determinados pontos do intervalo $[x_{i}, x_{i+1}]$ (os quais dependem da ordem dos métodos);
• expandindo-se $f(x,y)$ por Taylor em torno de $(x_i , y_i)$ e agrupando-se os termos em relação às potências de $h$, a expressão do método de Runge-Kutta coincide com a do método de Taylor de mesma ordem.

Na postagem "Solução Numérica de Equações Diferenciais Ordinárias usando o Método de Euler", vimos que o método de Euler para resolução do P.V.I, $y'=f(x,y)$, $y(x_0)=y_0$, consiste na aplicação das fórmulas:
\begin{equation}
x_{i+1}=x_i+h
\end{equation} e \begin{equation}
y_{i+1}=y_i+k_1, \  \ \mbox{para} \  \ i=0, 1, 2, \cdots,
\end{equation}
onde $k_1=hf(x_i, y_i)$ e $h$ chamado passo ou incremento de integração é próximo de $0$.

O método de Runge-Kutta é um aperfeiçoamento do método de Euler e consiste em somar ao $y_i$ não apenas um valor de $k_1$, mas uma média de vários valores de $k_1, \,k_2, \,k_3, \cdots $.

1.1.2 - Síntese do Método de Runge-Kutta de Ordem $p$

• Fórmula geral do Método de Runge-Kutta: $$y_{i+1}=y_i + h \phi(x_i, y_i; h)$$,
• $\phi(x_i, y_i; h)$ é chamada função incremento, e pode ser interpretada como a inclinação no intervalo considerado. 
• Fórmula geral da função incremento de ordem $p$: $$\phi(x, y; h)=b_1k_1+b_2k_2+ \cdots + b_pk_p , $$ 
•  $k_1=f(x, y)$
•  $k_2=f(x+c_2h, y+ha_{21}k_1)$,
•  $k_3=f(x+c_3h, y+h(a_{31}k_1+a_{32}k_2))$, 
•  $\cdots$ 
•  $k_p=f(x+c_ph, y+h(a_{p1}k_1+ \cdots + a_{p, p-1}k_{p-1}))$.  

• $b_i$, $c_i$ e $a_{ij}$: constantes obtidas igualando-se a fórmula geral de Runge-Kutta com os termos da expansão em série de Taylor de mesma ordem.
• $k_i$: relações de recorrência (cálculo computacional eficiente).

Uma demonstração completa do método de Runge-Kutta pode ser encontrada em livros como a referência bibliográfica [2]. No que segue, limitar-nos-emos a apresentar as fórmulas necessárias.

1.2 - Métodos de Runge-Kutta de $2^a$ Ordem

        O método de Runge-Kutta de $2^a$ ordem concorda com a precisão do método de série de Taylor até os termos de $2^a$ ordem em h: $$y_{i+1}=y_i+hf(x_i, y_i) + \frac{h^2}{2}\left[\frac{\partial f(x_i, y_i)}{\partial x} + f(x_i, y_i)\frac{\partial f(x_i, y_i)}{\partial y} \right]$$
Usando a notação simplificada, $f_i=f(x_i, y_i)$, temos:
\begin{equation}
y_{i+1}=y_i+hf_i + \frac{h^2}{2}\left[\frac{\partial f_i}{\partial x} + f_i \frac{\partial f_i}{\partial y} \right]
\end{equation}
Considerando, agora, a definição do método de Runge-Kutta de $2^a$ ordem:
\begin{equation}
y_{i+1}=y_i+h(b_{1}f(x_i, y_i) + b_{2}f(x_{i}+c_{2}h, y_{i}+ha_{21}k_1)
\end{equation}
Expandindo $f(x, y)$, em série de Taylor, em torno de $(x_i, y_i)$ (retendo somente os termos de derivada primeira), temos:
$$f(x_{i}+c_{2}h, y_{i}+ha_{21}k_1) \approx  f_i + c_{2}h \frac{\partial f_i}{\partial x}+a_{21}h f_i \frac{\partial f_i}{\partial y}$$
Substituindo esse resultado na equação anterior, e rearranjando, teremos:
\begin{equation}
y_{i+1}=y_i+h(b_1 + b_2)f_i + h^2 \left(b_{2}c_{2}\frac{\partial f_i}{\partial x} + b_{2}a_{21}f_i \frac{\partial f_i}{\partial y} \right)
\end{equation}
Assim, para haver concordância com o método de série de Taylor, resultado (4), é preciso que: $$b_1+b_2 = 1 \  \  \  \mbox{e} \  \  \ b_{2}c_{2}=\frac{1}{2} \ \  \  \mbox{e} \ \  \ b_{2}a_{21}=\frac{1}{2}.$$

• Resultado: Sistema não linear com 3 equações e 4 incógnitas, o que corresponde a uma variedade de métodos de segunda ordem. 

Apresentaremos aqui os dois mais conhecidos métodos de Runge-Kutta de $2^a$ ordem.

1.2.1 - Métodos de Euler Modificado ($b_2=1$, $b_1=0$, $c_2=a_{21}=1/2$)

Dado um P.V.I, $y' = f(x, y)$,  $y(x_0)=x_0$. A implementação deste método é dada pelas expressões:

• $y_{i+1} = y_i + h.k_2$
• $k_1 = f(x_i, y_i)$
• $k_2 = f(x_{i}+\frac{h}{2}, y_{i}+\frac{h}{2}k_{1})$
• com $x_{i+1} = x_i + h$,  para $h>0$ próximo de $0$,  $i=0, 1, 2, \cdots $ .

Para cada valor inteiro de $i$, a partir de $i=0$, calculam-se: $$x_{i+1}  \ \ \rightarrow \ \  k_1 \  \ \rightarrow \  \  k_2  \  \  \rightarrow \  \ y_{i+1}$$
Repete-se essa sequência de cálculos recursivamente, até chegar no valor de $y_i$ desejado.

1.2.2 - Métodos de Euler Melhorado ou de Heun ($b_2=1/2$, $b_1=1/2$, $c_2=a_{21}=1$)

Dado um P.V.I, $y' = f(x, y)$,  $y(x_0)=x_0$. A implementação deste método é dada pelas expressões:

• $y_{i+1} = y_i + \frac{h}{2}\left(k_1 + k_2 \right)$
• $k_1 = f(x_i, y_i)$
• $k_2 = f(x_{i}+h, y_{i}+hk_{1})$
• com $x_{i+1} = x_i + h$,  para $h>0$ próximo de $0$,  $i=0, 1, 2, \cdots $ .

Para cada valor inteiro de $i$, a partir de $i=0$, calculam-se: $$x_{i+1}  \ \ \rightarrow \ \  k_1 \  \ \rightarrow \  \  k_2  \  \  \rightarrow \  \ y_{i+1}$$
Repete-se essa sequência de cálculos recursivamente, até chegar no valor de $y_i$ desejado.


1.3 - Métodos de Runge-Kutta de Ordem Superiores

      De modo semelhante, podem ser deduzidas as fórmulas de Runge-Kutta de ordens superiores.
      Em cada ordem, também haverá uma infinidade de versões. A seguir serão fornecidas apenas as fórmulas mais conhecidas para os métodos de Runge-Kutta de $3^a$ e $4^a$ ordens:

1.3.1 - Métodos Runge-Kutta de $3^a$ Ordem

Dado um P.V.I, $y' = f(x, y)$,  $y(x_0)=x_0$. A implementação deste método é dada pelas expressões:

• $y_{i+1} = y_i + \frac{h}{6}\left(k_1 + 4k_2 + k_3 \right)$
• $k_1 = f(x_i, y_i)$
• $k_2 = f(x_{i}+\frac{h}{2}, y_{i}+\frac{h}{2}k_{1})$ 
• $k_3 = f(x_{i}+h, y_{i} - hk_{1} + 2hk_{2})$ 
• com $x_{i+1} = x_i + h$,  para $h>0$ próximo de $0$,  $i=0, 1, 2, \cdots $.

Para cada valor inteiro de $i$, a partir de $i=0$, calculam-se: $$x_{i+1}  \ \ \rightarrow \ \  k_1 \  \  \rightarrow \  \  k_2  \  \  \rightarrow  \  \  k_3  \  \  \rightarrow \  \ y_{i+1}$$
Repete-se essa sequência de cálculos recursivamente, até chegar no valor de $y_i$ desejado.


1.3.2 - Métodos Runge-Kutta de $4^a$ Ordem

Dado um P.V.I, $y' = f(x, y)$,  $y(x_0)=x_0$. A implementação deste método é dada pelas expressões:

• $y_{i+1} = y_i + \frac{h}{6}\left(k_1 + 2k_2 + 2k_3 + k_4 \right)$
• $k_1 = f(x_i, y_i)$
• $k_2 = f(x_{i}+\frac{h}{2}, y_{i}+\frac{h}{2}k_{1})$ 
• $k_3 = f(x_{i}+\frac{h}{2}, y_{i}+\frac{h}{2}k_{2})$
• $k_4 = f(x_{i}+h, y_{i} + hk_{3})$ 
• com $x_{i+1} = x_i + h$,  para $h>0$ próximo de $0$,  $i=0, 1, 2, \cdots $.

Para cada valor inteiro de $i$, a partir de $i=0$, calculam-se: $$x_{i+1}  \ \ \rightarrow \ \  k_1 \  \  \rightarrow \  \  k_2  \  \  \rightarrow  \  \  k_3  \  \  \rightarrow   \  \  k_4  \  \  \rightarrow \  \   y_{i+1}$$
Repete-se essa sequência de cálculos recursivamente, até chegar no valor de $y_i$ desejado.

 Nota:  O método de Runge-Kutta de $4^a$ ordem é o mais usado na solução numérica de problemas com equações diferenciais ordinárias.


1.4 - Exemplo usando o Método de Runge-Kutta de $2^a$ e $4^a$ Ordem

        A equação diferencial ordinária de $1^{a}$ ordem escolhida será (Boyce, W. E., Diprima, R. C. Equações Diferenciais Elementares e Problemas de Valor de Contorno. Ed. 7, pp. 420): 
\begin{equation}  
\left\{
\begin{array}{ll}
\frac{dy}{dx}=1-x+4.y\\
y(0)=1
\end{array}
\right.
\end{equation}
Esta equação será resolvida de $x=0$ a $x=2$. A solução da equação diferencial (7) com a condição inicial (P.V.I.) é conhecida:
\begin{equation}
y(x)=\frac{1}{4}.x-\frac{3}{16}+\frac{19}{16}e^{4.x}
\end{equation}
A solução numérica é encontrada com a avaliação das equações do Método de Euler (Runge-Kutta de $1^a$ ordem), Euler modificado (Runge-Kutta de $2^a$ ordem) e do Método de Runge-Kutta de $4^a$ ordem, apresentadas acima, com a condição inicial, $x_0=0$ e $y_0=1$.

Para efeitos de comparação, usaremos dois valores para o passo h=0.1 (Figura 2) e h=0.01 (Figura 3), juntamente com a solução algébrica (8).


1.4.1 - Resultados obtidos pelos Métodos de Euler, Euler modificado e de Runge-Kutta de $4^a$ ordem

Figura 2: Métodos de Runge-Kutta com passo de integração $h=0.1$.

Figura 3: Métodos de Runge-Kutta com passo de integração $h=0.01$.

Por torna-se um processo demasiadamente repetitivo pelo número de iterações, é comum utilizar o método com o auxílio de uma linguagem de programação.

 
1.5 - Métodos de Runge-Kutta implementado em linguagem de programação. 

Um programa escrito em C/C++ é apresentado a seguir. Uma parte do código-fonte pode ser retirada sem problemas, pois trata-se de rotinas voltadas para mostrar a interatividade entre o aplicativo gnuplot (para plotar gráficos) e a linguagem programação C/C++. Os dados obtidos são impressos em um arquivo de nome "runge_kutta.txt", podendo ser plotado usando o Origin, qtiplot, MatLab, SciLab, dentre outros.


/*=================== */
/*  Métodos de Runge-Kutta          */
/*=================== */

#include <stdio.h>
#include <stdlib.h>
#include <cmath>

FILE *fp;

static char nome[]="runge_kutta.txt";

void scrip_gnuplot(char str[25])
{
    fp=fopen("rk_2_4.plt", "w");
   
    fprintf(fp, "reset\n");
    fprintf(fp, "set title \"Metodos de Runge Kutta\" \n");
    fprintf(fp, "set xlabel \"Xn\" \n");
    fprintf(fp, "set ylabel \"Yn\" \n");
    fprintf(fp, "set key left \n");
    fprintf(fp, "set term jpeg \n" );
    fprintf(fp, "set output \"RK_2_4.jpeg\" \n");
    fprintf(fp, "plot \'%s\' u 1:2 t \" Euler \" w lp ls 7 lc 1 lw 1, \'%s\' u 1:3 t \" Euler Modificado \" w lp ls 7 lc 2 lw 1\n", str, str);
    fprintf(fp, "replot \'%s\' u 1:4 t \" Runge-Kutta 4a ordem \" w lp ls 7 lc 3 lw 1, \'%s\' u 1:5 t \" Analitica \" w l lc -1\n", str, str);
    fprintf(fp, "set output; set term wxt \n");
    fprintf(fp, "replot \n");
    fprintf(fp, "pause -1 \"Continuar?\" ");
    fclose(fp);
}
// Definição da equação diferencial: dy/dx=f(x,y)
double f(double x, double y)
{
    return (1-x+4*y);
}
// Definiçao da solução analitica da equação diferencial
double y(double x)
{
    return ((1.0/4.0)*x - (3.0/16.0) + (19.0/16.0)*exp(4.0*x));
}
// Definiçao do método de euler
double euler(double x, double y, double h)
{
double k1;

        k1 = f(x,y);
        y = y + h*k1;
       
return y;       
}
// Definiçao do método de euler modificado
double euler_modificado(double x, double y, double h)
{
double k1, k2;

        k1 = f(x,y);
        k2 = f(x+0.5*h, y+0.5*h*k1);
        y = y + h*k2;
       
return y;       
}
// Definiçao do método de runge-kutta de 4a ordem
double runge_kutta_4(double x, double y, double h)
{
double k1, k2, k3, k4;

        k1 = f(x,y);
        k2 = f(x+0.5*h, y+0.5*h*k1);
        k3 = f(x+0.5*h, y+0.5*h*k2);
        k4 = f(x+h, y+h*k3);
        y = y + h*(k1+2*k2+2*k3+k4)/6.0;

return y;       
}

int main(void)
{
double x0, y0, xmax, h;           // variaveis de entrada
double xn, yn_e, yn_em, yn_rk4;
double y_e, y_em, y_rk4;

// Dados de entrada:
x0=0.0;
xmax=2.0;
y0=1.0;
h=0.01;
// -----------------
    fp=fopen(nome,"w");
   
    xn=x0;
    yn_e=y0;
    yn_em=y0;
    yn_rk4=y0;
    while (xn<xmax)
    {
        fprintf(fp, "%f \t %f \t %f \t %f \t %f\n", xn, yn_e, yn_em, yn_rk4, y(xn));
       
        y_e = euler(xn, yn_e, h);
        y_em =euler_modificado(xn, yn_em, h);
        y_rk4=runge_kutta_4(xn, yn_rk4, h);
        xn=xn+h;
        yn_e = y_e;
        yn_em =y_em;
        yn_rk4=y_rk4;
    }
    fclose(fp);

    scrip_gnuplot(nome);
    system("gnuplot rk_2_4.plt");
    return 0;
}

1.6 - Conclusão

        A resolução numérica pode ter grau crescente de aproximação, de acordo com o método numérico escolhido e da escolha do valor do passo de integração. Para os valores do passo de integração analisados, observamos que o método de Euler não se mostra uma boa escolha, pois seus resultados não são nada precisos, quando comparado aos dois outros métodos. Já no primeiro caso, obtivemos uma boa aproximação por parte do método de Runge-Kutta de $4^a$ ordem. E, no segundo caso, já se verifica uma coincidência entre os métodos de Euler modificado e de Runge-Kutta de $4^a$ ordem com solução exata conhecida até onde se observou com passo de integração $h=0.01$, ou seja, o "erro" ou diferença entre a aproximação e a solução exata é praticamente nula para todos os pontos.

Assim, diante dos exemplos estudados, concluímos que o Método de Runge-Kutta de $4^a$ ordem é bastante eficaz para resolução de EDOs de primeira ordem, o que explica o porquê de ser o método de passo simples mais utilizado por engenheiros e cientistas na resolução numérica de equações diferenciais ordinárias.

1.7 - Referências

[1]  W. E. Boyce, C. R. Diprima, "Equações Diferenciais Elementares e Problemas de Valor de Contorno", LTC, Rio de Janeiro, 2011.

[2]  M. A. G. Ruggiero, V. L. R. Lopes, "Cálculo Numérico: Aspectos Teóricos Computacionais", Makron Books, São Paulo, 1996.

[3]  G. Baratto, "Solução de Equações Diferenciais Ordinárias Usando Métodos Numéricos", Departamento de Eletrônica e Computação - Estudo de Casos em Engenharia Elétrica - UFSM, 2007.
[4] L. N. de Andrade, "Cálculo Numérico: Introdução à Matemática Computacional",  Departamento de Matemática - UFPB, 2013.

Nesta postagem, o nosso objetivo, além de apresentar este importante tema da Matemática, é o de fornecer exemplos da interatividade entre as linguagem de programação C/C++ e o aplicativo gnuplot, que juntos minimizam as tarefas de uma pesquisa.