Comparando os métodos de Runge-Kutta de baixa ordem (RK-2, RK3 (A, B e C) e RK4)

Comparação entre RK2, RK3 e RK4

Os métodos de Runge-Kutta são uma das técnicas numéricas mais populares para resolver Equações Diferenciais Ordinárias (EDOs) da forma:

\( y' = f(t, y), \quad y(t_0) = y_0 \)

Por que usar Runge-Kutta?

Eles permitem obter alta precisão sem calcular derivadas de ordem superior, avaliando a função \(f(t,y)\) em vários pontos dentro de cada passo \(h\).

Comparação dos Métodos

Método	Ordem	Estágios	Erro Local	Erro Global	Custo Computacional
RK2 (Heun, Midpoint...)	2	2	O(h³)	O(h²)	Baixo
RK3	3	3	O(h⁴)	O(h³)	Médio
RK4 Clássico	4	4	O(h⁵)	O(h⁴)	Alto

Fórmulas dos Métodos

1. Runge-Kutta de 2ª Ordem (Método de Heun)

k₁ = f(tₙ, yₙ)
k₂ = f(tₙ + h, yₙ + h k₁)
yₙ₊₁ = yₙ + (h/2)(k₁ + k₂)
    

2. Runge-Kutta de 3ª Ordem (um deles, talvez o mais usado)

k₁ = f(tₙ, yₙ)
k₂ = f(tₙ + h/2, yₙ + (h/2)k₁)
k₃ = f(tₙ + h, yₙ - h k₁ + 2 h k₂)
yₙ₊₁ = yₙ + (h/6)(k₁ + 4k₂ + k₃)
    

3. Runge-Kutta de 4ª Ordem (Clássico)

k₁ = f(tₙ, yₙ)
k₂ = f(tₙ + h/2, yₙ + (h/2)k₁)
k₃ = f(tₙ + h/2, yₙ + (h/2)k₂)
k₄ = f(tₙ + h, yₙ + h k₃)
yₙ₊₁ = yₙ + (h/6)(k₁ + 2k₂ + 2k₃ + k₄)
    

O RK4 continua sendo o método mais usado devido ao excelente equilíbrio entre precisão e simplicidade. O RK2 é útil quando velocidade é prioritária, enquanto o RK3 serve como um meio-termo. O código Scilab a seguir mostra uma comparação prática entre os métodos. São usados 3 RK-3. O valor do passo tem forte influência na acurácia dos métodos. No código Scilab foi incluída a extrapolação de Richardson para melhorar ainda mais o resultado do RK-4. 

Implementação em Scilab:

clc; close(winsid());

function r=ff(t) // função
      // r = exp(-t).*(1 + t.*t);   //y(t)
      r = t.*t.*exp(-t).*cos(3*t);
endfunction

function r=gg(t, y) // derivada
     // r = -y + 2*t.*exp(-t);    // y' = -y + 2texp(-t)
     r = 2*t.*exp(-t).*cos(3*t)-y-3*t.*t.*exp(-t).*sin(3*t);
endfunction

function r=rk2(t, y, h) // RK-2
    k1 = h*gg(t,y);
    k2 = h*gg(t+h,y+k1);
    r = y + (k1+k2)/2;
endfunction

function r=rk3A(t, y, h) // RK-3A
    k1 = h*gg(t,y);
    c2 = 2/3;
    c3 = 2/3;
    a32 = 1;
    a31 = -1/3;
    k2 = h*gg(t+h*c2,y+k1*c2);
    k3 = h*gg(t+h*c3,y+k2*a32+k1*a31);
    r = y + (k1+2*k2+k3)/4;
endfunction
 
function r=rk3B(t, y, h) // RK-3B
    k1 = h*gg(t,y);
    k2 = h*gg(t+h/2,y+k1/2);
    k3 = h*gg(t+h,y-k1+2*k2);
    r = y + (k1+4*k2+k3)/6;
endfunction

function r=rk3C(t, y, h) // RK-3C
    k1 = h*gg(t,y);
    c2 = (15-sqrt(5))/24;
    c3 = (5+sqrt(5))/8;
    a32 = (15+sqrt(5))/11;
    a31 = (-65+3*sqrt(5))/88;
    k2 = h*gg(t+h*c2,y+k1*c2);
    k3 = h*gg(t+h*c3,y+k2*a32+k1*a31);
    r = y + (k1+3*k2+k3)/5;
endfunction

function r=rk4(t, y, h) // RK-4
    k1 = h*gg(t,y);
    k2 = h*gg(t+h/2,y+k1/2);  
    k3 = h*gg(t+h/2,y+k2/2);
    k4 = h*gg(t+h,y+k3);
    r = y + (k1+2*(k2+k3)+k4)/6; // RK-4
endfunction

t = 0; h = 0.05;
y = ff(t);
y2 = y; y3a = y; y3b = y; y3c = y; 
y4 = y; y42 = y; y4m = y;
vy2 = [];
vy3a = [];
vy3b = [];
vy3c = [];
vy4 = [];
vy42 = [];
vy4
while t<5
    vy2 = [vy2, y2];
    vy3a = [vy3a, y3a];
    vy3b = [vy3b, y3b];
    vy3c = [vy3c, y3c];
    vy4 = [vy4, y4];
    vy42 = [vy42, y42];
    
    y2 = rk2(t,y2,h);
    y3a = rk3A(t,y3a,h);
    y3b = rk3B(t,y3b,h);
    y3c = rk3C(t,y3c,h);
    y4 = rk4(t,y4,h);
    y42 = rk4(t,y42,h/2); /// passo h/2
    y42 = rk4(t+h/2,y42,h/2);  /// passo h/2
    t = t + h;
end

vy4m = (vy42*16-vy4)/15; // Extrapolação de Richardson

tempo = 1:max(size(vy4));
tempo = tempo - 1; tempo = tempo*h;

ya = ff(tempo);  // função analítica
title('Função no tempo e soluções numéricas');
plot(tempo,ya,tempo,vy2,tempo,vy3a,tempo,vy3b,tempo,vy4);

// erros:
e2 = abs(vy2-ya);   e2L = log10(e2(2:$));
e3a = abs(vy3a-ya); e3aL = log10(e3a(2:$));
e3b = abs(vy3b-ya); e3bL = log10(e3b(2:$));
e3c = abs(vy3c-ya); e3cL = log10(e3c(2:$));
e4 = abs(vy4-ya);   e4L = log10(e4(2:$));
e4m = abs(vy4m-ya);   e4Lm = log10(e4m(2:$));
figure; 
plot(tempo,e2,tempo,e3a,tempo,e3b,tempo,e3c,tempo,e4,tempo,e4m);
legend('Erro 2','Erro 3A','Erro 3B','Erro 3C','Erro 4','Erro 4 m');
title('Erros, h ='+string(h));

figure; 
tp = tempo(2:$);
plot(tp,e2L,tp,e3aL,tp,e3bL,tp,e3cL,tp,e4L,tp,e4Lm);
title('Erros em escala Log, h ='+string(h));
legend('Erro 2','Erro 3A','Erro 3B','Erro 3C','Erro 4','Erro 4 m',3);

 Resultados gráficos: