Resolvendo uma EDO usando Runge-Kutta de 3a., 4a., 5a., 6a. e 7a. ordem - comparação.

 

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.