.jpg)
Podemos resolver uma EDO (equação diferencial ordinária) usando métodos numéricos. Um dos métodos mais 'populares' são os da família Runge-Kutta. Quanto mais alta a ordem do método usado, melhor a precisão encontrada (menores os erros). Os métodos de ordem mais alta (5a., 6a, e 7a., ...) em geral, não são mencionados nos livros didáticos.
No código Scilab abaixo resolvemos a EDO
\[ \frac{dy}{dt} = 3e^{-t}\cos(3t)-y(t) \]
usando os métodos RK-3, RK-4, RK-5, RK-6 e RK-7.
Código:
close(winsid()); clc;
// Definir a função da equação diferencial dy/dt = f(t, y)
function dydt=ff(t, y)
//dydt = exp(-t) - y;
dydt = 3*exp(-t).*cos(3*t) - y;
endfunction
// Implementação do método Runge-Kutta de 6ª ordem
function [t, y3, y4, y5, y6, y7]=runge_kutta(t0, y0, tf, h)
// t0: tempo inicial
// y0: valor inicial de y
// tf: tempo final
// h: tamanho do passo
// Número de passos
N = ceil((tf - t0) / h);
// Inicializar vetores para armazenar os resultados
t = t0:h:(t0 + N*h);
y = zeros(1, N+1);
y4 = y;
y6 = y;
y5 = y;
y3 = y;
y7 = y;
y6(1) = y0;
y4(1) = y0;
y5(1) = y0;
y3(1) = y0;
y7(1) = y0;
// Coeficientes do método RK6 (Butcher tableau simplificado)
for i = 1:N
ti = t(i);
y6i = y6(i);
y4i = y4(i);
y5i = y5(i);
y3i = y3(i);
y7i = y7(i);
// Estágios do RK3:
k1 = h * ff(ti, y3i);
k2 = h * ff(ti + h/2, y3i + k1/2);
k3 = h * ff(ti + h, y3i + (2*k2-k1));
// Atualizar y usando a combinação dos estágios
y3(i+1) = y3i + (k1 + 4*k2 + k3)/6;
// Estágios do RK4:
k1 = h * ff(ti, y4i);
k2 = h * ff(ti + h/2, y4i + k1/2);
k3 = h * ff(ti + h/2, y4i + k2/2);
k4 = h * ff(ti + h, y4i + k3);
// Atualizar y usando a combinação dos estágios
y4(i+1) = y4i + (k1 + 2*k2 + 2*k3 + k4)/6;
// Estágios do RK5 (xx estágios)
k1 = h * ff(ti, y5i);
k2 = h * ff(ti + h/4, y5i + k1/4);
k3 = h * ff(ti + 3*h/8, y5i + 3*k1/32 + 9*k2/32);
k4 = h * ff(ti + 12*h/13, y5i + 1932*k1/2197 - ...
7200*k2/2197 + 7296*k3/2197);
k5 = h * ff(ti + h, y5i + 439*k1/216 - 8*k2 + ...
3680*k3/513 - 845*k4/4104);
k6 = h * ff(ti + h/2, y5i - 8*k1/27 + 2*k2 - ...
3544*k3/2565 + 1859*k4/4104 - 11*k5/40);
// Atualizar y usando a combinação dos estágios (fórmula de 5ª ordem)
y5(i+1) = y5i + (16*k1/135 + 6656*k3/12825 + 28561*k4/56430 - ...
9*k5/50 + 2*k6/55);
// Estágios do RK6
k1 = h * ff(ti, y6i);
k2 = h * ff(ti + h/5, y6i + k1/5);
k3 = h * ff(ti + 2*h/5, y6i + 2*k2/5);
k4 = h * ff(ti + h/5, y6i + (11*k1 - 4*k2+5*k3)/60);
k5 = h * ff(ti + 3*h/5, y6i + (-3*k1+4*k2+3*k3+8*k4)/20);
k6 = h * ff(ti + 4*h/5, y6i + (-12*k2-8*k3+22*k4+10*k5)/15);
k7 = h * ff(ti + h, y6i + (51*k1+140*k2+225*k3-280*k4-...
120*k5+60*k6)/76);
// Atualizar y usando a combinação dos estágios
y6(i+1) = y6i + (19*k1+50*k3+75*k4+50*k5+75*k6+19*k7)/288;
// // Estágios do RK7 (xx estágios)
k1 = h * ff(ti, y7i);
k2 = h * ff(ti + h/6, y7i + k1/6);
k3 = h * ff(ti + h/3, y7i + k2/3);
k4 = h * ff(ti + h/2, y7i + (k1+3*k3)/8);
k5 = h * ff(ti + 2*h/11, y7i + (148*k1+150*k3-56*k4)/1331);
k6 = h * ff(ti + 2*h/3, y7i -404*k1/243-170*k3/27 + ...
4024*k4/1701 + 10648*k5/1701);
k7 = h * ff(ti + 6*h/7, y7i+2466*k1/2401+1242*k3/343-...
19176*k4/16807-51909*k5/16807+1053*k6/2401);
b8=77/1440;
k8 = h * ff(ti, y7i+(k1/(576)+k4/(105)-1331*k5/(279552)-...
9*k6/(1024)+343*k7/(149760))/b8);
k9 = h * ff(ti + h, y7i+(-71/32-270*b8/11)*k1-195*k3/22+...
32*k4/7+29403*k5/3584-729*k6/512+1029*k7/1408+270*b8*k8/11);
// Atualizar y usando a combinação dos estágios (fórmula de 7ª ordem)
y7(i+1) = y7i+(77/1440 - b8)*k1+32*k4/105+1771561*k5/6289920+...
243*k6/2560+16807*k7/74880 + b8*k8+11*k9/270;
end
endfunction
// Parâmetros do problema
t0 = 0; // Tempo inicial
y0 = 0; // Condição inicial y(0) = 1
tf = 5; // Tempo final
h = 0.1; // Tamanho do passo
// Resolver a equação
[t, y3, y4, y5, y6, y7] = runge_kutta(t0, y0, tf, h);
// Solução analítica para comparação
//y_analitica = t.*exp(-t);
// eqs = 'dy/dt = exp(-t) - y';
y_analitica = exp(-t).*sin(3*t);
eqs = 'dy/dt = 3*exp(-t).*cos(3*t) - y';
e3 = y_analitica - y3;
e4 = y_analitica - y4;
e5 = y_analitica - y5;
e6 = y_analitica - y6;
e7 = y_analitica - y7;
// Exibir os resultados
disp("t y_RK4 y_RK6 y_analítico");
for i = 1:2:round(length(t)/2)
printf("%.2f %.4f %.4f %.4f\n",...
t(i), y4(i), y6(i), y_analitica(i));
end
// Plotar os resultados
scf(); // Nova janela gráfica
subplot(2,1,1);
plot(t, y_analitica, "r", "LineWidth", 2);
plot(t, y3, "g", "LineWidth", 2);
plot(t, y4, "k", "LineWidth", 2);
plot(t, y5, "b", "LineWidth", 2);
plot(t, y6, "y", "LineWidth", 2);
plot(t, y7, "m", "LineWidth", 2);
xlabel("tempo"); ylabel("y(t)");
title("Solução da equação "+eqs+" usando Runge-Kutta.");
legend(["Analítico",'Num. RK3','Num. RK4',...
'Num. RK5',"Num. RK6",'Num. RK7']);
subplot(2,5,6); plot(t,e3, "LineWidth", 2);
title('Erro - RK3');
subplot(2,5,7); plot(t,e4, "LineWidth", 2);
title('Erro - RK4');
subplot(2,5,8); plot(t,e5, "LineWidth", 2);
title('Erro - RK5');
subplot(2,5,9); plot(t,e6, "LineWidth", 2);
title('Erro - RK6');
subplot(2,5,10); plot(t,e7, "LineWidth", 2);
title('Erro - RK7');
/// Erros acumulados:
erros = [sum(abs(e3)), sum(abs(e4)), sum(abs(e5)), ...
sum(abs(e6)), sum(abs(e7))];
disp(erros);Resultado gráfico:
********
Comparação:
Erros acumulados (RK3, RK4, RK4, RK5, RK6 e RK7, passo h = 0.1):
0.0014108 0.0000166 0.0000008 9.048D-09 4.106D-11
Percebemos que o erro diminui com o aumento da ordem do método RK, aproximadamente duas ordens de grandeza de um $RK_i$ para o $RK_{i+1}$.
Obs: a complexidade computacional desses métodos cresce de forma linear, mas a dedução dos parâmetros desses métodos é bem mais complicada quando a ordem aumenta.
Bibliografia:
SÉKA, Hippolyte; KOUASSI, Assui Richard. A New Seventh Order Runge-kutta
Family: Comparison with the Method of Butcher and Presentation of a
Calculation Software. Mathematics and Computer Science, v. 4, n. 3, p. 68-75, 2019.
DE AQUINO, Francisco José Alves. Solução numérica de equações
diferenciais ordinárias usando Runge-Kutta: um estudo comparativo com
Scilab. CQD-Revista Eletrônica Paulista de Matemática, 2018.
.jpg)