Why ODE45 produces different answers for the same ODE function?

2 vues (au cours des 30 derniers jours)
JING JIAO
JING JIAO le 19 Mai 2019
Modifié(e) : JING JIAO le 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

Réponses (0)

Catégories

En savoir plus sur Programming dans Help Center et File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by