Usando a Série de Taylor para resolver uma EDO

Brook Taylor (18 agosto de 1685 - 29 dezembro 1731) foi um matemático inglês.

Podemos usar a famosa Série de Taylor (ver aqui) para resolver numericamente uma EDO (equação diferencial ordinária). Sabemos que para uma função $f$ continuamente diferenciável, em $x = a$, podemos escrever:
\[ f(x) = f(a) + f'(x)\frac{x-a}{1!} + f''(x)\frac{(x-a)^2}{2!} + f'''(x)\frac{(x-a)^3}{3!} + ... \]
Para uma EDO na forma $f'(x) = g(x,t)$, com a condição inicial $x = x_0$ em $t = t_0$, podemos escrever:
\[ x_{k+1} \cong x_k + hg(x,t)\] \[ t_{k+1} = t_k + h \] Usando a Série de Taylor: \[ x_{k+1} \cong x_k + hg(x,t) + \frac{h^2 g'(x,t)}{2} + \frac{h^3 g''(x,t)}{3!} + \frac{h^4 g'''(x,t)}{4!} + ... \] sendo o $h$ o passo de integração. Quanto mais termos usarmos, melhor será a precisão alcançada. Esse método, naturalmente, só pode ser usado vantajosamente se o cálculo de $g'(x,t)$, $g''(x,t)$, ..., etc, forem fáceis de calcular. Vejamos um exemplo simples. Seja a EDO \[ \frac{dx}{dt} = e^{-t} - x\] com $x(0) = 0$, cuja solução analítica é $x(t) = te^{-t}$. Então, podemos escrever como solução numérica: \[ x_{k+1} \cong x_k + h(e^{-t_k}-x_k) + h^2(x_k - 2e^{-t_k})/2 + h^3(3e^{-t_k} - x_k)/6 + ...\] O código abaixo apresenta a solução desta EDO usando a Série de Taylor (com 2, 3, 4 e 5 termos) descrita acima e compara com o método de Runge-Kutta de 4a. ordem. Código:

clc; close;

function f=fd(x, y)   //função - EDO
    f = exp(-x) - y;
endfunction

function g=ff(x)   //função - solução analítica
 g = x.*exp(-x);
endfunction

x0 = 0;
y0 = ff(x0);
h = 0.1;
xk = x0;
fim = round(5/h);
vx = zeros(1,fim); vx(1) = x0;
y4= zeros(1,fim); y4(1) = y0; y4k = y0; y = y0; // RK4
yt2 = y0; yt3 = y0; yt4 = 0;
yt2k = zeros(1,fim); yt3k = yt2k; yt4k = yt2k;
    
for k=2:fim
    // RK4:
    k1 = h*fd(xk,y4k);
    k2 = h*fd(xk+h/2,y4k+k1/2);
    k3 = h*fd(xk+h/2,y4k+k2/2);
    k4 = h*fd(xk+h,y4k+k3);  
    y4k = y4k + (k1 + 2*(k2 + k3) + k4)/6;
    y4(k) = y4k;

    ////////////// Taylor:
   yt2 = yt2 + h*(exp(-xk)-yt2)+ h*h*(yt2- 2*exp(-xk))/2; 
   yt3 = yt3 + h*(exp(-xk)-yt3)+ h*h*(yt3- 2*exp(-xk))/2 + h*h*h*(3*exp(-xk)-yt3)/6; 
   yt4 = yt4 + h*(exp(-xk)-yt4)+ h*h*(yt4- 2*exp(-xk))/2 + h*h*h*(3*exp(-xk)-yt4)/6+ h*h*h*h*((yt4-4*exp(-xk)))/24; 
   yt2k(k) = yt2;
   yt3k(k) = yt3;
   yt4k(k) = yt4;
      
    xk = xk + h;
    vx(k) = xk;
end;

ya = ff(vx); close; close; 
plot(vx,y4,vx,yt2k,vx,yt3k,vx,yt4k); title('Soluções analítica e numéricas');
figure; plot(vx, 1000*(ya-y4), vx, ya-yt2k, vx, 10*(ya-yt3k), vx, 100*(ya-yt4k));
legend('Erro RK4 x 1000','Erro T2','Erro T3','Erro T4');
title('Erros');

---

Gráficos: