MATLAB 代码主要用于研究传染病传播模型中接触追踪措施对疾病传播的影响-CSDN博客

clear
close all
step = 0.05; % 定义数值求解的步长，用于离散化时间变量。
K = 4; 
beta=0.568; 
gamma=0.2; 
g=0.3; 
R0 = beta*K/(beta+gamma);%基本再生数
bmax =8; % 最大传染期

%theta_a函数
theta_a = @(x) (exp(-(beta+gamma)* x));
epsilon_d = 1;
%症状比例 epsilon_s：表示出现症状的个体比例
epsilon_s = 0.85; % percentage of symptomatic individuals

%诊断参数：定义最大诊断时间 dmax、诊断延迟 delay_diagnosis 和被诊断的有症状个体比例 epsilon_d，并定义诊断概率密度函数 density_diagnosis 和诊断存活函数 surv_diagnosis_f
dmax = 15;
% 已知对数正态分布的均值和标准差
mean_incubation = 3.3091;
std_incubation = 2.3557;
% 计算 mu 和 sigma
mu = log(mean_incubation / sqrt(1 + (std_incubation^2 / mean_incubation^2)));
sigma = sqrt(log(1 + (std_incubation^2 / mean_incubation^2)));
% 定义对数正态分布的累计分布函数
cumulative_distribution_f = @(s) logncdf(s, mu, sigma);
surv_diagnosis_f = @(s) 1 - epsilon_s * cumulative_distribution_f(s);


%迭代以探索参数的不同组合，epsilon_c_vector 和 cmax_vector：分别定义接触追踪覆盖率和接触追踪窗口的取值范围
legendname = {};
epsilon_c_vector = 0:0.1:1;     
cmax_vector = 2:5;

%创建一个全零矩阵，矩阵的行数等于 epsilon_c_vector 的长度，列数等于 cmax_vector 的长度。
%后续将这个全零矩阵赋值给 Rct_matrix、rct_matrix、Theta_d_matrix 和 Thera_ct_matrix，这几个矩阵可能用于存储不同条件下（由 epsilon_c_vector 和 cmax_vector 组合确定）的计算结果。
Rct_matrix = zeros(length(epsilon_c_vector),length(cmax_vector));%zeros(length(epsilon_c_vector),length(cmax_vector))；
rct_matrix = Rct_matrix;
Theta_d_matrix = Rct_matrix;
Thera_ct_matrix = Rct_matrix;
%会生成一个包含 length(cmax_vector) 行的矩阵，每一行代表一种 RGB 颜色，后续可能会在绘图时使用这些颜色来区分不同的曲线
figure(5); hold on
plot([0 1],[1 1],'Color',[0.5 0.5 0.5],'LineStyle','--','LineWidth',1,'HandleVisibility','off'); % 参考线R=1
colorscode = lines(length(cmax_vector));%colorscode = lines(length(cmax_vector));

for index_x = 1:length(cmax_vector)
for index_y = 1:length(epsilon_c_vector)
%在每次循环中，从 cmax_vector 中取出第 index_x 个元素赋值给 cmax，从 epsilon_c_vector 中取出第 index_y 个元素赋值给 epsilon_c。这两个参数可能是模型中的关键参数，用于后续的计算。
    cmax = cmax_vector(index_x);
    epsilon_c = epsilon_c_vector(index_y);
 %线性系统的解
    nd = dmax/step;
    nc = cmax/step;
    nb = bmax/step;
    N = max([nb,nd,nc,nc+nb]);
    dgrid = step*(1:nd);
    bgrid = step*(1:nb);       %因此，bgrid 表示的是另一个等间距的值序列，起始于 step，终止于 step*nb，间距同样为 step。
    Ngrid = step*(1:N);
    %已知参数的初始化（函数离散化）
    
    theta = zeros(N,1);%用于存储每个时间步的传染性值，初始化为全零向量。
     for itau = 1:N
        tau = itau*step;
        theta(itau) = theta_a(tau);%通过循环遍历每个时间步
     end
    
    h_d = zeros(N,1);%诊断危险率，初始化为全零向量。
    surv_d = (1-epsilon_d*epsilon_s)*ones(N,1); % survival diagnosis%诊断存活函数，初始化为 (1 - epsilon_d * epsilon_s) 的全 1 向量。
    surv_d(1) = surv_diagnosis_f(step);
    h_d(1) = -log(surv_d(1))/step;%首先计算第一个时间步的诊断概率密度、存活函数和危险率。
     for itau = 2:(dmax/step)
        surv_d(itau) = surv_diagnosis_f(itau*step);
        h_d(itau) = - (log(surv_d(itau))-log(surv_d(itau-1)))/step;
    end
    surv_d(nd+1:end)=surv_d(nd);
    %通过公式计算再生数
    %trapz 是 MATLAB 中的梯形积分函数，用于对离散数据进行数值积分。它基于梯形法则，通过将相邻数据点之间的区域近似为梯形，然后计算这些梯形的面积之和来得到积分值。
