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