MATLAB中离散SIR模型的分岔图

编程问答 2022-06-01

问题描述

我的 MATLAB 代码出现问题，无法显示离散 SIR 模型的分岔图。

我的模型是：

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));

我尝试了下面的代码，但它让 MATLAB 忙碌了将近 30 分钟，并且没有显示任何数字。

MATLAB 代码：

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,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),I,'.');

end

有什么建议吗？

解决方法

它非常慢，因为您绘制的是 h 的单个值与具有 7200 个点的向量 S。我假设您只想绘制 S 与 h 的最后一个值。因此，在 S 命令中将 S(end) 替换为 plot 会改变一切。你真的不需要使用 hold 并且最好为每个轴调用一次 plot，所以我会这样做：

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,2);
plot(h,I_end,'.');xlabel ('h');
xlabel ('h');
ylabel ('I');

现在在我的电脑上运行时间为 0.2 秒。

matlabmodeling