modeling - Bifurcation diagram of discrete SIR model in MATLAB

Question

I have a problem with my MATLAB code to display a figure of a Bifurcation diagram of discrete SIR model.

My model is:

S(n+1) = S(n) - h*(0.01+beta*S(n)*I(n)+d*S(n)-gamma*R(n))

I(n+1) = I(n) + h*beta*S(n)*I(n)-h*(d+r)*I(n)

R(n+1) = R(n) + h*(r*I(n)-gamma*R(n));

I tried out the code below but it keeps MATLAB busy for almost 30 mins and showed up no figure.

MATLAB code:

close all;
clear all;
clc;

%Model parameters

beta = 1/300;
gamma = 1/100;
D = 30;                  % Simulate for D days
N_t = floor(D*24/0.1);   % Corresponding no of hours
d = 0.001;
r = 0.07;

%Time parameters

dt = 0.01;
N = 10000;

%Set-up figure and axes

figure;

ax(1) = subplot(2,1,1);
hold on
xlabel ('h');
ylabel ('S');

ax(2) = subplot(2,1,2);
hold on
xlabel ('h');
ylabel ('I');

%Main loop

for h = 2:0.01:3

    S = zeros(N,1);
    I = zeros(N,1);
    R = zeros(N,1);

    S(1) = 8;
    I(1) = 5;
    R(1) = 0;

    for n = 1:N_t
        S(n+1) = S(n) - h*(0.01+beta*S(n)*I(n)+d*S(n)-gamma*R(n)); 
        I(n+1) = I(n) + h*beta*S(n)*I(n)-h*(d+r)*I(n);
        R(n+1) = R(n) + h*(r*I(n)-gamma*R(n));
    end

    plot(ax(1),h,S,'color','blue','marker','.');
    plot(ax(2),h,I,'color','blue','marker','.');

end

Any suggestions?

question from:https://stackoverflow.com/questions/65540702/bifurcation-diagram-of-discrete-sir-model-in-matlab

Answer 1

It was very slow because you were plotting a single value of h versus a vector S which had 7200 points. I assumed that you only want to plot the last value of S versus h. So replacing S with S(end) in the plot command changes everything. You really didn't need to use hold and it's better call plot once for each axis, so here is how I would do it:

beta = 1/300;
gamma = 1/100;
D = 30;                 % Simulate for D days
N_t = floor(D*24/0.1);   % Corresponding no of hours
d = 0.001;
r = 0.07;

%%Time parameters
dt = 0.01;
N = 10000;

%%Main loop
h = 2:0.01:3;
S_end = zeros(size(h));
I_end = zeros(size(h));
for idx = 1:length(h)
    S = zeros(N_t,1);
    I = zeros(N_t,1);
    R = zeros(N_t,1);
    S(1) = 8;
    I(1) = 5;
    R(1) = 0;

    for n=1:(N_t - 1)
        S(n+1) = S(n) - h(idx)*(0.01+beta*S(n)*I(n)+d*S(n)-gamma*R(n)); 
        I(n+1) = I(n) + h(idx)*beta*S(n)*I(n) - h(idx)*(d+r)*I(n);
        R(n+1) = R(n) + h(idx)*(r*I(n)-gamma*R(n));
    end
    S_end(idx) = S(end);
    I_end(idx) = I(end);
end

figure(1)
subplot(2,1,1);
plot(h,S_end,'color','blue','marker','.');
xlabel ('h');
ylabel ('S');

subplot(2,1,2);
plot(h,I_end,'color','blue','marker','.');xlabel ('h');
xlabel ('h');
ylabel ('I');

This now runs in 0.2 seconds on my computer.