Hi Stacy.
Basically, you seem to be willing to solve your system of ODE by using the method of forward Euler. If so, there is an error in your fomula (see % ***). The forward Euler method writes:
It seems that you programmed:
Whith this modification (and this time step h) both methods provide quite similar results.
function ode()
global alpha_0 alpha beta n
alpha_0 = 0.03;
alpha = 298.2;
beta = 0.2;
n = 2;
h = 0.01; % = time step
Tmax = 10; % = max. time
time = 0:h:Tmax;
init_cond = [10 10 0 0 0 0];
% Solve by using forward Euler
x = zeros(Tmax,6);
x(1,1) = 10;
x(1,2) = 10;
for t = 1:numel(time)-1
x(t+1, 1) = x(t, 1) + h * (beta*(x(t,4) - x(t,1))); % ***
x(t+1, 2) = x(t, 2) + h * (beta*(x(t, 5) - x(t,2))); % ***
x(t+1, 3) = x(t, 3) + h * (beta*(x(t,6) - x(t,3))); % ***
x(t+1, 4) = x(t, 4) + h * (alpha_0 + (alpha/(1+(x(t,3))^n)) - (x(t,4))); % ***
x(t+1, 5) = x(t, 5) + h * (alpha_0 + (alpha/(1+(x(t,1)^n))) - (x(t,5))); % ***
x(t+1, 6) = x(t, 6) + h * (alpha_0 + (alpha/(1+(x(t,2))^n))-(x(t,6))); % ***
end
plot (time,x(:,1),'--r',time,x(:,2),':b',time,x(:,3),'-.k');
hold on
% Solve by using Runge-Kutta
t = (0:1:10);
[t, output] = ode45(@repressilator, t, init_cond);
plot(t, output(:,1),'ro', t, output(:,2),'bo', t, output(:,3),'ko');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d = repressilator(t, ic)
global alpha_0 alpha beta n
pA = ic(1);
pB = ic(2);
pC = ic(3);
mA = ic(4);
mB = ic(5);
mC = ic(6);
d = zeros(6, 1);
d(1) = beta*mA - beta*pA;
d(2) = beta*mB - beta*pB;
d(3) = beta*mC - beta*pC;
d(4) = alpha_0 + (alpha/(1+(pC)^n)) - mA;
d(5) = alpha_0 + (alpha/(1+(pA)^n)) - mB;
d(6) = alpha_0 + (alpha/(1+(pB)^n)) - mC;
end
For an unknown reason I cannot run directly in the web browser. On my computer I get this figure (points are computed by using ode45, dashed lines are computed by using forward Euler):
Mind that this depends on the time step: by taking a higher time step the results are much different. That's normal because ode45 implements a couple of Runge-Kutta formulas of order 4 and 5 (= a couple of formulas with different accuracy are used together to estimate the error, so as to adapt the time step), and therefore it is highly accurate. Conversely, the forward Euler method (which in this implementation has a fixed time step) converges only at the first order.
Both methods converges, so they are bound to converge to the same result as the time step is reduced.
In principle, both methods are only conditionnally stable, that is they diverges for long Tmax, depending on the time step. In the case of forward Euler, the method will diverge if the time step h is higher that the double of the smallest time constant of the system you want to simulate. In the case of ode45, the fact that time step is adaptative can "save" the situation. However, in the case of stiff ode (= in shorts, ode in which you have very different time constants) the number of required step will be very high, and so the computational time. Luckly, this doesn't seem to be your case.
