matlab传染病模型分段

Question

Matlab中传染病模型，对在不同时间区间中传染概率赋值，但程序一直报错，不知道什么原因，想麻烦帮看一下。

function dY= Fsisd(t,Y,b,g,m)
  
S=Y(1);%S为易感状态
I=Y(2);%I为感染状态
D=Y(3);%D为死亡状态

dY = zeros(3,1);
%b为传染概率，m为恢复率，g为死亡率
dY(1) = -b*S*I + m*I;
dY(2) =  b*S*I - g*I-m*I;
dY(3) =  g*I

end


clc;
clear;
close all;
warning off;
addpath(genpath(pwd));
%-- initial conditions
t=0:0.1:120;
b=0.21*(t<3)+0.16*(t>=3);%在时间t<3时，传染概率b为为0.21；在时间t>=3时，传染概率b为为0.16；
g=0.24*(t>=0); %死亡率的取值 
m=0.06*(t>=0); % 恢复率的取值
i0=0.2;    %  I状态初始占比
s0=0.8; %  S状态初始占比
d0=0;    %  R状态初始占比

T=120;   % evaluation time

 S0I0D0=[s0 i0 d0];    % initial condictions Vector
  Tspam=[0:0.1:T]; % time interval
  %-- Numerical Integration
  
  [T,Y] = ode45(@(t,Y) Fsisd(t,Y,b,g,m),Tspam,S0I0D0);

S=Y(:,1); % Solution S
I=Y(:,2); % Solution I
D=1-S-I; % Solution D 


%-----   plots -----
plot(T,S,'k');
hold on;
grid on;
plot(T,I,'r--');
plot(T,D,'b-.');
axis([0 120,0,1])
title('典当行-SISD model ')
xlabel('Time')
ylabel('各状态比例')
legend('S','I','D','Location','best')

Answer 1

该回答引用ChatGPT GPT-4

运行结果

代码如下

clc;
clear;
close all;
warning off;
addpath(genpath(pwd));
%-- initial conditions
t=0:0.1:120;
b_func=@(t) (0.21*(t<3)+0.16*(t>=3)); %在时间t<3时，传染概率b为为0.21；在时间t>=3时，传染概率b为为0.16；
g_func=@(t) (0.24*(t>=0)); %死亡率的取值 
m_func=@(t) (0.06*(t>=0)); % 恢复率的取值
i0=0.2;    %  I状态初始占比
s0=0.8; %  S状态初始占比
d0=0;    %  R状态初始占比
 
T=120;   % evaluation time
 
 S0I0D0=[s0 i0 d0];    % initial condictions Vector
  Tspam=[0:0.1:T]; % time interval
  %-- Numerical Integration
  
  [T,Y] = ode45(@(t,Y) Fsisd(t,Y,b_func,g_func,m_func),Tspam,S0I0D0);
 
S=Y(:,1); % Solution S
I=Y(:,2); % Solution I
D=1-S-I; % Solution D 

%-----   plots -----
plot(T,S,'k');
hold on;
grid on;
plot(T,I,'r--');
plot(T,D,'b-.');
axis([0 120,0,1])
title('典当行-SISD model ')
xlabel('Time')
ylabel('各状态比例')
legend('S','I','D','Location','best')

function dY= Fsisd(t,Y,b_func,g_func,m_func)
  
S=Y(1);%S为易感状态
I=Y(2);%I为感染状态
D=Y(3);%D为死亡状态
 
dY = zeros(3,1);
%b为传染概率，m为恢复率，g为死亡率
b=b_func(t);
g=g_func(t);
m=m_func(t);

dY(1) = -b*S*I + m*I;
dY(2) =  b*S*I - g*I-m*I;
dY(3) =  g*I;
 
end