Efeito do passo para solução de uma EDO

A solução numérica de uma equação diferencial (ver as postagens anteriores aqui e aqui) requer o uso de um passo de integração. Se esse passo for muito grande, os erros podem ser intoleráveis ou levar a uma instabilidade da solução; se o passo for muito pequeno o custo computacional pode ser muito elevado e tempo gasto também. Além disso, uma redução do passo não leva necessariamente a um erro final menor. Logo, a escolha desse passo é uma tarefa não trivial. Existem métodos para fazer um ajuste automático do passo (solução de passo variável) que conseguem uma boa aproximação numérica dentro de um tempo de cálculo razoável.

Nesta postagem evidenciamos, através de um exemplo numérico, que a escolha do passo pode reduzir o erro entre a solução numérica e a solução analítica a um valor mínimo. Isso ocorre porque o erro é composto por duas parcelas: erros de truncamento (devido à fórmula usada na solução da equação diferencial) e erros de arredondamento (devido o número finito de bits do computador). Quanto menor o valor do passo, menor o erro de truncamento, mas os erros de arredondamento vão se acumulando. Assim, vai existir um ponto (passo ótimo) onde essa soma é mínima (ver figura abaixo). O passo ótimo também depende do método usado para solucionar a EDO (equação diferencial ordinária).

A escolha do passo pode minimizar o erro total da solução numérica.

Exemplo: solucionar a equação diferencial: $$\frac{dy}{dt} = 2te^{-t} - y(t),$$
com $y(0) = 0$. A solução analítica é $y(t) = t^{2}e^{-t}$. A solução analítica será usada para calcular os erros das soluções numéricas. Solução e cálculo dos erros usando o método Runge-Kutta de segunda ordem e o Runge-Kutta de 4a. ordem (código Scilab):

function f=fty(t, y)
    f = 2*t*exp(-t) - y;
endfunction

close; close;  close; close;
clc;

dt = 0.1;
y2 = 0;
y4 = 0;
t = 0;
yt = 0;

ve2 = zeros(1,10);
ve4 = ve2;
vdt = ve2;
vtempo2 = ve2;
vtempo4 = ve2;

for k=0:12
    dt = 1e-7*(10^(k/2));
    t = 0;
    y2 = 0;
    disp(k);
    t1=clock();
    while (t<0.5)
      k1 = fty(t,y2);
      k2 = fty(t+dt,y2+k1*dt);
      y2 = y2 + dt*(k1+k2)/2;
        
      t = t + dt;
    end;
    t2 = clock();
    yt = t*t*exp(-t);
    ve2(k+1) = abs(y2-yt);
    vdt(k+1) = dt;
    vtempo2(k+1) = etime(t2,t1);
end

for k=0:12
    dt = 1e-7*(10^(k/2));
    t = 0;
    y4 = 0;
    disp(k);
    t1=clock();
    while (t<0.5)

      k1 = fty(t,y4);
      k2 = fty(t+dt/2,y4+k1*dt/2);
      k3 = fty(t+dt/2,y4+k2*dt/2);
      k4 = fty(t+dt,y4+k3*dt);        
      y4 = y4 + dt*(k1+2*k2+2*k3+k4)/6;
        
      t = t + dt;
    end;
    t2 = clock();
    yt = t*t*exp(-t);
    ve4(k+1) = abs(y4 - yt);
    vdt(k+1) = dt;
    vtempo4(k+1) = etime(t2,t1);
end

plot(vdt,ve2,'-om',vdt,ve4,'-*b');
legend('Erro RK2','Erro RK4');
xlabel('Tamnaho do passo (dt)'); ylabel('Erro absoluto: abs(yRK - yt)');
a=gca();
a.log_flags = "lln";

figure; 
vtempo2 = vtempo2 + 1e-6;
vtempo4 = vtempo4 + 1e-6;
plot(vdt,vtempo2,'-om',vdt,vtempo4,'-*b'); title('Tempo de execução');
xlabel('Tamnaho do passo (dt)'); ylabel('Tempo (s)');
legend('Tempo para o RK2','Tempo para o RK4');
a=gca();
a.log_flags = "lln";

---------------------------------------------------

Gráficos:

Fica evidente que o Runge-Kutta de 4a. ordem apresenta um desempenho muito superior ao RK-2, e mesmo efetuando um maior números de cálculos por iteração, o erro final pode ser muito menor. O posso ótimo varia com o método usado, sendo $10^{-3}$ para o RK-4. Como esperado, o tempo de execução é maior para o RK-4 para um mesmo número de iterações. Fica também evidente que o custo total (tempo computacional mais erros) leva à escolha de um passo "moderado", um valor menor ou igual a $0,01$ para o RK-4 já deve ser satisfatório.