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.
采纳的回答
Riccardo Scorretti
2022-4-11
编辑：Riccardo Scorretti
2022-4-11
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 个）
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
类别
在 帮助中心 和 File Exchange 中查找有关 General Applications 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
