sábado, 18 de abril de 2015

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

1.0 INTRODUÇÃO


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}

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.
Observação: Na solução numérica não se determina a expressão literal da função $y(x)$, mas sim uma solução aproximada do P.V.I. num conjunto discreto de pontos.


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} 
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}$. 
Figura 1: Ilustração do método de Euler
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} 
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
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)
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.


/*==========================*/
/*  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;
}

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