我试图重现图5.6（附件）摘自教科书《人类和动物传染病建模（official code repo）”（Keeling 2008）验证了我实施的季节性强制SEIR是否（流行病学模型）是正确的。教科书中执行季节强迫的官方程序表明，β 1的值过大会导致数值错误，但是如果图中的Beta值为1，而这并不导致数值错误，那么原则上这不应该是问题的原因，我的实现可以正确地生成第0行、第0列和第1行的图，图5.6的第0列，但由于感染细胞分数的数值解，我的图中没有剩余细胞的输出（参见底部的代码）生成0（和ln（0）--〉-inf）。
我确实收到以下警告：
ODEint警告：在此呼叫上完成的额外工作
第68章：我的意思是：我的意思是：我的意思是：我的意思是：我的意思是：我的意思是：运行时警告：在记录infected = np.log（odeint（
第68章：我的意思是：我的意思是：我的意思是：我的意思是：我的意思是：我的意思是：运行时警告：在日志infected = np.log（odeint（
以下是教科书上的数字：

我的数字在同一时间范围内（990 - 1000年）。自然对数采取的部分感染：

我的数字，但时间范围较短（0 - 100年）。感染率的自然对数。感染人口的数值解似乎在5年和20年之间失败的大多数季节性参数（β 1和R 0）：

我的图与上面的时间范围较短，但没有感染分数的自然对数。

代码重现我的身材：

# Code to minimally reproduce figure

import itertools
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def seir(y, t, mu, sigma, gamma, omega, beta_zero, beta_one):
    """System of diff eqs for epidemiological model.

    SEIR stands for susceptibles, exposed, infectious, and
    recovered populations. 

    References:
        [SEIR Python Program from Textbook](http://homepages.warwick.ac.uk/~masfz/ModelingInfectiousDiseases/Chapter2/Program_2.6/Program_2_6.py)
        [Seasonally Forced SIR Program from Textbook](http://homepages.warwick.ac.uk/~masfz/ModelingInfectiousDiseases/Chapter5/Program_5.1/Program_5_1.py)
    """
    s, e, i = y
    beta = beta_zero * (1 + beta_one * np.cos(omega * t))
    sdot = mu - (beta*i + mu)*s
    edot = beta*s*i - (mu  + sigma)*e
    idot = sigma*e - (mu + gamma)*i
    return sdot, edot, idot

def solve_beta_zero(basic_reproductive_rate, gamma):
    """Defined in the last paragraph of pg. 159 of textbook Keeling 2008."""
    return gamma * basic_reproductive_rate

# Model parameters (see Figure 5.6 description)

mu = 0.02 / 365
sigma = 1/8
gamma = 1/5
omega = 2 * np.pi / 365  # frequency of oscillations per year

# Seasonal forcing parameters

r0s = [17, 10, 3]
b1s = [0.02, 0.1, 0.225]

# Permutes params to get tuples matching row i column j params in figure

# e.g., [(0.02, 17), (0.02, 10) ... ]

seasonal_params = [p for p in itertools.product(*(b1s, r0s))]

# Initial Conditions: I assume these are proportions of some total population

s0 = 6e-2
e0 = i0 = 1e-3
initial_conditions = [s0, e0, i0]

# Timesteps

nyears = 1000
days_per_year = 365
ndays = nyears * days_per_year
timesteps = np.arange(1, ndays+1, 1)

# Range to slice data to reproduce my figures

# NOTE: CHange the min slice or max slice for different ranges

min_slice = 990  # or 0
max_slice = 1000 # or 100
sliced = slice(min_slice * days_per_year, max_slice * days_per_year)
x_ticks = timesteps[sliced]/days_per_year

# Define figure

nrows = 3
ncols = 3
fig, ax = plt.subplots(nrows, ncols, sharex=True, figsize=(15, 8))

# Iterate through parameters and recreate figure

for i in range(nrows):
    for j in range(ncols):

        # Get seasonal parameters for this subplot
        beta_one = seasonal_params[i * nrows + j][0]
        basic_reproductive_rate = seasonal_params[i * nrows + j][1]

        # Compute beta zero given the reproductive rate
        beta_zero = solve_beta_zero(
            basic_reproductive_rate=basic_reproductive_rate, 
            gamma=gamma)

        # Numerically solve the model, extract only the infected solutions,
        # slice those solutions to the desired time range, and then natural 
        # log scale them
        solutions = odeint(
            seir, 
            initial_conditions, 
            timesteps, 
            args=(mu, sigma, gamma, omega, beta_zero, beta_one))

        infected_solutions = solutions[:, 2]

        log_infected = np.log(infected_solutions[sliced])

        # NOTE: To inspect results without natural log, uncomment the
        # below line
        # log_infected = infected_solutions[sliced]

        # DEBUG: For shape and parameter printing
        # print(
        #     infected_solutions.shape, 'R0=', basic_reproductive_rate, 'B1=', beta_one)

        # Plot results
        ax[i,j].plot(x_ticks, log_infected)

        # label subplot
        ax[i,j].title.set_text(rf'$(R_0=${basic_reproductive_rate}, $\beta_{1}=${beta_one})')

fig.supylabel('NaturalLog(Fraction Infected)')
fig.supxlabel('Time (years)')

免责声明：
我的短期解决方案是简单地将季节性参数的列表更改为能够产生该范围内数据的值，这充分说明了季节性强迫的影响，但关键是要重现该图，如果作者能够做到这一点，其他人也应该能够做到。