SIS模型Python实现：β参数优化与二阶导数求解问询

1. 明确二阶导数为0的核心条件 #

先从SIS连续模型的基础方程推导：
总人数恒定 ( S + I = N )，感染人数的一阶导数为：

dI/dt = β*I*(N-I)/N - γ*I = I*(β*(1 - I/N) - γ)

对一阶导数求导得到二阶导数：

d²I/dt² = (dI/dt) * [β*(1 - 2I/N) - γ]

要让 ( d²I/dt²=0 )，有两种情况：

• 情况1：( dI/dt=0 )，即系统达到稳态，此时 ( I^* = N*(1 - γ/β) )，但这是长期稳态，大概率不是第10时间步的目标状态。
• 情况2：( β*(1 - 2I/N) - γ = 0 )，整理得 ( β = γ / (1 - 2I(t)/N) )。这里的 ( I(t) ) 是第10时间步的感染人数，而 ( I(t) ) 本身是β的函数（解析解含指数项），因此这是一个隐式方程，无法直接代数求解，必须用数值方法。

2. 数值求解β的两种可行方法 #

方法一：基于解析解的根查找 #

SIS模型的解析解（当 ( β≠γ ) 时）为：

I(t) = (I0 * N * (β - γ)) / (I0*(β - γ) + (β*N - γ*I0)*exp(-(β - γ)*t))

将 ( t=10 ) 代入情况2的等式，得到关于β的隐式方程：

β = γ / (1 - 2*[I0*N*(β-γ)/(I0*(β-γ)+(βN - γI0)*exp(-(β-γ)*10))]/N)

可以用牛顿-拉夫逊法求解该方程的根，步骤如下：

1. 定义目标函数 ( f(β) = β - γ / (1 - 2*I(β,10)/N) )，其中 ( I(β,10) ) 是代入β后的第10时间步感染人数。
2. 选择合理的初始β值（需满足 ( β>γ )，否则感染人数会单调衰减，无有效零点）。
3. 迭代更新β：( β_{new} = β_{old} - f(β_{old})/f’(β_{old}) )，直到 ( |f(β)| < 1e-6 )（精度可调整）。

方法二：基于离散迭代的数值优化 #

如果你的代码是用离散时间步（欧拉法、Runge-Kutta法）计算 ( I(t) )，可以用非线性最小二乘法或梯度下降法优化β：

1. 定义损失函数 ( L(β) = [d²I/dt²(10)]² )，目标是最小化该函数（当 ( L(β)→0 ) 时满足条件）。
2. 用中心差分近似离散二阶导数：

d²I/dt² ≈ (I(t+Δt) - 2*I(t) + I(t-Δt)) / (Δt²)

其中 ( t=10 ) 对应第10时间步，需计算第9、10、11步的 ( I ) 值。
3. 调用Python优化库（如scipy.optimize的root或minimize函数）求解：

• 用root时，直接将目标设为 ( d²I/dt²(10)=0 )；
• 用minimize时，最小化 ( [d²I/dt²(10)]² )。

3. Python实现示例（基于解析解+scipy优化） #

import numpy as np
from scipy.optimize import root

def calculate_I(β, t, N, γ, I0):
    """计算t时刻的感染人数（解析解）"""
    if β == γ:
        return I0 / (1 + (γ*I0/N)*t)
    numerator = I0 * N * (β - γ)
    denominator = I0*(β - γ) + (β*N - γ*I0)*np.exp(-(β - γ)*t)
    return numerator / denominator

def target_function(β, N, γ, I0, t_target):
    """目标函数：满足d²I/dt²=0的隐式方程"""
    I = calculate_I(β, t_target, N, γ, I0)
    return β - γ / (1 - 2*I/N)

# 给定参数
N = 1000  # 总人数
γ = 0.2   # 恢复率
I0 = 10   # 初始感染人数
t_target = 10  # 目标时间步

# 初始猜测β（需大于γ）
initial_β = 0.3

# 求解最优β
result = root(target_function, initial_β, args=(N, γ, I0, t_target))

if result.success:
    optimal_β = result.x[0]
    print(f"最优β值：{optimal_β:.4f}")
    # 验证二阶导数
    dt = 0.01
    I_prev = calculate_I(optimal_β, t_target - dt, N, γ, I0)
    I_curr = calculate_I(optimal_β, t_target, N, γ, I0)
    I_next = calculate_I(optimal_β, t_target + dt, N, γ, I0)
    d2I = (I_next - 2*I_curr + I_prev) / (dt**2)
    print(f"t=10时二阶导数：{d2I:.4f}")
else:
    print("求解失败：", result.message)

4. 关键注意事项 #

• 初始β必须大于γ，否则感染人数单调递减，二阶导数无有效零点（除非初始就是稳态）。
• 若出现 ( 2*I(t)/N ≥1 )，会导致β为负数，无实际意义，需调整初始β确保 ( I(t) < N/2 )。
• 离散迭代时，时间步长Δt不宜过大，否则会引入较大数值误差。

内容的提问来源于stack exchange，提问作者Syed Ali Mohsin Bukhari