手搓传染病模型（SEIA-拓展）

继上文手搓SEIA基础模型后，想要加两个创新点

1.传播过程中加入干预措施

假设传染率（β）随时间以某种函数形式变化，例如指数衰减：

另一种可能的调整是分段常数形式：

2.暴露者转化为有症状感染者的概率（ p ）可能随时间变化

可以将p定义为一个动态函数：

这里有一个难点，就是如何去刻画分段函数？

我们话不多说，开始手搓

% 模型参数
N = 21858000; % 总人数
I0 = 1000; % 初始有症状感染人数
A0 = 500; % 初始无症状感染人数
E0 = 0; % 初始暴露人数
S0 = N - I0 - A0 - E0; % 初始易感人数
beta0 = 0.5; % 初始传染率
w = 0.3; % 暴露转化速率
p0 = 0.6; % 初始发病概率
num_days = 160; % 模拟天数
t_intervention = 30; % 干预措施开始时间（天）
reduction_ratio = 0.5; % 传播率降低的比例
lambda = 0.05; % 传播率下降的速度
adjustable_ratio = 0.8; % 发病概率的调整比例% 定义时间依赖的传播率beta(t)
beta = @(t) (t <= t_intervention) * beta0 + (t > t_intervention) * beta0 * reduction_ratio * exp(-lambda * (t - t_intervention));
% 定义时间依赖的发病概率p(t)
p = @(t) (t <= t_intervention) * p0 + (t > t_intervention) * p0 * adjustable_ratio;
% x(1):初始有症状感染人群I，
% x(2):易感人群S， 
% x(3):无症状感染人群A，
% x(4):初始暴露人群E
dxdt = @(t, x) [p(t) * w * x(4); % dIdt-beta(t) * x(2) * (x(1) + x(3)) / N; % dSdt(1 - p(t)) * w * x(4); % dAdtbeta(t) * x(2) * (x(1) + x(3)) / N - w * x(4); % dEdt];[t, y] = ode45(dxdt, 1: num_days, [I0, S0, A0, E0]);
hold on
plot(t, y(:, 1));
plot(t, y(:, 2));
plot(t, y(:, 3));
plot(t, y(:, 4));
xline(t_intervention, '--b', sprintf('干预措施开始 (%d 天)', t_intervention));
legend('感染人数I', '易感人数S', '无症状感染人群A', '暴露人群E');

看下效果

Over！