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:
