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

Question

This is a post piggy backing off a post made earlier yesterday, I have made considerable progress towards completion but I'm having difficulty with this part, and I'm sure it's last hurdle I need to get by to have my code running smoothly. Once I can see the a graph that is similar it should be no problem working out any other kinks my code may have or adjustments needed for initial conditions, scalars, etc. Using an ODE solver like ODE15 does provide the necessary graphs but this is an exercise that tests whether you can implement the numerical scheme without it. I've been told this question belongs in a coding forum but I disagree due to the amount of numerical analysis that is involved.

The time step is increasing as we move from from 0.0001 to 100,000. So the step size, h, should be increasing in size as it goes along. I would think that I would need an adaptive time step, but don't know how to implement unless I have to use the fancier Runge-Kutta Fehlberg 4/5 pair, but I still have no idea how they came up with the error control conditions and why these values were chosen as necessary. A lot of hand wavy if you ask me, and I would have no idea what values to use for the local error and safety factor since I don't truly understand it conceptually, so to accurately code this in the smoothest way possible, I would prefer to bypass. If I have to tough it out and do it if it's my only option then I would like your input as to why, but before resorting to this way I would like to gather other possible options that are available to me. I did post earlier but was left hanging with regards to this question, I'm trying to do this without using an ODE solver like ODE15 or ODE45. I have h=10^-4 as a dummy variable, just a place holder for now.

This is the ODE system, we take N=5 for simplicity

$\begin{align}\dot c_0=-m(t)c_0+\epsilon c_1\\ \dot c_k=-m(t)c_k-\epsilon c_k+m(t)c_{k-1}+\epsilon c_{k+1}\\ \dot c_N=-\epsilon c_N+m(t)c_{N-1}\end{align}$

Due to constraint condition, $m(t)=M-\sum_{k=1}^5kc_k$

Butcher table

$\begin{array} {c|cccc} 0\\ \frac{1}{2} & \frac{1}{2}\\ \frac{1}{2} &0 &\frac{1}{2} \\ 1& 0& 0& 1\\ \hline & \frac{1}{6} &\frac{1}{3} &\frac{1}{3} &\frac{1}{6} \end{array}$

So, we get

\begin{array}{lr} \vec{Y}_1=\vec{y}_{n-1},\,\\ \vec{Y}_2=\vec{y}_{n-1}+\frac{h}{2}\vec{f}(t_{n-1},\vec{Y}_1),\, \\ \vec{Y}_3=\vec{y}_{n-1}+\frac{h}{2}\vec{f}(t_{n-1}+\frac{h}{2},\vec{Y}_2),\, \\ \vec{Y}_4=\vec{y}_{n-1}+h\vec{f}(t_{n-1}+\frac{h}{2},\vec{Y}_3),\, \\ W4=\vec{y}_{n-1}+h(\frac{1}{6}\vec{f}(t_{n-1},\vec{Y}_1)+\frac{1}{3}\vec{f}(t_{n-1}+\frac{1}{2}h,\vec{Y}_2)+\frac{1}{3}\vec{f}(t_{n-1}+\frac{1}{2}h,\vec{Y}_3)+\frac{1}{6}\vec{f}(t_{n-1}+h,\vec{Y}_4))\end{array}

Here is my code:

function cells

%constants
epsilon=0.0001;
M=30;
N=4;%double check if it plays role in code, mentioned in legend for graph
Ns1=20;
Ns2=10;
m=M;%constraint condition

%computing time variable
tstart=10^-4;
tend=10^5;
t=tstart;
n=500;%total time steps
logstart=log10(tstart);
logend=log10(tend);
lpm=[logstart:(logend-logstart)/n:logend];%log plot mesh
iplot=1;%counter
tarray=zeros(1,n+1);
tarray(iplot)=t;
%h=(10^-4);%%the problem%%

%Numerical scheme for sigma<1/2
c=[Ns1,0,0,0,0];%initial conditions
carray=zeros(5,n+1);
carray(:,iplot)=c;%adding initial condition to array of values

while t<tend %loop for top graph implemented with Classical RK4
    z=RK4(c,h,epsilon,m);
    tm=10^(lpm(iplot));
    if t>tm
        iplot=iplot+1;
        carray(:,iplot)=z;
        tarray(iplot)=t;
    end
    c=z;
    t=t+h;
end

logt=log10(tarray);
figure(1)
tiledlayout(2,1) %initiating display of multiple axes in one figure
nexttile%top plot
plot(logt,carray(2,:),'g')
hold on
plot(logt,carray(3,:),'y')
plot(logt,carray(4,:),'b')
plot(logt,carray(5,:),'k')
hold off

%%Numerical scheme for sigma>1/2
c=[Ns2,0,0,0,0];
carray=zeros(5,n+1);
carray(:,1)=c;

