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

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. 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);


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

• Dude, you did too much but thank you. It's very insightful. I've been rereading the article over and over again, I'm still picking up something new every time I read it. The learning curve for understanding these peer reviewed articles is amazing. I guess for a student to grasp all the material in a research paper usually doesn't happen in a first outing. Dec 7, 2021 at 14:04
• How did you graph the step sizes in the bottom two graphs? That's very interesting, I was thinking the graph would be smoother and continuous. But thank you for providing a step size that I can use, but I am considering using an implicit multistep method now since you mentioned the differences. I would assume the solution would be more robust Dec 7, 2021 at 18:50
• The lower plots are the semilogy plots, the step size is constructed using the differences of the time points, the returned points are the internally selected points for the optimal step size given the error level. Dec 7, 2021 at 19:53
• Your method of implementing the model and figuring out a way to incorporate m(t) without having to resort in making it a constant like I did really gives me a sigh of relief. So I was right about my initial inclination to code m(t) within the framework of the ODE system. I was hoping that setting m(t)=M would satisfy things since it was proving too difficult to code the former but I see that's not going to work. Dec 7, 2021 at 23:19
• To fully understand what you did, mind explaining how you got T, the length of integration, and what it's used for? I see you passing on small t to the model function but was there an initial value that it was equal to before doing this? And what role does it have after you pass it on? I don't see it being used once this happens Dec 7, 2021 at 23:23