Uma
Equação Diferencial é uma equação que envolve derivadas
de uma ou mais funções. Elas servem para descrever o comportamento
de sistemas dinâmicos e possuem enorme aplicação em diversas
áreas, como em engenharia (no estudo do comportamento de um circuito
elétrico ou do movimento oscilatório de estruturas), Física
(lançamento de foquetes, oscilações forçadas, ...), Biologia
(crescimento de populações de bactérias) ou economia (aplicações
financeiras). Elas são classificadas de acordo com o seu tipo,
ordem e grau. Se uma equação diferencial envolve derivadas de uma
função de uma única variável independente, ela é dita ser
Equação Diferencial Ordinária. Caso
a equação diferencial envolva as derivadas parciais de uma função
de duas ou mais variáveis independentes, é uma Equação
Diferencial Parcial.
Uma equação diferencial ordinária (ou E.D.O.) de ordem n pode ser expressa na seguinte forma:
$$\frac{d^n\,y}{dx^n}=F \left(x, y,\frac{dy}{dx}, \frac{d^{2}y}{dx^{2}}, \cdots, \frac{d^{n-1}y}{dx^{n-1}} \right)$$
onde $x$ é a variável independente, $y$ é uma função desta variável independente e $\frac{d^k\,y}{dx^k}$, com $k=1,2,3, \cdots, n$ são as derivadas de $y$ em relação a $x$.
Há vários métodos que resolvem analiticamente uma E.D.O., entretanto nem sempre é possível obter uma solução analítica. Neste caso, os métodos numéricos são uma saída para se encontrar uma solução aproximada.
Nesta postagem, vamos 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 equação diferencial ordinária (ou E.D.O.) de ordem n pode ser expressa na seguinte forma:
$$\frac{d^n\,y}{dx^n}=F \left(x, y,\frac{dy}{dx}, \frac{d^{2}y}{dx^{2}}, \cdots, \frac{d^{n-1}y}{dx^{n-1}} \right)$$
onde $x$ é a variável independente, $y$ é uma função desta variável independente e $\frac{d^k\,y}{dx^k}$, com $k=1,2,3, \cdots, n$ são as derivadas de $y$ em relação a $x$.
Há vários métodos que resolvem analiticamente uma E.D.O., entretanto nem sempre é possível obter uma solução analítica. Neste caso, os métodos numéricos são uma saída para se encontrar uma solução aproximada.
Nesta postagem, vamos 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}
onde $y_0$ é um número dado. Os problemas de valor inicial (P.V.I) de ordem superior podem ser reduzidos a sistemas de primeira ordem através de variáveis auxiliares, o que permite a utilização do método numérico aqui apresentado.
Para resolver numericamente uma E.D.O. com P.V.I. (equação 1), supõem-se que ela satisfaz às condições de existência e unicidade. Esta solução numérica será encontrada para um conjunto finito de pontos (um intervalo fechado $[a,b]$) no eixo das abscissas.
Tomando-se $m$ subintervalos deste intervalo $[a, b]$, sendo $m \geq 1$, é possível determinar $m+1$ pontos onde as soluções numéricas devem ser calculados. Estes pontos $x_j \in [a, b]$ são igualmente espaçados entre si por um fator $h$, onde $x_j=x_0+j.h$, sendo que $x_0=a, \ \ x_m=b, \ \ h=(b-a)/m \ \ \mbox{e} \ \ j=0, 1, 2, \cdots, m$. O conjunto $\left\{x_0, x_1, x_2, \cdots, x_m \right\}$ obtido denomina-se rede ou malha de $[a, b]$.
Para facilitar a interpretação dos métodos, convenciona-se a seguinte notação:
- $y(x_j)$ é a solução exata do P.V.I., obtida analiticamente;
- $y_j$ é a solução numérica.
1.1 - MÉTODO DE EULER
O método de Euler, também conhecido como método da reta tangente, é um dos métodos mais antigos que se conhece para solução de equações diferenciais ordinárias. O método é muito atraente por sua simplicidade, mas não muito usado em problemas práticos, dentre muitos outros métodos, pois apesar de simples, para conseguir "boas" aproximações é necessário um número maior de cálculos.
1.1.1 - Derivação da Fórmula de Euler
Seja uma E.D.O. com P.V.I. dada pela equação 1. O que se deseja é encontrar as aproximações $y_1, y_2, \cdots, y_m$ para as soluções exatas $y(x_1), y(x_2), \cdots, y(x_m)$.
Sendo que o ponto inicial $(x_0, y_0)$ é fornecido pelo problema, o primeiro passo então é a busca de $y_1$. Para isto, aproximando-se a solução $y(x)$ por uma Série de Taylor no ponto $x=x_0$ e truncando no segundo termo, tem-se:
\begin{equation}
y(x)=y(x_0)+y'(x_0).(x-x_0)
\end{equation}
Para $x=x_1$ tem-se:
\begin{equation}
y(x_1)=y(x_0)+(x_1-x_0).y'(x_0)
\end{equation}
Lembrando que, como os valores exatos $y(x_n)$ são desconhecidos, são utilizados os valores aproximados $y_n$ e que $x_1-x_0=h$ e $y'(x_0)=f(x_0,y_0)$, onde $h$ é a distância entre os pontos consecutivos $x_n$ e $x_{n+1}$, então:
\begin{equation}
y_1=y_0+h.f(x_0, y_0)
\end{equation} Para encontrar $y_2$ na abscissa $x=x_2$ adota-se o mesmo procedimento (tomando-se ($x_1, y_1$) como ponto de partida). Assim, a solução aproximada é:
\begin{equation}
y_2=y_1+h.f(x_1, y_1)
\end{equation}
Se utilizarmos esse processo recursivas vezes chegaremos a uma fórmula geral, e esse procedimento é chamado método de Euler:
\begin{equation}
y_{n+1}=y_n+h.f(x_n, y_n)
\end{equation}
\begin{equation}
y_{n+1}=y_n+h.f(x_n, y_n)
\end{equation}
onde $x_n=x_{n-1}+h=x_0+n.h$, com $n=0, 1, 2, \cdots, m-1$.
Assim a solução numérica fornecida pelo método de Euler é a poligonal com vértices $(x_0, y_0), (x_1, y_1), (x_2, y_2), \cdots, (x_{n+1}, y_{n+1})$, onde os segmentos de reta no intervalo $(x_k, x_{k+1})$, $k=0, 1, 2, \cdots, n$ têm como equação $L_k(x)=y_k+f(x_k,y_k)(x-x_k)$.
Na figura 1, ilustramos a idéia do método de Euler. Veja que quanto menor o valor da diferença entre $x_{n+1}$ e $x_n$ menor o erro da estimativa para $y_{n+1}$.
1.1.2 - Exemplo usando o Método de Euler
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 (6):
\begin{equation}
x_{n+1}=x_n+h
\end{equation}
Figura 1: Ilustração do método de Euler |
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 (6):
\begin{equation}
x_{n+1}=x_n+h
\end{equation}
e
\begin{equation}
y_{n+1}=y_n+h.f(x_n, y_n)=y_n+h.(1-x_n+4.y_n)
\end{equation}
Com a condição inicial, $x_0=0$ e $y_0=1$, os próximos valores são calculados com o uso recursivo das equações (9) e (10). Para efeitos de comparação, usaremos dois valores para o passo h=0.01 (Figura 2) e h=0.001 (Figura 3), juntamente com a solução algébrica (8).
1.1.3 - Resultados
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.2 - Método de Euler 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 "metodo_euler.txt", podendo ser plotado usando o Origin, qtiplot, MatLab, SciLab, dentre outros.
/*==========================*/
\begin{equation}
y_{n+1}=y_n+h.f(x_n, y_n)=y_n+h.(1-x_n+4.y_n)
\end{equation}
Com a condição inicial, $x_0=0$ e $y_0=1$, os próximos valores são calculados com o uso recursivo das equações (9) e (10). Para efeitos de comparação, usaremos dois valores para o passo h=0.01 (Figura 2) e h=0.001 (Figura 3), juntamente com a solução algébrica (8).
1.1.3 - Resultados
Figura 2: Gráfico apresentando a solução numérica com h=0.01 (método de Euler) e análitica (eq. 8) |
Figura 3: Gráfico apresentando a solução numérica com h=0.001 (método de Euler) e análitica (eq. 8) |
1.2 - Método de Euler 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 "metodo_euler.txt", podendo ser plotado usando o Origin, qtiplot, MatLab, SciLab, dentre outros.
/*==========================*/
/* Método de Euler */
/*==========================*/
#include <stdio.h>
#include <stdlib.h>
#include <cmath>
FILE *fp;
static char nome[]="metodo_euler.txt";
// Função para plotar o gráfico via o aplicativo gnuplot (opcional no programa)
void scrip_gnuplot(char str[25])
{
fp=fopen("metodo_euler.plt", "w");
fprintf(fp, "reset\n");
fprintf(fp, "set title \"Metodo de Euler\" \n");
fprintf(fp, "set xlabel \"Xn\" \n");
fprintf(fp, "set ylabel \"Yn\" \n");
fprintf(fp, "set term jpeg \n" );
fprintf(fp, "set output \"Metodo_Euler.jpg\" \n");
fprintf(fp, "plot \'%s\' u 1:2 t \" Numerica \" w lp ls 7 lc 1 lw 1, \'%s\' u 1:3 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));
}
int main(void)
{
double x0, y0, xmax, h; // variaveis de entrada
double xn, yn, xn1, yn1;
// Dados de entrada:
x0=0.0;
xmax=2.0;
y0=1.0;
h=0.001; // passo h
// -----------------
fp=fopen(nome,"w");
xn=x0;
yn=y0;
while (xn<xmax)
{
fprintf(fp, "%f \t %f \t %f\n", xn, yn, y(xn));
yn1=yn+h*f(xn,yn);
xn1=xn+h;
xn=xn1;
yn=yn1;
}
fclose(fp);
scrip_gnuplot(nome); // chamada opcional no programa
system("gnuplot metodo_euler.plt"); // chamada opcional no programa
return 0;
}
#include <stdio.h>
#include <stdlib.h>
#include <cmath>
FILE *fp;
static char nome[]="metodo_euler.txt";
// Função para plotar o gráfico via o aplicativo gnuplot (opcional no programa)
void scrip_gnuplot(char str[25])
{
fp=fopen("metodo_euler.plt", "w");
fprintf(fp, "reset\n");
fprintf(fp, "set title \"Metodo de Euler\" \n");
fprintf(fp, "set xlabel \"Xn\" \n");
fprintf(fp, "set ylabel \"Yn\" \n");
fprintf(fp, "set term jpeg \n" );
fprintf(fp, "set output \"Metodo_Euler.jpg\" \n");
fprintf(fp, "plot \'%s\' u 1:2 t \" Numerica \" w lp ls 7 lc 1 lw 1, \'%s\' u 1:3 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));
}
int main(void)
{
double x0, y0, xmax, h; // variaveis de entrada
double xn, yn, xn1, yn1;
// Dados de entrada:
x0=0.0;
xmax=2.0;
y0=1.0;
h=0.001; // passo h
// -----------------
fp=fopen(nome,"w");
xn=x0;
yn=y0;
while (xn<xmax)
{
fprintf(fp, "%f \t %f \t %f\n", xn, yn, y(xn));
yn1=yn+h*f(xn,yn);
xn1=xn+h;
xn=xn1;
yn=yn1;
}
fclose(fp);
scrip_gnuplot(nome); // chamada opcional no programa
system("gnuplot metodo_euler.plt"); // chamada opcional no programa
return 0;
}
1.3 - Conclusão
Através do Método Numérico de Euler, obtemos uma boa aproximação para a solução $y(x)$ de uma E.D.O. de primeira ordem. No primeiro caso, obtivemos uma aproximação considerada "boa" em boa parte do intervalo considerado, $[0, 2]$. Já no segundo caso, obtivemos uma aproximação coincidente com a solução exata conhecida até onde se observou com passo de integração $h=0.001$, 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 Numérico de Euler é eficaz para resolução de EDOs de primeira ordem, apresentando em alguns casos erros pequenos que podem ser melhores ajustados através da escolha do passo de integração.
1.4 - 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", Depto de Eletrônica e Computação - Estudo de Casos em Engenharia Elétrica, 2007.
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.
Nenhum comentário:
Postar um comentário