% 假设前期已完成参数定义（K, beta, gamma, step, N等）
% 以及 surv_d、dens_d、h_d 等变量的初始化和循环计算
  Rd =beta* K* step*trapz(theta .* surv_d);
      r0 =beta*(K-1)-gamma;
      % 定义匿名函数表示目标函数
objectiveFunction = @(x) 1 - integral(@(t) K * beta * exp(-(beta + gamma + x) * t), 0, Inf);
r00 = fsolve(objectiveFunction, 1.5); % 初始猜测值为1.5
     rd = fsolve(@(x) 1- K*beta*step*trapz(theta.*surv_d.*exp(-x*step*(1:N)')), r00);%考虑诊断因素后的增长率 。
   %初始化接触者追踪概率
    x0 = zeros(2*N+1,1); %最后一项表示指数增长率
    x0(1:nc)=ones(1,nc);%将 x0 的前 nc 个元素赋值为 1。nc 是接触追踪窗口的离散化步数，这里这样初始化可能是假设在接触追踪窗口内，接触追踪相关的概率初始值为 1（可能是一种初始的默认状态）。
    x0(N+1:N+nc)=ones(1,nc);%将 x0 的前 nc 个元素赋值为 1。nc 是接触追踪窗口的离散化步数，这里这样初始化可能是假设在接触追踪窗口内，接触追踪相关的概率初始值为 1（可能是一种初始的默认状态）。
    x0(end)=rd;%将 x0 的最后一个元素赋值为 rd，即考虑诊断因素后的增长率，以此作为考虑接触追踪后增长率的初始猜测值。
    % linear_contact_tracing(x0(1:N),x0(N+1),step,nc,nd,epsilon_c,beta_mat,h_d,surv_d)
    options = optimoptions('fsolve','Display','none','MaxIter',100000);
  %使用 optimoptions 函数为 fsolve 求解器设置选项。'Display','none' 表示不显示求解过程的输出信息，'MaxIter',100000 表示最大迭代次数设置为 100000，以确保有足够的迭代次数来找到合适的解。
    Sol = fsolve(@(x) [x(1:N);x(N+1:2*N);1] - linear_contact_tracing(x(1:N),x(N+1:2*N),x(2*N+1),step,K,beta,gamma,g,nc,nd,epsilon_c,h_d,surv_d), x0);%使用 fsolve 求解上述非线性方程组，初始猜测值为 x0，求解结果存储在 Sol 中。
    h_1ct = Sol(1:N);
    h_2ct = Sol(N+1:2*N); 
    rct = Sol(2*N+1);%求解非线性方程组得到满足条件的解向量 Sol，然后从该向量中提取最后一个元素作为考虑接触追踪后的增长率。

    
    % 计算考虑接触追踪措施后的再生数 
    Rct = 0;
    surv_ct = zeros(N,1);
    for itau = 1:N
        surv_ct(itau) = exp(-step*K*sum(h_2ct(1:itau)))*exp(step*K*sum(exp(step*K*sum(-beta*(1:itau)))*h_2ct(1:itau)))*exp(-step*sum(h_1ct(1:itau)));
        Rct = Rct + step*K*beta*exp(-beta*itau)*surv_ct(itau)*surv_d(itau)*exp(-g*(itau));
    end
        
    Rct_matrix(index_y,index_x) = Rct;
    rct_matrix(index_y,index_x) = rct;

 
    % 防止传播
    Theta_d(index_y,index_x) = Rd/R0;
    Theta_ct(index_y,index_x) = Rct/R0;
end

 % 这部分代码的主要功能是创建一个新的图形窗口（编号为 5），并在该窗口中绘制有效再生数（Rct_matrix）随接触追踪覆盖率（epsilon_c_vector）的变化曲线。
    figure(5); hold on
    plot(epsilon_c_vector,Rct_matrix(:,index_x),'LineWidth',2,'Color',colorscode(index_x,:))
    legendname{index_x} = [num2str(cmax),'d tracing'];%将图例名称存储在 legendname 元胞数组中。num2str(cmax) 将变量 cmax 转换为字符串，然后与 'd tracing' 拼接成一个完整的图例名称。
    end