Comparando métodos Runge-Kutta de ordem 2, 3, 4 e 5.

Uso de um passo muito grande para solução de uma EDO, o  método RK-2 apresenta um erro significativo.

Em muitas situações, o método de Runge-Kutta (RK) de 4a. ordem tem precisão suficiente para resolver numericamente uma EDO (equação diferencial ordinária) - ver postagens antigas aqui neste mesmo blog (aqui) - de forma satisfatória. Entretanto, se for necessária uma maior precisão, podemos tentar um método RK de ordem mais elevada. Em casos extremos, um método ordem Runge-Kutta-Nyström ( RKN1210 12a./10a.) pode ser tentado. Nessa postagem, apresentamos uma comparação entre os métodos RK de ordem 2, 3, 4 (forma clássica) e 5 (método sugerido por Butcher - ver aqui a página desse pesquisador aqui). Existem várias fórmulas possíveis para esses métodos, aqui apresentamos as mais comuns. São elas:

RK-2:
\begin{align*} K_1 & = hf(t_k,y(t_k)) \\ K_2 & = hf(t_k + h,y(t_k) + K_1) \\ y(t_{k+1}) & = y(t_k) + \frac{K_1+K_2}{2}\\ t_{k+1} & = t_k + h \end{align*}
RK-3:
\begin{align*} K_1 & = hf(t_k,y_k) \\ K_2 & = hf(t_k + h/2,y_k + K_1/2) \\ K_3 & = hf(t_k + h,y_k + K_2) \\ y_{k+1} & = y_{k} + \frac{1}{4}(K_1 + 2K_2 + K_3)\\ t_{k+1} & = t_k + h \end{align*}
RK-4: \begin{align*} K_1 & = hf(t_k,y_k) \\ K_2 & = hf(t_k + h/2,y_k + K_1/2) \\ K_3 & = hf(t_k + h/2,y_k + K_2/2) \\ K_4 & = hf(t_k + h,y_k + K_3) \\ y_{k+1} & = y_{k} + \frac{1}{6}(K_1 + 2K_2 + 2K_3 + K_4) \\ t_{k+1} & = t_k + h \end{align*}
RK-5:
\begin{align*} K_1 &= h f(t_k,y_k);\\ K_2 &= h f(t_k+h/4,y_k+K_1/4);\\ K_3 &= h f(t_k+h/4,y_k+(K_1+K_2)/8);\\ K_4 &= h f(t_k+h/2,y_k+(K_3 - K_2/2));\\ K_5 &= h f(t_k+3*h/4,y_k+(3 K_1 + 9 K_4)/16); \\ K_6 &= h f(t_k+h,y5kb+(-3 K_1+2 K_2+12 K_3 - 12 K_4 + 8 K_5)/7); \\ y_{k+1} &= y_k + (7 K_1 + 32 K_3 + 12 K_4 + 32 K_5 + 7 K_6)/90; \end{align*}

Iremos testar a EDO:
\begin{equation*} \frac{dy}{dx} = \exp(-x^2/2) - yx; \end{equation*}
cuja solução analítica é \begin{equation*} y = (x+1)\exp(-x^2/2); \end{equation*}
usando os métodos RK descritos acima. Para saber mais sobre os métodos RK para solução de EDOs ver aqui e aqui.

Gráficos:

Comparação da solução analítica com as soluções numéricas. 

Menor erro para os métodos RK de ordem mais alta.

Com um passo $h$ menor, menores os erros de todos os métodos.

Os métodos de ordem mais alta apresentam uma redução do erro ainda melhor que os métodos de baixa ordem.

Erro acumulado em função do passo $h$.

"Zoom" da figura anterior: se o passo for muito pequeno, erros de arrendondamento predominam e diminuem a acurácia.

Código Scilab:

//// Runge-Kutta: comparações
clc; clear; xdel(winsid()); // limpando e fechando janelas

function f=fd(x, y)   //função - EDO
//    f = 9.8 - 0.1*y.*y;
   f = exp(-x.*x/2) - y.*x;
//    f = -y +2*exp(-x).*cos(2*x);
//    f = 2*y./(exp(x)+1) + (2.0)./(exp(2*x)+2*exp(x)+1);
//    f = -2*y + 2*x.*exp(-2*x);
endfunction

function g=ff(x)   //função - solução analítica
    k = 7*sqrt(2)/5;
//    g = sqrt(98)*(1-exp(-k*x))./(1+exp(-k*x));
    g = (x+1).*exp(-x.*x/2);
