ode45 solver SIR model

조회 수: 63 (최근 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 CenterFile 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