阻滞增长模型（Logistic模型）

年份	人口（百万）	年份	人口（百万）
1790	3.93	1900	76.21
1800	5.31	1910	92.23
1810	7.24	1920	106.02
1820	9.64	1930	123.20
1830	12.87	1940	132.16
1840	17.07	1950	151.33
1850	23.19	1960	179.32
1860	31.44	1970	203.30
1870	38.56	1980	226.54
1880	50.19	1990	248.71
1890	62.98	2000	281.42

计算过程

取 \( t = 0 \) 对应1790年，则时间变量 \( t \) 以10年为单位。

步骤一：利用非线性最小二乘法，以 \( K = 400 \)、\( r = 0.3 \)、\( A = 100 \) 为初始猜测，迭代求解得到：

\( K \approx 395.26 \)（百万人）
\( r \approx 0.3107 \)（每10年）
\( A \approx 99.56 \)

步骤二：验证拟合效果。拐点时间为：

\[ t^* = \frac{\ln 99.56}{0.3107} \approx 14.81 \text{（个10年）} \]

对应年份约为 \( 1790 + 148 = 1938 \) 年，即美国人口增速在1930-1940年代达到峰值，这与历史事实基本吻合。

步骤三：预测2010年人口：

\[ N(22) = \frac{395.26}{1 + 99.56 \times e^{-0.3107 \times 22}} \approx 311.7 \text{（百万人）} \]

实际2010年美国人口约308.7百万人，误差约1%。

Python代码实现

曲线拟合

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# 美国人口数据
years = np.array([1790, 1800, 1810, 1820, 1830, 1840, 1850, 1860, 1870,
                  1880, 1890, 1900, 1910, 1920, 1930, 1940, 1950, 1960,
                  1970, 1980, 1990, 2000])
population = np.array([3.93, 5.31, 7.24, 9.64, 12.87, 17.07, 23.19, 31.44,
                       38.56, 50.19, 62.98, 76.21, 92.23, 106.02, 123.20,
                       132.16, 151.33, 179.32, 203.30, 226.54, 248.71, 281.42])

# 时间变量（以10年为单位，1790年为起点）
t_data = (years - 1790) / 10


# 定义Logistic函数
def logistic(t, K, r, A):
    """
    Logistic增长模型
    参数:
        t: 时间
        K: 环境容纳量
        r: 内禀增长率
        A: 与初始条件相关的参数
    """
    return K / (1 + A * np.exp(-r * t))


# 非线性最小二乘拟合
p0 = [400, 0.3, 100]  # 初始猜测值
popt, pcov = curve_fit(logistic, t_data, population, p0=p0, maxfev=10000)

K_fit, r_fit, A_fit = popt
print(f"拟合结果:")
print(f"  环境容纳量 K = {K_fit:.2f} 百万人")
print(f"  内禀增长率 r = {r_fit:.4f} (每10年)")
print(f"  参数 A = {A_fit:.2f}")

# 计算拐点
t_inflection = np.log(A_fit) / r_fit
year_inflection = 1790 + t_inflection * 10
print(f"  拐点时间: {year_inflection:.1f} 年")
print(f"  拐点处人口: {K_fit/2:.2f} 百万人")

# 计算拟合优度 R^2
N_pred = logistic(t_data, *popt)
SS_res = np.sum((population - N_pred) ** 2)
SS_tot = np.sum((population - np.mean(population)) ** 2)
R_squared = 1 - SS_res / SS_tot
print(f"  拟合优度 R² = {R_squared:.6f}")

# 绘图
t_plot = np.linspace(0, 25, 500)
N_plot = logistic(t_plot, *popt)

plt.figure(figsize=(10, 6))
plt.scatter(1790 + t_data * 10, population, color='red', s=60,
            zorder=5, label='实际数据')
plt.plot(1790 + t_plot * 10, N_plot, 'b-', linewidth=2,
         label=f'Logistic拟合 (K={K_fit:.1f}, r={r_fit:.3f})')
plt.axhline(y=K_fit, color='gray', linestyle='--', alpha=0.7,
            label=f'容纳量 K={K_fit:.1f}')
plt.axhline(y=K_fit/2, color='green', linestyle=':', alpha=0.7,
            label=f'拐点 K/2={K_fit/2:.1f}')
plt.axvline(x=year_inflection, color='green', linestyle=':', alpha=0.7)

plt.xlabel('年份', fontsize=12)
plt.ylabel('人口（百万人）', fontsize=12)
plt.title('美国人口Logistic增长模型拟合', fontsize=14)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('logistic_fit.png', dpi=150)
plt.show()

