Método de Euler

$\newcommand{\R}{\mathbb{R}}$

O método de Euler é construído para aproximar a solução do PVI \begin{equation}\label{standard} y'(t) = f(t,y), \quad t>t_0, \qquad y(t_0)=y_0, \end{equation} onde $y:\R\rightarrow \R^m,$ $f:\R\times\R^m\rightarrow \R^m$ e $y_0\in\R^m.$ Essa aproximação será obtida apenas em alguns valores discretos de $t,$ ou seja apenas em $t_0 \lt t_1 \lt \cdots \lt t_N.$ Esses pontos são ditos a malha numérica. Um aproximação para a função $y$ em $t=t_n$ será representada por $y_n,$ ou seja $$ y_n \approx y(t_n), \quad n=0,1\ldots, N. $$

O método de Euler é um método de passo simples, ou seja, um método que para determinar a aproximação em um ponto da malha utiliza apenas informações do passo anterior. Em oposição a métodos de passo simples, há os métodos de passos múltiplos, nos quais informações obtidas em dois ou mais passos anteriores são utilizadas para definir a aproximação no passo seguinte.

A relação entre os valores da função $y$ em pontos consecutivos da malha é feita através da polinômio de Taylor de primeira ordem, ou seja, \begin{align} y(t_{n+1}) &= y(t_n) + y'(t_n)(t_{n+1}-t_n) + y''(\tau){(t_{n+1}-t_n)^2\over 2}\nonumber \\ & = y(t_n) + f(t_n,y(t_n))h_n + y''(\tau){h_n^2\over 2}, \label{eq2} \end{align} onde $\tau$ é algum valor entre $t_n$ e $t_{n+1},$ $h_n \equiv (t_{n+1}-t_n)$ e usamos a equação diferencial \eqref{standard} para identificar $y'(t_n)$ como $f(t_n,y(t_n)).$ A fórmula de iteração é obtida da identidade \eqref{eq2}, quando $y(t_n)$ é substituída por sua aproximação $y_n$ e o termo de segunda ordem é desprezado. Ficamos assim com \begin{equation}\label{euler} y_{n+1} = y_n + h_nf(t_n,y_n). \end{equation} Repare que, se em \eqref{eq2} temos uma identidade, em \eqref{euler} o sinal de igual é usado no sentido de que estamos definindo o valor de $y_{n+1}$ desta forma.

O valor exato de $y$ em $t_0$ é conhecido, uma vez que a condição inicial em \eqref{standard} o prescreve como $y_0.$ Desta forma, $r(t) = y_0+(t-t_0)f(t_0,y_0)$ é exatamente a reta tangente à função $y$ em $t_0.$ Portanto, a aproximação $y_1$ é obtida pelo valor desta reta tangente em $t=t_1,$ como esquematizado no gráfico a seguir.

Do ponto de vista prático é muito simples utilizar o método de Euler. Basicamente, é necessário apenas definir a malha. No caso de malhas regulares, basta definir o passo de discretização $h.$ Por exemplo, considere o PVI $$ y' = -2ty, \quad t\gt 0, \qquad y(0)=1. $$ A solução deste PVI é $y(t) = e^{-t^2}.$ Considere uma malha regular para o intervalo $[0,1]$ com passo de discretização $h = 0.1.$ No Octave, as aproximações $y_n,$ para $n=0,1,\ldots,10$ são calculadas com:

f = @(t,y) -2*t*y;      # função da equação diferencial
N = 11;                 # quantidade de pontos na malha
t = linspace(0,1,N)';   # pontos da malha
h = t(2)-t(1)           # passo de discretização
h =  0.10000
y = zeros(N,1);         # vetor para as aproximações
y(1) = 1;               # condição inicial

# Iteração de Euler
for n=1:N-1; y(n+1) = y(n) + h*f(t(n),y(n)); endfor
>
Y = @(t) exp(-t.^2);    # solução exata apenas para conferência
[t y abs(Y(t)-y)]
ans =

   0.00000   1.00000   0.00000
   0.10000   1.00000   0.00995
   0.20000   0.98000   0.01921
   0.30000   0.94080   0.02687
   0.40000   0.88435   0.03221
   0.50000   0.81360   0.03480
   0.60000   0.73224   0.03457
   0.70000   0.64437   0.03175
   0.80000   0.55416   0.02687
   0.90000   0.46550   0.02064
   1.00000   0.38171   0.01383

norm(Y(t)-y,inf)
ans =  0.034803

Repare que no Octave os vetores são indexados a partir de 1. Como pode ser visto, na terceira coluna acima, a diferença entre $y_n$ e $y(t_n)$ foi grande. Veja que a $\max_n|y(t_n)-y_n| = 0.034803.$ Isto era esperado, uma vez que o método de Euler é um método de primeira ordem. Para reduzir este erro, só temos como opção reduzir passo de discretização. Escolhendo $h=0.005,$ e portanto $N=201$, temos que o maior erro absoluto entre $y(t_n)$ e $y_n$ cai para $0.0016325.$

Na próxima aula veremos como construir métodos melhores que o método de Euler, de modo que para obter uma aproximação mais acurada não seja necessário reduzir tão drasticamente o valor de $h,$ o que resultaria em um aumento proporcional na quantidade de iterações realizadas.

1. Para os problemas de valor inicial abaixo, estime $y(1)$ pelo método de Euler usando (i) $h=0.1$ e (ii) $h=0.001.$ Compute o erro absoluto nas duas aproximações, usando a solução exata apresentada.

1. $y'=t^2-y,$ $y(0)=1.$ [ $y(t) = -e^{-t}+t^2-2t+2$ ]
2. $y'=3y+3t,$ $y(0)=1.$ [ $y(t) = {4\over 3}e^{3t}-t-{1\over 3}$ ]

2. Considere o sistema de equações não lineares $$ \left\{ \begin{array}{l} x'(t) = \sigma (y - x), \\ y'(t) = x (\rho -z ) -y, \\ z'(t) = xy - \beta z, \end{array} \right. $$ onde $\sigma,$ $\rho$ e $\beta$ são parâmetros constantes. Denote por $Y(t) = (x(t),y(t),z(t))^T.$

1. Para $\sigma = 10,$ $\rho=28$ e $\beta=8/3$, estime $Y(t)$, partindo de $Y_0 = (-6,-6,15)^T,$ para $0 \le t \le 100.$
2. Faça um gráfico exibindo a trajetória computada, ou seja, plotando os pontos $Y(t).$ (Veja o comando plot3 do Octave.)
3. Resolva novamente o problema, utilizando porém uma condição inicial perturbada por um fator da ordem de $10^{-5},$ ou seja partindo de $\tilde{Y}_0 = Y_0 + R,$ com $||R|| \le 10^{-5}.$ Plote a trajetória $\tilde{Y}(t)$ computada. A figura que surgiu se parece com a obtida no item anterior?
4. Plote a diferença entre as duas trajetórias, ou seja, plot $Y(t)-\tilde{Y}(t).$ O que você observa?
5. Pesquise sobre o efeito borboleta e comente a relação com o experimento acima.