Sistemas de Equações Diferenciais Ordinárias: Métodos de Runge-Kutta de Ordem p

1.0 Introdução

Em duas postagens anteriores, apresentamos os métodos numéricos de Runge-Kutta com o intuito de resolver equações diferenciais de primeira ordem, mais especificamente em problemas de valor inicial (P.V.I) de primeira ordem. Entretanto, a maioria das equações diferenciais com importância prática, são de ordem maior que $1$ ou então são sistemas de equações diferenciais. Nesta postagem, nosso intuito é o de apresentar o procedimento de como resolver um sistema de equações diferenciais de primeira ordem e, também, o de resolver numericamente uma equação diferencial de ordem elevada, usando os métodos de Runge-Kutta.

1.1 - Sistemas de Equações Diferenciais Ordinárias

Consideremos um sistema de $n$ equações diferenciais de primeira ordem:

\begin{equation}
\begin{array}{llll}
y'_1=f_{1}(x, y_1, y_2, \cdots, y_n)\\
y'_2=f_{2}(x, y_1, y_2, \cdots, y_n)\\
\vdots\\
y'_n=f_{n}(x, y_1, y_2, \cdots, y_n)
\end{array}
\end{equation}

o qual pode ser escrito, como: $$\mathbf{y' =  f(x, y)} , $$ onde $\mathbf{y}$, $\mathbf{y'}$ e $\mathbf{f}$ são vetores com componentes $\mathbf{y_i}$, $\mathbf{y'_i}$ e $\mathbf{f_{i}(i=1, 2, \cdots, n)}$, respectivamente. Para que esse sistema possua uma única solução é necessário impormos uma condição adicional sobre $\mathbf{y}$. Esta condição é usualmente da forma: $$\mathbf{y(x_{0}) = y_{0}} , $$ para um dado número $\mathbf{x_0}$ e um vetor $\mathbf{y_0}$

Agora descreveremos como os métodos apresentados nas postagens anteriores para a solução de equações diferenciais de primeira ordem podem ser aplicados para resolver sistemas de equações diferenciais de $1^a$ ordem. Para efeito de simplicidade, e sem perda de generalidade, consideramos apenas o caso em que $n=2$, isto é, o sistema possui apenas duas equações, e para maior clareza usaremos a notação:

\begin{equation}
\left\{
\begin{array}{llll}
y' =  f(x, y, z)\\
z' = g(x, y, z)\\
y(x_{0})=y_{0}\\
z(x_{0})=z_{0} \  \  \  x \in \left[x_{0}, b \right]
\end{array}
\right.
\end{equation}

Assim, se desejamos resolver o sistema $(2)$ pelo método de Euler (Runge-Kutta de ordem 1), teremos:

\begin{equation}
\begin{array}{ll}
y_{i+1}=y_i + hf(x_i, y_i, z_i)\\
z_{i+1}=z_i + hg(x_i, y_i, z_i)\\
\end{array}
\end{equation}
que será aplicado passo a passo.

No caso de escolhermos o método de Runge-Kutta de $4^a$ ordem para resolvermos o sistema $(2)$, teremos as seguintes expressões:

\begin{equation}
\begin{array}{ll}
y_{i+1}=y_i + \frac{h}{6}(k_1 + 2k_2 + 2k_3 +k_4)\\
z_{i+1}=z_i + \frac{h}{6}(l_1 + 2l_2 + 2l_3 + l_4)\\
\end{array}
\end{equation}
onde

• $k_1 = f(x_i, y_i, z_i)$
• $l_1 = g(x_i, y_i, z_i)$
• $k_2 = f(x_{i}+\frac{h}{2}, y_{i}+\frac{h}{2}k_{1}, z_i + \frac{h}{2}l_1)$ 
• $l_2 = g(x_{i}+\frac{h}{2}, y_{i}+\frac{h}{2}k_{1}, z_i + \frac{h}{2}l_1)$ 
• $k_3 = f(x_{i}+\frac{h}{2}, y_{i}+\frac{h}{2}k_{2}, z_i + \frac{h}{2}l_2)$
• $l_3 = g(x_{i}+\frac{h}{2}, y_{i}+\frac{h}{2}k_{2}, z_i + \frac{h}{2}l_2)$
• $k_4 = f(x_{i}+h, y_{i} + hk_{3}, z_i + l_3)$  
• $l_4 = g(x_{i}+h, y_{i} + hk_{3}, z_i + l_3)$  

Nota: A aplicação de um método numérico a um sistema de equações ordinárias de primeira ordem se processa como no caso de uma única equação, só que aqui devemos aplicar o método numérico a cada uma das componentes do vetor.

1.2 - Equações Diferenciais de Ordem Elevada

Finalmente, mostraremos como equações de ordem superiores a 1 podem ser escritas e portanto resolvidas como um sistema de equações de primeira ordem. Consideremos a equação diferencial de ordem $n$: $$y^{(n)} = f(x, y, y', \cdots, y^{(n-1)}).$$

com as condições iniciais: $$y(x_0) = y_0, \  \  \  y'(x_0) = y'_0,  \  \  \cdots, \  \ y^{(n-1)}(x_0)=y^{(n-1)}_0.$$

Novamente, para simplicidade, mas sem perda de generalidade, consideremos a equação diferencial de segunda ordem:

\begin{equation}
\left\{
\begin{array}{lll}
\frac{d^2\,y}{dx^2}+q(x)\frac{dy}{dx} = r(x,y)\\

y(x_{0})=y_{0}\\
\frac{d\,y}{dx}|_{x_0}=y'_{0}
\end{array}
\right.
\end{equation}

Podemos resolver qualquer equação diferencial de ordem elevada reduzindo-a a um sistema de equações diferenciais de primeira ordem. Para tanto basta fazer uma mudança adequada de variável, para a equação de segunda ordem, (5). Fazendo a seguinte mudança de variável: $$z(x, y) = \frac{d\,y}{dx}, $$

a equação (5) é escrita como um conjunto de duas equações diferenciais de $1^a$ ordem, ou seja, o P.V.I, definido por (5), se reduz a:
\begin{equation}
\left\{
\begin{array}{llll}
\frac{d\,y}{dx} =  z(x, y)\\
\frac{d\,z}{dx}= r(x, y) - q(x)z(x, y)\\
y(x_{0})=y_{0}\\
z(x_{0})=z_{0}
\end{array}
\right.
\end{equation}

Observe que o conjunto de equações em (6) pode ser escrito na forma:

\begin{equation}
\left\{
\begin{array}{llll}
y' = f(x, y, z)\\
z' = g(x, y, z)\\
y(x_{0})=y_{0}\\
z(x_{0})=z_{0}
\end{array}
\right.
\end{equation}

Nota: Ao reduzir o P.V.I. definido pela equação $(5)$, a um sistema de E.D.O. de primeira ordem, o problema passa a ser análogo ao P.V.I definido pelo sistema (2).

Portanto, para a solução da equação diferencial ordinária de $2^a$ ordem $(5)$, é necessário que as equações diferenciais ordinárias de $1^a$ ordem, apresentadas em $(7)$ sejam resolvidas. A seguir, apresenta-se a solução destas equações, usando o método de Euler e o Método de Runge-Kutta de $4^a$ ordem.

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

        A equação diferencial ordinária de $2^{a}$ ordem escolhida será:

\begin{equation}
\frac{d^{2}y}{dx^{2}}-y=e^{x}
\end{equation}
com $x$ variando de $0$ a $1$, com as seguintes condições iniciais:
\begin{equation}
\left\{
\begin{array}{ll}
y(0)=1\\
\frac{dy}{dx}|_{x_0}=0
\end{array}
\right.
\end{equation}

Primeiro passo: é reduzir esta equação diferencial de $2^a$ ordem para um conjunto de $2$ equações diferenciais ordinárias de $1^a$ ordem:

\begin{equation}
\left\{
\begin{array}{ll}
\frac{dy}{dx} = z \\
\frac{dz}{dx} =y+e^{x}
\end{array}
\right.
\end{equation}

De acordo com as condições iniciais em $(9)$, tem-se $x_0 = 0$ e $y_0 = 1$. 

A solução da equação diferencial ordinária de $2^a$ ordem, $(8)$, com a condição inicial $(9)$ é:

\begin{equation}
y(x)=\frac{1}{4}[e^{x}(1+2x)+3e^{-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), e do Método de Runge-Kutta de $4^a$ ordem, apresentadas respectivamente pelas equações (3) e (4), com a condição inicial, $x_0=0$ e $y_0=0$.

Nota: vale lembrar que, a solução aproximada da equação diferencial de segunda ordem encontra-se na primeira componente do vetor, isto é, apenas nos interessa o valor de $y_i$, apesar de termos de calcular, a cada passo, todas as componentes do vetor.

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

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

Figura 1: Solução da E.D.O. de $2^a$ ordem com passo $h=0.1$

Figura 2: Solução da E.D.O. de $2^a$ ordem com passo $h=0.01$

Figura 3: Solução da E.D.O. de $2^a$ ordem com passo $h=0.001$

1.4 - 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 - E.D.O. de 2ª ordem   */
/*=============================*/
 
#include <stdio.h>
#include <stdlib.h>
#include <cmath>

FILE *fp;
static char nome[]="runge_kutta.txt";

void scrip_gnuplot(char str[25])
{
    fp=fopen("rk4_2_ordem.plt", "w");
    fprintf(fp, "reset\n");
    fprintf(fp, "set title \"Métodos de Runge Kutta\" \n");
    fprintf(fp, "set xlabel \"Xn\" \n");
    fprintf(fp, "set ylabel \"Yn\" \n");
    fprintf(fp, "set key bottom \n" );
    fprintf(fp, "set term jpeg \n" );
    fprintf(fp, "set output \"Runge_Kutta.jpg\" \n");
    fprintf(fp, "plot \'%s\' u 1:2 t \" Euler \" w lp ls 7 lc 1 lw 1, \'%s\' u 1:3 t \" Runge-Kutta 4a ordem \" w lp ls 7 lc 2 lw 1\n", str, str);
    fprintf(fp, "replot \'%s\' u 1:4 t \" Analitica \" w l lc -1\n", str);
    fprintf(fp, "set output; set term wxt \n");
    fprintf(fp, "replot \n");
    fprintf(fp, "pause -1 \"Continuar?\" ");
    fclose(fp);
}

double f(double x, double y, double z)
{
    return z;                        // equação diferencial: dy/dx=f(x,y,z)
}

double g(double x, double y, double z)
{                                       
    return (y+exp(x));        // equação diferencial: dz/dx=g(x,y,z)
}

double y(double x)           // solução analitica da equação diferencial
{
    return (0.25*(exp(x)*(1.0+2.0*x)+3.0*exp(-x)));
}

double euler(double x, double y, double *z, double h)
{
    y = y + h*f(x,y,*z);        // método de euler
    *z = *z + h*g(x,y,*z);
return y;       
}

double runge_kutta_4(double x, double y, double *z, double h)
{
double k1, k2, k3, k4, l1, l2, l3, l4;

    k1 = h*f(x, y, *z);         // método de runge-kutta de ordem 4
    l1 = h*g(x, y, *z);
    k2 = h*f(x+0.5*h, y+0.5*k1, *z+0.5*l1);
    l2 = h*g(x+0.5*h, y+0.5*k1, *z+0.5*l1);
    k3 = h*f(x+0.5*h, y+0.5*k2, *z+0.5*l2);
    l3 = h*g(x+0.5*h, y+0.5*k2, *z+0.5*l2);
    k4 = h*f(x+h, y+k3, *z+l3);
    l4 = h*g(x+h, y+k3, *z+l3);
    y = y + (k1+2*k2+2*k3+k4)/6.0;
    *z = *z + (l1+2*l2+2*l3+l4)/6.0;
return y;       
}

int main(void)
{
double x0, y0, z0, xmax, h;                          // variaveis de entrada

double xn, yn_e, yn_rk4, *zn_e, *zn_rk4;   // variaveis incrementadas
double z0_e, z0_rk4;

// Dados de entrada:
x0=0.0;
y0=1.0;
z0=0.0;
xmax=1.0;
h=0.01;
// -----------------
   
    fp=fopen(nome,"w");

    xn = x0;
    yn_e = yn_rk4 = y0;
   
    z0_e = z0_rk4 = z0;
   
    zn_e = &z0_e;
    zn_rk4 = &z0_rk4;
    while (xn<xmax)
    {
        fprintf(fp, "%f \t %f \t %f \t %f\n", xn, yn_e, yn_rk4, y(xn));
       
        yn_e = euler(xn, yn_e, zn_e, h);
        yn_rk4 = runge_kutta_4(xn, yn_rk4, zn_rk4, h);
        xn = xn+h;
    }
    fclose(fp);
    scrip_gnuplot(nome);
    system("gnuplot rk4_2_ordem.plt");
    return 0;
}

1.5 - 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. Comparando os resultados da solução numérica usando o método de Euler com os resultados da solução usando o método de Runge-Kutta de $4^a$ ordem, observa-se que neste segundo método a precisão é sempre maior, mesmo com o uso de um passo $10$ vezes maior.

Assim, diante dos exemplos estudados, concluímos que o Método de Runge-Kutta de $4^a$ ordem é bastante eficaz para resolução de sistemas 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.6 - 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.