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.
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.
[4] L. N. de Andrade, "Cálculo Numérico: Introdução à Matemática Computacional", Departamento de Matemática - UFPB, 2013.