数值求解（使用scipy.integrate）

from scipy.integrate import odeint, solve_ivp

# 定义Logistic微分方程
def logistic_ode(N, t, r, K):
    """Logistic微分方程右端函数"""
    return r * N * (1 - N / K)


# 参数设置
r = 0.3107
K = 395.26
N0 = 3.93  # 1790年人口

# 方法一：使用odeint求解
t_span = np.linspace(0, 25, 1000)
N_numerical = odeint(logistic_ode, N0, t_span, args=(r, K))

# 方法二：使用solve_ivp求解（更现代的接口）
def logistic_ode_ivp(t, N, r, K):
    """适用于solve_ivp的形式（注意t和N的顺序）"""
    return r * N * (1 - N / K)

sol = solve_ivp(
    logistic_ode_ivp,
    t_span=[0, 25],
    y0=[N0],
    args=(r, K),
    method='RK45',
    t_eval=np.linspace(0, 25, 1000),
    rtol=1e-8,
    atol=1e-10
)

# 与解析解对比
N_analytical = K / (1 + ((K - N0) / N0) * np.exp(-r * t_span))

# 计算数值误差
error = np.max(np.abs(N_numerical.flatten() - N_analytical))
print(f"数值解与解析解的最大误差: {error:.2e}")

# 绘制对比图
plt.figure(figsize=(10, 6))
plt.plot(1790 + t_span * 10, N_analytical, 'b-', linewidth=2,
         label='解析解')
plt.plot(1790 + sol.t * 10, sol.y[0], 'r--', linewidth=1.5,
         label='数值解 (RK45)')
plt.scatter(years, population, color='black', s=40,
            zorder=5, label='实际数据')
plt.xlabel('年份', fontsize=12)
plt.ylabel('人口（百万人）', fontsize=12)
plt.title('Logistic模型：解析解与数值解对比', fontsize=14)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('logistic_comparison.png', dpi=150)
plt.show()

参数敏感性分析

# 分析参数r和K对模型的影响
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# r的影响
t_plot = np.linspace(0, 25, 500)
r_values = [0.2, 0.3, 0.4, 0.5]
for r_val in r_values:
    N_val = K / (1 + 99.56 * np.exp(-r_val * t_plot))
    axes[0].plot(t_plot, N_val, linewidth=2, label=f'r = {r_val}')

axes[0].set_xlabel('时间（10年为单位）', fontsize=11)
axes[0].set_ylabel('人口（百万人）', fontsize=11)
axes[0].set_title('内禀增长率 r 的影响', fontsize=13)
axes[0].legend(fontsize=10)
axes[0].grid(True, alpha=0.3)

# K的影响
K_values = [300, 400, 500, 600]
for K_val in K_values:
    A_val = (K_val - N0) / N0
    N_val = K_val / (1 + A_val * np.exp(-0.3107 * t_plot))
    axes[1].plot(t_plot, N_val, linewidth=2, label=f'K = {K_val}')

axes[1].set_xlabel('时间（10年为单位）', fontsize=11)
axes[1].set_ylabel('人口（百万人）', fontsize=11)
axes[1].set_title('环境容纳量 K 的影响', fontsize=13)
axes[1].legend(fontsize=10)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('logistic_sensitivity.png', dpi=150)
plt.show()

置信区间估计

# 基于协方差矩阵估计参数的置信区间
perr = np.sqrt(np.diag(pcov))
print("参数估计的标准误差:")
print(f"  sigma_K = {perr[0]:.2f}")
print(f"  sigma_r = {perr[1]:.4f}")
print(f"  sigma_A = {perr[2]:.2f}")

# 95%置信区间
from scipy.stats import t as t_dist
n = len(population)
p = 3  # 参数个数
dof = n - p
t_critical = t_dist.ppf(0.975, dof)

print(f"\n95%置信区间 (自由度={dof}):")
print(f"  K: [{K_fit - t_critical*perr[0]:.2f}, {K_fit + t_critical*perr[0]:.2f}]")
print(f"  r: [{r_fit - t_critical*perr[1]:.4f}, {r_fit + t_critical*perr[1]:.4f}]")
print(f"  A: [{A_fit - t_critical*perr[2]:.2f}, {A_fit + t_critical*perr[2]:.2f}]")

# 预测区间（蒙特卡罗方法）
np.random.seed(42)
n_sim = 5000
params_sim = np.random.multivariate_normal(popt, pcov, n_sim)

t_future = np.linspace(0, 30, 300)
N_sims = np.array([logistic(t_future, *p) for p in params_sim
                   if p[0] > 0 and p[1] > 0 and p[2] > 0])

