ODE solver vs discrete integration

3 views (last 30 days)
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.

Accepted Answer

Riccardo Scorretti
Riccardo Scorretti on 11 Apr 2022
Edited: Riccardo Scorretti on 11 Apr 2022
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.
  1 Comment
Stacy Genovese
Stacy Genovese on 12 Apr 2022
I just woke up and said to myself "You need to multiply by dt!" LOL. Thank you so much.

Sign in to comment.

More Answers (1)

Torsten
Torsten on 11 Apr 2022
Edited: Torsten on 11 Apr 2022
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

Tags

Products


Release

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!