//    g = exp(-x).*sin(2*x);
//    g = (exp(x)-1)./(exp(x)+1);
//      g = x.*x.*exp(-2*x);
endfunction

h = 0.4; // passo inicial
vh = zeros(10,1); verros = zeros(10,4);
v1 = 1-sqrt(1/2); v2 = 1+sqrt(1/2);
for kh=1:10
    h = h/2; vh(kh) = h;
    fim = round(5/h);
    x0 = 0; y0 = ff(x0);  xk = x0;
    x = zeros(1,fim); x(1) = x0;
    y2 = zeros(1,fim); y2(1) = y0; y2k = y0; // RK-2
    y3 = zeros(1,fim); y3(1) = y0; y3k = y0; // RK-3
    y4 = zeros(1,fim); y4(1) = y0; y4k = y0; // RK-4
    y5b = zeros(1,fim); y5b(1) = y0; y5kb = y0; // RK-5b

for k=2:fim
    // RK2:
    k1 = h*fd(xk,y2k);
    k2 = h*fd(xk+h,y2k+k1);
    y2k = y2k + (k1 + k2)/2;
    y2(k) = y2k;
 
    // RK3:
    k1 = h*fd(xk,y3k);
    k2 = h*fd(xk+h/2,y3k+k1/2);
    k3 = h*fd(xk+h,y3k+k2);
    y3k = y3k + (k1 + 2*k2 + k3)/4;
    y3(k) = y3k;

    // RK4:
    k1 = h*fd(xk,y4k);
    k2 = h*fd(xk+h/2,y4k+k1/2);
    k3 = h*fd(xk+h/2,y4k+k2/2);
    k4 = h*fd(xk+h,y4k+k3);  
    y4k = y4k + (k1 + 2*(k2 + k3) + k4)/6;
    y4(k) = y4k;

    // RK5 - Butcher:
    k1 = h*fd(xk,y5kb);
    k2 = h*fd(xk+h/4,y5kb+k1/4);
    k3 = h*fd(xk+h/4,y5kb+(k1+k2)/8);
    k4 = h*fd(xk+h/2,y5kb+(k3 - k2/2));  
    k5 = h*fd(xk+3*h/4,y5kb+(3*k1 + 9*k4)/16);  
    k6 = h*fd(xk+h,y5kb+(-3*k1+2*k2+12*k3 - 12*k4+8*k5)/7);  
    y5kb = y5kb + (7*k1 + 32*k3 + 12*k4 + 32*k5 + 7*k6)/90;
    y5b(k) = y5kb;

    xk = xk + h;
    x(k) = xk;
end

ya = ff(x);  // função analítica - teórica.
// Erros e somas dos erros absolutos:
e2 = abs(y2-ya); s2 = sum(e2);
e3 = abs(y3-ya); s3 = sum(e3);
e4 = abs(y4-ya); s4 = sum(e4);
e5b = abs(y5b-ya); s5b = sum(e5b);
verros(kh,:) = [s2, s3, s4, s5b];

if kh==1 then  // Gráficos:
    figure; plot(x,ya,x,y2,x,y3,x,y4,x,y5b);
    title(['Passo h: ',string(h)]);    
    xlabel('Eixo - tempo'); ylabel('Função');
    legend('y(t)','RK2','RK3','RK4','RK5b');
    figure; plot(x,e2,x,e3,x,e4*100,x,e5b*1e3);
    legend('RK2','RK3','RK4 x 100','RK5b x 1e3');
    title(['Erro acumulado: ',string([s2, s3, s4, s5b])]);
    xlabel('Eixo - tempo'); ylabel('Erro');
end;
if kh==8 then  // Gráficos:
    figure; plot(x,ya,x,y2,x,y3,x,y4,x,y5b);
    title(['Passo h: ',string(h)]);    
    xlabel('Eixo - tempo'); ylabel('Função');
    legend('y(t)','RK2','RK3','RK4','RK5b');
    figure; plot(x,e2,x,e3,x,e4*1e5,x,e5b*1e6);
    legend('RK2','RK3','RK4 x 1e5','RK5b x 1e6');
    title(['Erro acumulado: ',string([s2, s3, s4, s5b])]);
    xlabel('Eixo - tempo'); ylabel('Erro');
end;
end;

// Erros em função do passo h:
verros = log10(verros);
figure; plot(vh,verros,'-o');
legend('RK2','RK3','RK4','RK5b',4);
title('Erro acumulado'); xlabel('Tamanho do passo h');
ylabel('Erro acumulado - escala log.');