ode45 solver SIR model

40 次查看（过去 30 天）
Ollie Stringer
Ollie Stringer 2021-1-5
编辑： Mischa Kim 2021-1-6
I've used ode45 to solve a simple SIR model, I've got the graph to work as I wish but I'm strugling to output any numerical values to discuss.
I'm wanting to find where dI/dt = 0, for the time where the pandemic will be at its peak and the area under the curve for the total number of infected.
Here's my (probably very messy!) code:
%user parameters
N = 45742000; %total population
I0 = 1; %initial infected population
r0 = 12.9 %reproductive value
R0 = 0;%initial recovered population
i_period = 9; %duration of infectious period
tspan = [1, 50]; %length of time to run simulation over
%rate constant calculations
mu = 1 / i_period %recovery rate
beta = r0 * mu / N %infection rate
S0 = N - I0 - R0 %initial susceptible population
N0 = I0 + R0 + S0; %total population
%---------------------------------------------------
%feeding parameters into function
pars = [beta, mu, N, r0];
y0 = [S0 I0 R0];
Running SIR model function
%using the ode45 function to perform intergration
[t,y] = ode45(@sir_rhs, tspan, y0, [], pars);
figure()
plot(t,y(:,2), 'r');
xlim(tspan);
ylabel('Population (n)');
xlabel('Time (days)');
legend('Infected','Location','SouthEast');
function f = sir_rhs(t,y,pars)
f = zeros(3,1);
f(1) = -pars(1)*y(1)*y(2);
f(2) = pars(1)*y(1)*y(2) - pars(2)*y(2);
f(3) = pars(2) * y(2);
end
---------------------------------------------------------------------------------------------------
function f = sir_rhs(t, y, pars); %function contains differential equations
f = zeroes(4, 1); %creates an empty matrix which will be filled with values for susceptible, infected, recovered and total population, respectively
beta = pars(1);
mu = pars(2);
N = pars(3);
r = pars(4);
S= y(1);
I= y(2);
R= y(3);
n = y(4);
f(1) = -beta*S*I; %susceptible population differential equation
f(2) = beta*S*I-mu*I; %infected populatin differential equation
f(3) = mu * I; %recovered population differential equation
end
Im only interested in solving f(2) where it's equal to zero, which I think will be where the disease is at maximum or where it ends (I don't think it reaches zero more than once here because there's a single peak. I could use the sum of y(2) at each time point to calculate the total infected throughout the time course but there must be a simpler alternative!

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

采纳的回答

Mischa Kim
Mischa Kim 2021-1-5
编辑：Mischa Kim 2021-1-6
Hi Ollie,
to get the area under the curve, just add another integration variable that integrates over the curve, y(2), I called it f(4). To locate the peak of the curve use an events function as shown below. ve contains the time stamp of the peak.
y0 = [S0 I0 R0 0];
%using the ode45 function to perform intergration
options = odeset('Events',@events);
[t,y,ve,~,~] = ode45(@sir_rhs, tspan, y0, options, pars);
figure()
plot(t,y(:,2),t,y(:,4));
xlim(tspan);
ylabel('Population (n)');
xlabel('Time (days)');
legend('Infected','Location','SouthEast');
function f = sir_rhs(~,y,pars)
f = zeros(4,1);
f(1) = -pars(1)*y(1)*y(2);
f(2) = pars(1)*y(1)*y(2) - pars(2)*y(2);
f(3) = pars(2) * y(2);
f(4) = y(2);
end
function [value,isterminal,direction] = events(~,y,pars)
value = pars(1)*y(1)*y(2) - pars(2)*y(2);
isterminal = 0; % Do not stop the integration when zero is detected
direction = 0; % Detect all zeros
end

更多回答（0 个）

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by