Why ODE45 produces different answers for the same ODE function?

1 Ansicht (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)

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by