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!
Best Answer