ODE solver vs discrete integration

I'm wondering if anyone can help me out. And my confusion may be because my understanding of the math is weak and not necessarily a MATLAB issue. I have a set of ODEs and I'm trying to solve them using ODE45 and discretely and I'm getting different results. This is what I have:

alpha_0 = 0.03;
alpha = 298.2;
beta = 0.2;
n = 2;
time = [0:300];
method = 1;
init_cond = [10 10 0 0 0 0];
t = (0:1:10);

if method == 1
    x = zeros(Tmax,6);
    x(1,1) = 10;
    x(1,2) = 10;

    for t = 1:300
        x(t+1, 1) = beta*(x(t,4) - x(t,1))
        x(t+1, 2) = beta*(x(t, 5) - x(t,2));
        x(t+1, 3) = beta*(x(t,6) - x(t,3));
        x(t+1, 4)=alpha_0 + (alpha/(1+(x(t,3))^n)) - (x(t,4));
        x(t+1, 5)=alpha_0 + (alpha/(1+(x(t,1)^n))) - (x(t,5));
        x(t+1, 6)=alpha_0 + (alpha/(1+(x(t,2))^n))-(x(t,6));
    end
    figure()
    plot (time,x(:,1),'--r',time,x(:,2),':b',time,x(:,3),'-.k')
    
else
    [t, output] = ode45(@repressilator, t, init_cond);
    figure()
    plot(t, output(:,1), t, output(:,2), t, output(:,3))
    
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

I feel like I'm doing something really silly. Any ideas? Thanks in advance.

Answer 1

[Riccardo Scorretti]

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.

@Torsten

[Stacy Genovese]

I just woke up and said to myself "You need to multiply by dt!" LOL. Thank you so much.

Answer 2

[Torsten]

alpha_0 = 0.03;
alpha = 298.2;
beta = 0.2;
n = 2;
time = [0:0.1:300];
dt = 0.1;
method = 1;
init_cond = [10 10 0 0 0 0];

if method == 1
    x = zeros(numel(time),6);
    x(1,1) = 10;
    x(1,2) = 10;
    for t = 1:3000
        x(t+1,1) = x(t,1) + dt*beta*(x(t,4) - x(t,1));
        x(t+1,2) = x(t,2) + dt*beta*(x(t,5) - x(t,2));
        x(t+1,3) = x(t,3) + dt*beta*(x(t,6) - x(t,3));
        x(t+1,4)=  x(t,4) + dt*(alpha_0 + (alpha/(1+(x(t,3))^n)) - (x(t,4)));
        x(t+1,5)=  x(t,5) + dt*(alpha_0 + (alpha/(1+(x(t,1))^n)) - (x(t,5)));
        x(t+1,6)=  x(t,6) + dt*(alpha_0 + (alpha/(1+(x(t,2))^n)) - (x(t,6)));
    end
    figure(1)
    plot (time,x(:,1),'--r',time,x(:,2),':b',time,x(:,3),'-.k')
    
else
    [t, output] = ode45(@(t,y)repressilator(t,y,alpha_0,alpha,beta,n), time, init_cond);
    figure(2)
    plot(t, output(:,1), t, output(:,2), t, output(:,3))
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d = repressilator(t, ic,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