Why ODE45 produces different answers for the same ODE function?
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
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
0 Kommentare
Antworten (0)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!