Resolvendo numericamente uma equação diferencial de terceira ordem

Em uma postagem anterior (ver aqui) mostramos como resolver numericamente uma equação diferencial ordinária (EDO) de 2a. ordem. Agora, iremos mostrar detalhadamente como resolver uma de 3a. ordem. A ideia é basicamente uma extensão do problema de 2a. ordem, mas existe um detalhe importante: técnicas menos precisas que o método Runge-Kutta 4 (RG4) podem ser instáveis.Uma equação diferencial de primeira ordem pode ser escrita como:
$$\frac{dy}{dt} = f(t,y)$$
A solução de Euler (mais simples) é dada por:
$$\begin{align*} y_{n+1} & = y_n + \Delta t f(t_n,y_n) \\  t_{n+1} & = t_n + \Delta t \end{align*}$$
Sendo que $y_0$ deve ser conhecido para um instante $t_0$ e passo $\Delta t$ não pode ser muito grande. Se $\Delta t$ for um valor muito grande o algoritmo terá problemas de estabilidade. A solução pelo método RG4 é dada pelo seguinte algoritmo:
$$\begin{align*} k_1 &= f(t_n, y_n)\\ k_2 &= f(t_n + \Delta t/2, y_n + k_1 \Delta t/2)\\  k_3 &= f(t_n + \Delta t/2, y_n + k_2 \Delta t /2)\\  k_4 &= f(t_n + \Delta t, y_n + k_3 \Delta t)\\  y_{n+1} &= y_n + \Delta t(k_1 + 2k_2 + 2k_3 + k_4)/6\\  t_{n+1} &= t_n + \Delta t \\ \end{align*}$$
Esse algoritmo leva a resultados muito melhores que o método de Euler. Já uma equação de 2a. ordem é da forma:
$$\frac{d^2y}{dt^2} + f(t,y) \frac{dy}{dt} = g(t,y)$$
E sua solução requer o uso de uma variável auxiliar:
  $$z(t,y) = \frac{dy}{dt}$$
Assim, teremos que resolver simultaneamente
$$\frac{dy}{dt} = z(t,y)$$ e
$$\begin{align*} \frac{dz}{dt} = & g(t,y) - f(t,y) z(t,y) \\ = & h(t,y,z) \end{align*}$$  
Já uma equação de 3a. ordem é da forma:
$$\frac{d^3y}{dt^3} + f(t,y) \frac{d^2y}{dt^2} + g(t,y) \frac{dy}{dt} = h(t,y)$$
E sua solução requer o uso de duas variáveis auxiliares:
  $$z(t,y) = \frac{dy}{dt}, \text{ } w(t,y,z) = \frac{dz}{dt}$$
Assim, teremos que resolver simultaneamente
$$\begin{align*} \frac{dy}{dt} = & z(t,y) \\ \frac{dz}{dt} = & w(t,y,z) \\ \text{e } &\\ \frac{dw}{dt} = & h(t,y) - g(t,y)z - f(t,y)w \\ = & s(t,y,z,w)\\ \end{align*}$$  
Um exemplo. Seja a EDO descrita por:
$$\frac{d^3y}{dt^3} + \frac{1}{3}\frac{d^2y}{dt^2} + \frac{dy}{dt} = \frac{1}{3}y + \frac{1}{9}exp(-t/3)\cos(t)$$
 Com as seguintes condições inicias: $y(0) = 0$, $z = dy/dt = 1$, $w = dz/dt = -2/3$, $t_0 = 0$, $\Delta t = 0.05$. Então
$$\begin{align*} \frac{dy}{dt} = & z\\ \frac{dz}{dt} = & w\\ \frac{dw}{dt} = & -w/3 - z + y/3 + (1/9)exp(-t/3)cos(t) \end{align*}$$
 Logo, o algoritmo necessário é dado por:
$$\begin{align*} ky1 = &\Delta t*z; kz1 = w*\Delta t; kw1 = \Delta t*ftyzw(t,y,z,w);\\ ky2 = &\Delta t*(z + kz1/2); kz2 = \Delta t*(w + kw1/2); \\ kw2 = &\Delta t*(ftyzw(t+dt/2,y+ky1/2,z+kz1/2,w+kw1/2));\\ ky3 = &\Delta t*(z + kz2/2); kz3 = \Delta t*(w + kw2/2); \\ kw3 = &\Delta t*(ftyzw(t+dt/2,y+ky2/2,z+kz2/2,w+kw2/2));\\ ky4 = &\Delta t*(z + kz3); kz4 = \Delta t*(w + kw3); \\ kw4 = &\Delta t*(ftyzw(t+dt,y+ky3,z+kz3,w+kw3));\\ y = &y + (ky1 + 2*ky2 + 2*ky3 + ky4)/6; \\ z = &z + (kz1 + 2*kz2 + 2*kz3 + kz4)/6; \\ w = &w + (kw1 + 2*kw2 + 2*kw3 + kw4)/6;\\ t = &t + \Delta t; \end{align*}$$
com $ftyzw(t,y,z,w) = -w/3 - z + y/3 + (1/9)exp(-t/3)cos(t)$. Então, sua solução numérica usando RG4 é dada pelo código Scilab:
ky1 = dt*z; kz1 = w*dt; kw1 = dt*ftyzw(t,y,z,w);
ky2 = dt*(z + kz1/2); kz2 = dt*(w + kw1/2);
kw2 = dt*(ftyzw(t+dt/2,y+ky1/2,z+kz1/2,w+kw1/2));
ky3 = dt*(z + kz2/2); kz3 = dt*(w + kw2/2);
kw3 = dt*(ftyzw(t+dt/2,y+ky2/2,z+kz2/2,w+kw2/2));
ky4 = dt*(z + kz3); kz4 = dt*(w + kw3);
kw4 = dt*(ftyzw(t+dt,y+ky3,z+kz3,w+kw3));
y = y + (ky1 + 2*ky2 + 2*ky3 + ky4)/6;
z = z + (kz1 + 2*kz2 + 2*kz3 + kz4)/6;
w = w + (kw1 + 2*kw2 + 2*kw3 + kw4)/6;
t = t + dt;
onde ftywz é uma função descrita no código.

