Finding appropriate step size for ODE system numerically solved without ODE solver in Matlab

Solution 1:

Reading the paper, $c_k(t)$ is the number (averaged, per some unit volume, more of a continuous density with particle number interpretation) of particles or nuclei with $k$ bound ligands, a maximum of $N$ binding places are available per nucleus, there is a number of $N_s$ nuclei that stays constant, and a total amount of $M$ ligands in the whole system, staying also constant.

In a first consequence, at time $t$ there are $$ m(t)=M-\sum_{k=1}^Nkc_k(t) $$ free ligands available to bind to one of the free binding places.

While the dynamic starts out as $$ \dot c_k(t)=v_k(t)-v_{k-1}(t),~~~~ \left\{\begin{aligned} v_{-1}(t)&=v_N(t)=0,~\text{ else }\\v_k(t)&=-p_km(t)c_k(t)+q_kc_{k+1}(t) \end{aligned}\right.$$ for the numerical experiments it gets simplified to $p_k=p$, $q_k=q$, and then further by rescaling the time to $p=1$ and $q=\varepsilon$, $$ v_k(t)=-m(t)c_k(t)+εc_{k+1}(t) $$

So the model can be implemented as

function dotc = model(t,c,M,eps)
    N = length(c)-1;
    m = M - sum([0:N]'.*c);
    v = -m*c;
    v(1:N) += eps*c(2:N+1);
    v(N+1)=0;
    % the v formula is now complete, from now on dotc gets computed in the place of v
    v(2:N+1)-=v(1:N)
    dotc = v 
end

This is then called with the parameters and built-in solver (octave version)

  % Parameters
  Ns = 20; % number of nuclei a) 20, b) 10
  N = 4; % max binding number
  M = 30; % total number of ligands
  eps = 1e-4; % detachment rate
  T=5e5; % length of integration
  
  % initial state
  c_init = zeros(N+1,1);
  c_init(1) = Ns;
  
  % solver parameters  
  opts = odeset('AbsTol',1e-5,'RelTol',1e-6);
  % call solver in step-reporting mode 
  [t_dyn,c_dyn] = ode15s(@(t,c)model(t,c,M,eps),[0,T],c_init,opts);
  
  % plots
  clf;
  subplot(2,1,1)
  plot(log10(1e-4+t_dyn),c_dyn(:,2:end));
  subplot(2,1,2)
  semilogy(log10(1e-4+t_dyn(1:end-1)),t_dyn(2:end)-t_dyn(1:end-1));

resulting in the plots

$N_s=20$	$N_s = 10 $

The upper plots looks like the ones of the paper. In the step size plots a wide variation can be seen, from $10^{-5}$ to $10^4$. This is only possible with an implicit solver like ode15s using implicit multi-step methods for stiff systems.

The explicit RK solver ode45 stabilizes the step size around $h=0.2$, with the corresponding time taken to reach $T=10^5$ due to an excessive number of steps. Larger step sizes are not possible, as they would lead outside the stability region and induce explosive oscillations.

So I'm sure that the attempt with RK4 will be painfully slow. From the step size behavior I would try to steer the step size as $$ h = \max(10^{-5},0.1·\min(1,t)). $$