while t<tend%loop for bottom graph implemented with Classical RK4
    z=RK4(c,h,epsilon,m);
    tm=10^(lpm(iplot));
    if t>tm
        iplot=iplot+1;
        carray(:,iplot)=z;
        tarray(iplot)=t;
    end
    c=z;
    t=t+h;
end

logt=log10(tarray);
nexttile%bottom plot
plot(logt,carray(2,:),'g');
hold on
plot(logt,carray(3,:),'y');
plot(logt,carray(4,:),'b');
plot(logt,carray(5,:),'k');
hold off

%Classical RK4
function W4=RK4(c,h,epsilon,m) 
X1=c;
X2=c+h*(1/2)*rhs(X1,epsilon,m);
X3=c+h*(1/2)*rhs(X2,epsilon,m);
X4=c+h*rhs(X3,epsilon,m);
W4=c+h*((1/6)*rhs(X1,epsilon,m)+(1/3)*rhs(X2,epsilon,m)+(1/3)*rhs(X3,epsilon,m)+(1/6)*rhs(X4,epsilon,m));

%System of ODE's
function dcdt=rhs(c,epsilon,m)
dcdt(1)=-m*c(1)+epsilon*c(2);
dcdt(2)=-m*c(2)-epsilon*c(2)+m*c(1)+epsilon*c(3);
dcdt(3)=-m*c(3)-epsilon*c(3)+m*c(2)+epsilon*c(4);
dcdt(4)=-m*c(4)-epsilon*c(4)+m*c(3)+epsilon*c(5);
dcdt(5)=-epsilon*c(5)+m*c(4);

Here is picture of graph I need to display, observe figure description underneath for more info.

Here is link to paper

http://www.csun.edu/~dorsogna/mwebsite/papers/PRE-5tom.pdf

Update: I was able to finally produce the graph but it's not necessarily the same as the one in the paper. Need help debugging or an explanation as to why this is. I believe it might be the step size, h, but I need to be sure. It might be too big. I am only focusing on part a. of the figure for now.

function cells

%constants
epsilon=0.0001; %detachment rate
M=30; %total number of ligands
N=4; %max binding number
Ns1=20; %number of nuclei for first set
Ns2=10; %number of nuclei for second set

%computing time variable
T=500000; %length of integration
tstart=10^-4;
tend=10^5;
t=tstart;
%n=500; %total time steps
%logstart=log10(tstart);
%logend=log10(tend);
%lpm=[logstart:(logend-logstart)/n:logend]; %log plot mesh
iplot=1; %counter
tarray=zeros(1,1);
tarray(iplot)=t;
h=max(10^(-5),0.1*min(1,t));

%Numerical scheme for sigma<1/2
c=[Ns1;0;0;0;0]; %initial conditions
carray=zeros(N+1,1); %removed fixed number of columns for matrix in order to amend column vectors in loop dependent on t 
carray(:,iplot)=c; %adding initial condition to first slot of array 
m=M-sum([1 2 3 4 5]'.*carray(:,iplot));
while t<=tend %loop for top graph implemented with Classical RK4
    z=RK4(c,h,epsilon,m);
    iplot=iplot+1;
    carray=[carray z];
    m=M-sum([1 2 3 4 5]'.*carray(:,iplot));
    t=t+h;
    tarray=[tarray t];
    c=z;
    h=max(10^(-5),0.1*min(1,t));
end

length(tarray)
length(carray(2,:))
logt=log10(tarray);
figure(1)
plot(logt,carray(2,:),'g')
hold on
plot(logt,carray(3,:),'y')
plot(logt,carray(4,:),'b')
plot(logt,carray(5,:),'k')
hold off

%Classical RK4
function W4=RK4(c,h,epsilon,m) 
X1=c;
X2=c+h*(1/2)*rhs(X1,epsilon,m);
X3=c+h*(1/2)*rhs(X2,epsilon,m);
X4=c+h*rhs(X3,epsilon,m);
W4=c+h*((1/6)*rhs(X1,epsilon,m)+(1/3)*rhs(X2,epsilon,m)+(1/3)*rhs(X3,epsilon,m)+(1/6)*rhs(X4,epsilon,m));

%System of ODE's
function dcdt=rhs(c,epsilon,m)
dcdt=zeros(5,1);
dcdt(1)=-m*c(1)+epsilon*c(2);
dcdt(2)=-m*c(2)-epsilon*c(2)+m*c(1)+epsilon*c(3);
dcdt(3)=-m*c(3)-epsilon*c(3)+m*c(2)+epsilon*c(4);
dcdt(4)=-m*c(4)-epsilon*c(4)+m*c(3)+epsilon*c(5);
dcdt(5)=-epsilon*c(5)+m*c(4);

Answer 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)). $$