Erro de Aproximação
Pro sua própria natureza, os métodos numéricos resultam em soluções aproximadas. Um certo grau de erro ocorre em casa iteração e esses tendem a se acumular. Em alguns casos, após várias iterações, a distância da solução real poderá ser tão grande que novos cálculos se tornam inúteis. Os erros cometidos são de duas naturezas: erros de arredondamentos e erros de truncamento.

Erros de Arredondamentos:
Quando usamos um computador ou uma calculadora para aproximar a solução de uma equação diferencial, por apresentarem casas decimais finitas trabalha-se com, erro de arredondamento.

Erros de truncamento:
Erros de truncamentos são causados pelo tipo de técnica empregada para a atualização do valor de $y$. Os erros de truncamento de Euler são maiores que os de Runge-Kutta.

Sobre Runge e Kutta

"Carl David Runge (1856-1927), matemático e físico alemão, trabalhou muitos anos em espectroscopia. A análise de dados o levou a considerar problemas em computação numérica e o método de Runge-Kutta tem origem em seu artigo sobre soluções em 1901 por M. Wilhelm Kutta (1867-1944). Kutta era um matemático alemão que trabalhava com aerodinâmica e é, também, muito conhecido por suas contribuições importantes à teoria clássica de aerofólio." (Valle, 2012)

Bibliografia

BOYCE, E. William, DIPRIMA, C. Richard. Equações diferenciais elementares e problemas de valores de contorno. Sétima edição, editora LTC.
BARROSO, Leonidas Conceição; et al. Calculo numérico (com aplicações). 2. Ed. Minas Gerais, editora HARBRA.
VALLE, Karine Nayara F. Métodos Numéricos de Euler e Runge-Kutta. Monografia (Especialização) - Programa de Pós-graduação em Matemática para Professores com Ênfase em Cálculo da UFMG. 2012.

Código Scilab:

function ydot=ftyzw(t, y, z, w)
    ydot = -w/3 - z + y/3 + (1/9)*exp(-t/3).*cos(t);
    //ydot = -6*y + cos(4*t)
endfunction

clc;
close;

/////////////////////
// y = exp(-t/3)sin(t)
// y''' + (1/3)y'' + y' = (1/3)y + (1/9)exp(-t/3)cos(t)
// y' = z
// z' = w
// w' = f(t,y,z,w)
////////////////////

y = 0; z = 1; w = -2/3;
t = 0; dt = 0.05;
vy = zeros(1,200);
for k=1:200
    vy(k)=y;
    ky1 = dt*z; kz1 = w*dt;
    kw1 = dt*ftyzw(t,y,z,w);
    ky2 = dt*(z + kz1/2);
    kz2 = dt*(w + kw1/2);
    kw2 = dt*(ftyzw(t+dt/2,y+ky1/2,z+kz1/2,w+kw1/2));
    ky3 = dt*(z + kz2/2);
    kz3 = dt*(w + kw2/2);
    kw3 = dt*(ftyzw(t+dt/2,y+ky2/2,z+kz2/2,w+kw2/2));
    ky4 = dt*(z + kz3);
    kz4 = dt*(w + kw3);
    kw4 = dt*(ftyzw(t+dt,y+ky3,z+kz3,w+kw3));
    y = y + (ky1 + 2*ky2 + 2*ky3 + ky4)/6;
    z = z + (kz1 + 2*kz2 + 2*kz3 + kz4)/6;
    w = w + (kw1 + 2*kw2 + 2*kw3 + kw4)/6;
    t = t + dt;
end
ty=0:0.05:(k*dt-0.05);
yy = exp(-ty/3).*sin(ty);
erro3 = abs(yy-vy);
figure; plot(ty,vy,ty,yy,ty,1e6*erro3);
legend('Solução analítica','Solução numérica','|Erro x 10e6|');
title('Método de Runge-Kutta 4 - EDO de 3a. ordem');
xlabel('tempo (s)'); ylabel('Amplitude');