N_mean = np.mean(N_sims, axis=0)
N_lower = np.percentile(N_sims, 2.5, axis=0)
N_upper = np.percentile(N_sims, 97.5, axis=0)

plt.figure(figsize=(10, 6))
plt.fill_between(1790 + t_future * 10, N_lower, N_upper,
                 alpha=0.3, color='blue', label='95%预测区间')
plt.plot(1790 + t_future * 10, N_mean, 'b-', linewidth=2,
         label='均值预测')
plt.scatter(years, population, color='red', s=50,
            zorder=5, label='实际数据')
plt.xlabel('年份', fontsize=12)
plt.ylabel('人口（百万人）', fontsize=12)
plt.title('Logistic模型预测及置信区间', fontsize=14)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('logistic_confidence.png', dpi=150)
plt.show()

模型的推广与变体

含时滞的Logistic模型

实际种群的反馈往往存在时间延迟，引入时滞 \( \tau \) 后：

\[ \frac{dN}{dt} = rN(t)\left(1 - \frac{N(t-\tau)}{K}\right) \]

时滞模型可能产生振荡行为，当 \( r\tau > \pi/2 \) 时平衡点失稳。

广义Logistic模型（Richards模型）

引入形状参数 \( \nu \) 来调整曲线的不对称性：

\[ \frac{dN}{dt} = rN\left[1 - \left(\frac{N}{K}\right)^\nu\right] \]

当 \( \nu = 1 \) 时退化为标准Logistic模型。\( \nu < 1 \) 时拐点位于 \( K/2 \) 之下，\( \nu > 1 \) 时拐点位于 \( K/2 \) 之上。

Gompertz模型

另一种常用的S形增长模型：

\[ \frac{dN}{dt} = rN\ln\frac{K}{N} \]

其解析解为 \( N(t) = K\exp\left(-Ae^{-rt}\right) \)。Gompertz模型的拐点在 \( N = K/e \approx 0.368K \)，适用于不对称增长过程。

应用注意事项与局限性

适用场景

人口预测：区域人口的中长期预测
传染病传播：累计感染人数的增长趋势
产品扩散：新产品的市场渗透率（Bass模型的特殊情形）
生物培养：微生物在有限培养基中的生长
学习曲线：技能习得过程的建模
肿瘤生长：肿瘤体积随时间的变化

使用注意事项

环境容纳量的合理性：\( K \) 的估计高度依赖数据范围。若数据仅覆盖增长早期，\( K \) 的估计可能极不可靠。建议数据至少覆盖到拐点之后。
参数初值的选取：非线性拟合对初始值敏感。可先用线性化方法或三点法获得粗略估计作为初值。
时间尺度的选择：应根据具体问题选择合适的时间单位，避免参数值过大或过小导致的数值问题。
数据质量：异常值会严重影响拟合结果。建议先进行数据清洗和探索性分析。
模型选择：并非所有S形增长都适合用标准Logistic模型。若数据表现出明显的不对称性，应考虑Richards模型或Gompertz模型。

局限性

单一种群假设：Logistic模型仅考虑种内竞争，忽略了种间相互作用（捕食、共生等）。
环境恒定假设：模型假设 \( K \) 为常数，但实际中环境容纳量可能随时间变化（如技术进步、政策调整）。
连续性假设：模型假设种群变化是连续的，对于世代不重叠的种群，应使用离散Logistic映射： \[ N_{t+1} = rN_t(1 - N_t/K) \]
确定性假设：模型不包含随机性。对于小种群，随机涨落可能起决定性作用，需考虑随机微分方程模型。
对称性约束：标准Logistic曲线关于拐点严格对称，而许多实际增长过程是不对称的。
预测外推风险：模型在数据覆盖范围之外的预测可靠性迅速下降，尤其是对 \( K \) 的远期预测。

建模竞赛中的建议

在论文中明确阐述选择Logistic模型的依据和适用条件
进行残差分析，验证模型假设是否合理
进行参数敏感性分析，评估结果的稳健性
若条件允许，与其他模型（指数模型、Gompertz模型等）进行对比，论证所选模型的优越性
给出预测结果时附带置信区间或不确定性分析

小结

Logistic模型是从指数增长到有限增长的自然推广，通过引入阻滞因子捕获了资源约束下种群增长的本质特征。其解析解为S形曲线，具有明确的生物学和经济学解释。在数学建模实践中，掌握Logistic模型的推导、求解、参数估计和局限性分析，是处理各类增长问题的基础能力。

Keyboard shortcuts

数学建模指南