RK4 - SBPC2017 (Octave)

% Código utilizado no minicurso de computação científica na 69 SBPC, UFMG - Belo Horizonte

% 'Entre com o tamanho inicial da populacao (10): N(0)= '
N_0 = 10;

%'Taxa (razao) de crescimento (0.1 equivale a 10%), r= '
r = 0.1;

% taxa de morte (0.05 equivale a 5%) d =  (0.05)
d = 0.05;

% Capacidade de suporte
k = 1000;

% 'Tamanho do passo de integracao (0.00001 - 100), h= '
dt  = 1;

% 'Numero de geracoes (200): '
N_t = 200;

% Numero de passos adaptado ao tamnho do passo para dar o numero desejado de geracoes
N_t = N_t / dt;

% constroi um vetor de 0 a N_t com N_t elementos igualmente espacados
t = linspace(0, (N_t)*dt, N_t);

% constroi um vetor de N com 0 em cada posicao e tamnho N_t+2
N = zeros(N_t, 1);

% primeira posicao e o valor inicial
N(1) = N_0;

% primeiro t e zero
t(1) = 0;

% time step
h = dt;

% modelo exponencial
fpop = @(t,x) r*x;

% modelo com morte
% fpop = @(t,x) r*x - (x*d);

% modelo logistico
% fpop = @(t,x) (r*x)*(k-x)/k;

% laco do runge-kutta RK4 de 1 a N_t
for n = 1:N_t

    % calcula os 4 coeficientes
    k1 = h*fpop(t(n),N(n));
    k2 = h*fpop(t(n)+.5*h,N(n)+.5*k1);
    k3 = h*fpop(t(n)+.5*h,N(n)+.5*k2);
    k4 = h*fpop(t(n)+h,N(n)+.5*k3);

    % aplica a formula de RK
    N(n+1) = N(n) + (k1 + 2*k2 + 2*k3 +k4)/6;

    % calcula o proximo tempo
    t(n+1) = t(n) + h;
end

plot(t, N);
xlabel('t'); ylabel('N(t)');

legend('runge-kutta');
filestem = strcat('rk4_', num2str(N_t), '_', num2str(h));
print(filestem, '-dpng');