I am trying to solve this problem:
SIR Model for Spread of Disease
In the winter of 1968, New York City was swept by a new strain of influenza, the Hong Kong flu.
The SIR model represents the relationship between the population of the susceptible, infected and
recovered with ordinary differential equations given as:
dS
dt = f = − βSI
N (6)
dI
dt = g = βSI
N − γI (7)
dR
dt = h = γI, (8)
where S(t), I(t), R(t) are the susceptible, infected and recovered population at time t. β is the
transmission rate, γ represents the recovery rate, and N = S(t) + I(t) + R(t) is the fixed population.
The objective of this problem is to solve the system of differential equations and observe the change
in population over time, using the 4th order Runge-Kutta method. The constants for the 4th order
Runge Kutta method are defined as:
c1 = ∆t · f (tk, yk) (9)
c2 = ∆t · f
tk + 1
2 ∆t, yk + 1
2 c1
(10)
c3 = ∆t · f
tk + 1
2 ∆t, yk + 1
2 c2
(11)
c4 = ∆t · f (tk + ∆t, yk + c3) (12)
yk+1 = yk + 1
6 c1 + 1
3 c2 + 1
3 c3 + 1
6 c4 (13)
(a) Write a function that estimates the populations S(tk+1), I(tk+1) and R(tk+1). The inputs and
outputs of the function are given as follows:
function [Skp1, Ikp1, Rkp1] = RK4(f, g, h, dt, Sk, Ik, Rk)
% Computes the susceptible, infected and recovery population at tkp1
% using 4th order Runge-Kutta method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inputs:
% f: function handle of first ODE
% g: function handle of second ODE
% h: function handle of third ODE
% dt: time step
% Sk: Susceptible population at tk
% Ik: Infected population at tk
% Rk: Recovered population at tk
%Output:
% Skp1: Susceptible population at tkp1
% Ikp1: Infected population at tkp1
% Rkp1: Recovered population at tkp1
(Hint: Define the three anonymous functions as the following, in the main script)
f = @(s, i, r) -beta*s*i/N
g = @(s, i, r) beta*s*i/N -gamma*i
h = @(s, i, r) gamma*i
(b) Write a script that calculates the susceptible, infected and recovered populations at each timestep
from t = 0 to t = 200. Use the following simulation parameters and initial conditions: dt = 1,
β = 0.5, γ = 1/3, S(0) = 7, 900, 000, I(0) = 10 and R(0) = 0.
(c) Using results from (b), write code that plots the ratio of the susceptible, infected and recov-
ered populations as a function of time on a single figure. The figure should have a title, axis
labels, and legend. Save the plot using the saveas function and include the figure in the report.
Screenshots are not accepted. Discuss the findings of the plot, e.g., “Fig. 1 shows Y quantity
as a function of X quantity. We observe that [a few sentences on the findings].”
I wrote this code: (part b and c are written before a because they have to be before the function)
%% Part a
f = @(s, i, r) -beta*s*i/N;
g = @(s, i, r) beta*s*i/N -gamma*i;
h = @(s, i, r) gamma*i;
%% Part b
% Define the parameters
beta = 0.5;
gamma = 1/3;
N = 7.9e6; % Total population
dt = 1; % Time step
t_end = 200; % End time
% Define the initial conditions
S(1) = 7.9e6 - 10; % Susceptible population
I(1) = 10; % Infected population
R(1) = 0; % Recovered population
% Simulate the SIR model using RK4
for t = 1:t_end
[S(t+1), I(t+1), R(t+1)] = RK4(@(s, i, r) -beta*s*i/N, @(s, i, r) beta*s*i/N - gamma*i, @(s, i, r) gamma*i, dt, S(t), I(t), R(t));
end
%% Part c
% Calculate the ratios
Sr = S ./ N;
Ir = I ./ N;
Rr = R ./ N;
% Plot the results
plot(0:t_end, Sr, 'b', 0:t_end, Ir, 'r', 0:t_end, Rr, 'g', 'LineWidth', 2);
title('SIR Model for the Hong Kong Flu (1968) in NYC');
xlabel('Time (days)');
ylabel('Ratio of Population');
legend('Susceptible', 'Infected', 'Recovered');
saveas(gcf, 'SIR_model.png');
function [Skp1, Ikp1, Rkp1] = RK4(f, g, h, dt, Sk, Ik, Rk)
% Computes the susceptible, infected and recovery population at tkp1
% using 4th order Runge-Kutta method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inputs:
% f: function handle of first ODE
% g: function handle of second ODE
% h: function handle of third ODE
% dt: time step
% Sk: Susceptible population at tk
% Ik: Infected population at tk
% Rk: Recovered population at tk
%Output:
% Skp1: Susceptible population at tkp1
% Ikp1: Infected population at tkp1
% Rkp1: Recovered population at tkp1
% Define the constants
c1 = dt * f(Sk, Ik, Rk);
c2 = dt * g(Sk + 0.5*c1, Ik + 0.5*c1, Rk);
c3 = dt * g(Sk + 0.5*c2, Ik + 0.5*c2, Rk);
c4 = dt * h(Sk + c3, Ik + c3, Rk);
% Calculate the new values
Skp1 = Sk + 1/6 * (c1 + 2*c2 + 2*c3 + c4);
Ikp1 = Ik + 1/6 * (c1 + 2*c2 + 2*c3 + c4);
Rkp1 = Rk + 1/6 * (c1 + 2*c2 + 2*c3 + c4);
end
As you can see, i defined f at the beginning, but i keep getting the error: Unrecognized fiunction or variable f no matter what i do, i really need help. thank you.
