Why ODE45 produces different answers for the same ODE function?

2 Ansichten (letzte 30 Tage)
JING JIAO
JING JIAO am 19 Mai 2019
Bearbeitet: JING JIAO am 22 Mai 2019
Here I used ode45 to solve one simple ODE function with three variables x(1), x(2) and x(3). I saved the final output 'pop' as 'Pop_tracking.csv' that includes four columns: certain time t, x(1), x(2) and x(3). The output inside the function 'tracking_dpop2-x2-5-22-19.csv' includes three columns: some time t, dPop(2) and x(2). I assumed that if the time in the above csv files match, the x(2) values in both files should be equal to each other, right? However, they are quiet different: in 'tracking_dpop2-x2-5-22-19.csv', once dPop(2) is negative and x(2)<5, x(2) is automatically set to 0, while in 'Pop_tracking.csv', the similar time point, x(2) is equal to some values < 5. I This should not happen based on the 'if' loop in SI_eq function.
I wonder why this happens? Is there anything wrong about my coding? Thanks for your attention.
Here is my coding:
function [t] = Test_of_ODE45(r,K,mu,beta,gamma,alpha)
r=0.05;
K=100;
mu=0.05;
alpha=0.01;
beta=0.03;
gamma=0.05;
x0=[5 1 0];
tspan = 0:1:100;
options =odeset('RelTol', 1e-8);
[t, pop]=ode45(@SI_eq, tspan, x0,options, r,K,mu,beta,gamma,alpha);
dlmwrite('Pop_tracking_5-22-19.csv',[t pop],'-append');
S=pop(:,1);
I=pop(:,2);
R=pop(:,3);
% plot everything
subplot(4,1,1)
h = plot(t,S);
xlabel 'Time';
ylabel 'S'
subplot(4,1,2)
h=plot(t,I);
xlabel 'Time';
ylabel 'I'
subplot(4,1,3)
h=plot(t,R);
xlabel 'Time';
ylabel 'R'
subplot(4,1,4)
h=plot(t,S+I+R);
xlabel 'Time';
ylabel 'S+I+R'
end
function dPop = SI_eq(t, pop, r, K, mu, beta, gamma,alpha)
x = pop(1:3);
dPop = zeros(3,1);
N=x(1)+x(2)+x(3);
dPop(1)=r*N*(1-N/K)-mu*x(1)-beta*x(1)*x(2);
dPop(2)= beta*x(1)*x(2)-mu*x(2)-gamma*x(2)-alpha*x(2);
dPop(3)= gamma*x(2)-mu*x(3);
if dPop(2)<0 & x(2)<0.5
x(2)=0;
end
dlmwrite('tracking_dpop2-x2-5-22-19.csv',[t dPop(2) x(2)],'-append');
end

Antworten (0)

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by