How do you plot the solution of a system of odes against a parameter of the function and not time?

4 Ansichten (letzte 30 Tage)
Hello,
I am trying to recreate some figures from a mathematical modelling paper, however I have arrived at one figure and unsure on how to plot it. This is the code that I have done for the function, it is an SIS model made up of two susceptible compartments.
function f=f(t,y)
f=zeros(3,1);
m=0.8;
A=200;
theta=0.004;
%lambda=0.00002; %R0>1
%lambda=0.00000989; %R0<1
n=0.7;
xi=0.012;
mu=0.035;
delta=0.3;
lambda1=0.1;
phi=0.01;
f(1)=m*A-theta*y(1)-lambda*y(1).*y(3)+n*xi*y(3)-mu*y(1);
f(2)=(1-m)*A+theta*y(1)-lambda*(1+delta*lambda1)*y(2).*y(3)+(1-n)*xi*y(3)-mu*y(2);
f(3)=lambda*y(1).*y(3)+lambda*(1+delta*lambda1)*y(2).*y(3)-(xi+phi+mu)*y(3);
This is the figure that I am wanting to replicate:
My main obstacle is trying to plot delta against the solution of y(:3), since it is not time. I have considered using a method, such as implict euler, however I am unsure on how to specify which is the argument of a function when the Matlab function has more inputs. I have tried using ode45 with this code
[t,u]=ode45(@f,[0 3000], [500 300 10]);
plot(t,u(:,3),'LineWidth',2)
but seem to have no luck in trying to change t to delta. Any help would be greatly appreciated.

Akzeptierte Antwort

Cris LaPierre
Cris LaPierre am 5 Aug. 2021
Bearbeitet: Cris LaPierre am 5 Aug. 2021
In your code, delta and m are constants. It is likely they ran multiple simulations to create that figure, varying the contants and capturing the desired value of y.
m = 0.2:0.2:0.8;
delta = 0:.05:2;
for a = 1:length(m)
for b = 1:length(delta)
[t,u]=ode45(@f,[0 3000], [500 300 10],[],m(a),delta(b));
val(b,a) = u(end,3);
end
end
plot(delta,val)
legend("m="+m','Location','Northwest')
function f=f(t,y,m,delta)
f=zeros(3,1);
% m=0.8;
A=200;
theta=0.004;
lambda=0.00002; %R0>1
%lambda=0.00000989; %R0<1
n=0.7;
xi=0.012;
mu=0.035;
% delta=0.3;
lambda1=0.1;
phi=0.01;
f(1)=m*A-theta*y(1)-lambda*y(1).*y(3)+n*xi*y(3)-mu*y(1);
f(2)=(1-m)*A+theta*y(1)-lambda*(1+delta*lambda1)*y(2).*y(3)+(1-n)*xi*y(3)-mu*y(2);
f(3)=lambda*y(1).*y(3)+lambda*(1+delta*lambda1)*y(2).*y(3)-(xi+phi+mu)*y(3);
end

Weitere Antworten (0)

Kategorien

Mehr zu Loops and Conditional Statements finden Sie in Help Center und File Exchange

Tags

Produkte


Version

R2018b

Community Treasure Hunt

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

Start Hunting!

Translated by