I have a custom implementation of some ODE solvers for a class. I'm having a problem where if I make the time step dt smaller than .2, the program will stall. However, if I comment out one of the Runge-Kutta solvers, it will execute very quickly, and I can switch which one is commented out so I can get the solution from both solvers. I'm wondering how to fix this. I've been trying to find some way that one solver might be interfering with the other, but I don't see how this could be happening.
Implementation:
global dt;
dt = .5; % going below ~.25 makes the program take a very long time to exit
g = 1;
c_d = 2;
m = 3;
tf = 15;
dudt = @(t, u) g - (c_d/m) * u.^2
[t_euler, u_euler] = euler(dudt, [0, tf], 0);
[t_rk4, u_rk4] = rk4(dudt, [0, tf], 0); % either of this one or rk2
% can be commented out to make the program
% run quickly, but rk2 and rk4 cannot be
% run at the same time
[t_rk2, u_rk2] = rk2(dudt, [0, tf], 0);
%% rk4.m %%
function [t, u] = rk4(odefun, tspan, u0)
t0 = tspan(1);
t = [ t0 ];
t_new = t0;
global dt;
tf = tspan(2);
u_new = u0;
u = [ u0 ];
while (t_new < tf)
if (t_new + dt > tf)
dt = tf - t_new;
end
k1 = dt * odefun(t_new, u_new);
k2 = dt * odefun(t_new + dt/2, u_new + k1/2);
k3 = dt * odefun(t_new + dt/2, u_new + k2/2);
k4 = dt * odefun(t_new + dt, u_new + k3);
u_new = u_new + 1/6 * (k1 + 2*k2 + 2*k3 + k4);
% rk2.m is the same as rk4.m, except u_new = u_new + k2
% euler.m is also the same, except u_new = u_new + k1
u = [ u, u_new ];
t_new = t_new + dt;
t = [ t, t_new ];
end
end
When computing the interval length for the last step to meet tf exactly, you reduce the global variable dt, and in the other method the now smaller value of dt is used. With your data, you might have the new dt as small as 1e-16 or even 0.
You want to treat the global dt as a constant, perhaps use a local variable h instead of dt. Or pass dt as parameter of the functions.