Plotting the derivative of infected population SI model

1 次查看（过去 30 天）
Amal Matrafi
Amal Matrafi 2023-12-16
评论： Sam Chak 2023-12-16
Hello;
I'm trying to draw the following model
I tried the code ode45 but it didn't work.
Is there a specific way to link to the same image? Thank you.
  3 个评论
显示 1更早的评论
Amal Matrafi
Amal Matrafi 2023-12-16
clc;
clear all
close all;
N=1000;
tend=300;
I0=10;
tspan = [0,tend];S0 = N - I0;y0 = [S0; I0];
opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
beta=0.1;
[t,y] = ode45(@(t,y) SIRfunc(t,y,beta,N), tspan, y0,opts);
plot(t,y(:,2),'b','LineWidth',1);
hold on
beta=0.2;
[t,y] = ode45(@(t,y) SIRfunc(t,y,beta,N), tspan, y0,opts);
plot(t,y(:,2),'r','LineWidth',1);
hold on
beta=0.25;
[t,y] = ode45(@(t,y) SIRfunc(t,y,beta,N), tspan, y0,opts);
plot(t,y(:,2),'LineWidth',1);
hold off
legend('\beta_{1}=0.1','\beta_{1}=0.2','\beta_{1}=0.25')
xlim([0 140])
ylim([0 1000])
function dydt = SIRfunc(~,y,beta,N)
dydt = [-beta/N*y(2)*y(1);
beta/N*y(2)*y(1)];
end

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Star Strider
Star Strider 2023-12-16
You need to plot the derivatives, not the solved equations.
One option is to use the gradient function:
dy2dt = gradient(y(:,2),t)
There are other, more direct (and probably more accurate) ways of calculating it from the original differential equation function. That requires a loop.
.
Sam Chak
Sam Chak 2023-12-16
Ran in:
Hi @Amal Matrafi
You can use deval() to obtain the first derivative. Alternatively, as suggested by @Star Strider, you can also use the gradient() approach to obtain the first derivative.
beta = [0.1, 0.2, 0.25];
N = 1000;
tend = 150;
I0 = 10;
tspan = [0,tend];
S0 = N - I0;
y0 = [S0; I0];
opts = odeset('RelTol', 1e-2, 'AbsTol', 1e-4);
for j = 1:numel(beta)
sol = ode45(@(t, y) SIRfunc(t, y, beta(j), N), tspan, y0, opts);
t = linspace(0, 150, 1501);
[y, yp] = deval(sol, t);
plot(t, yp(2,:)), hold on
end
grid on
hold off
xlabel('t'), ylabel('dI/dt')
legend('\beta_{1} = 0.1','\beta_{2} = 0.2', '\beta_{3} = 0.25')
%% SI Model
function dydt = SIRfunc(t, y, beta, N)
dydt = [-beta/N*y(2)*y(1);
beta/N*y(2)*y(1)];
end
  1 个评论
Sam Chak
Sam Chak 2023-12-16
Ran in:
Before learning to use deval(), I utilized the right-hand side of the state equation by directly substituting the solution from ode45(). This is pure math stuff!
beta = [0.1, 0.2, 0.25];
N = 1000;
tend = 150;
I0 = 10;
tspan = linspace(0, tend, 10*tend+1);
S0 = N - I0;
y0 = [S0; I0];
opts = odeset('RelTol', 1e-2, 'AbsTol', 1e-4);
for j = 1:numel(beta)
[t, y] = ode45(@(t, y) SIRfunc(t, y, beta(j), N), tspan, y0, opts);
dIdt = beta(j)/N*y(:,2).*y(:,1);
plot(t, dIdt), hold on
end
grid on
hold off
xlabel('t'), ylabel('dI/dt')
legend('\beta_{1} = 0.1','\beta_{2} = 0.2', '\beta_{3} = 0.25')
%% SI Model
function dydt = SIRfunc(t, y, beta, N)
dydt = [-beta/N*y(2)*y(1);
beta/N*y(2)*y(1)];
end

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by