2026年数学建模美赛 常用模型算法 稳定性分析在数学建模中的应用：理论、方法与案例

2026美赛期间会持续更新相关内容，所有内容会发布到专栏内，会结合最新的chatgpt发布，只需订阅一次，赛后两天半价，内容达不到所有人预期，请勿盲目订阅！！！无论文！无论文！！！

摘要

稳定性分析是微分方程理论、控制理论和系统工程中的核心概念，在数学建模竞赛中具有广泛应用价值。本文系统阐述了稳定性分析的核心思想与数学模型，详细介绍了Lyapunov稳定性、结构稳定性等基本理论。通过典型场景分析、建模步骤解析、求解工具示例，并结合美国大学生数学建模竞赛简化案例，全面展示了稳定性分析在解决实际问题中的应用。最后，本文探讨了该方法的优缺点及改进方向，为数学建模竞赛参与者提供了实用的参考框架。

关键词：稳定性分析；Lyapunov方法；数学建模；微分方程；动力系统

1. 核心思想与数学模型

1.1 稳定性分析的基本哲学

稳定性分析研究的是系统在受到微小扰动后，能否恢复到原始状态或平衡状态的能力。这一概念源于对自然和工程系统行为的观察：许多系统在平衡点附近表现出"抗拒变化"的特性。在数学建模中，稳定性不是系统的固有属性，而是相对于特定平衡状态和特定类型扰动的特性。

稳定性分析的核心问题可以表述为：给定一个动力系统及其平衡点，当系统状态偏离平衡点时，系统是否会返回该平衡点？如果会，以何种方式返回？回答这些问题需要建立严格的数学框架。

1.2 数学基础与分类

1.2.1 动力系统的一般形式

考虑自治动力系统：

text

dx/dt = f(x), x ∈ ℝⁿ, f: ℝⁿ → ℝⁿ

其中x是状态变量，f是向量场。系统的平衡点（或不动点）是满足f(x*) = 0的点x*。

1.2.2 稳定性的严格定义

定义1（Lyapunov稳定性）：平衡点x*是Lyapunov稳定的，如果对于任意ε > 0，存在δ > 0，使得当‖x(0) – x*‖ < δ时，对于所有t ≥ 0，有‖x(t) – x*‖ < ε。

定义2（渐近稳定性）：平衡点x*是渐近稳定的，如果它是Lyapunov稳定的，并且存在η > 0，使得当‖x(0) – x*‖ < η时，lim_{t→∞} x(t) = x*。

定义3（指数稳定性）：平衡点x*是指数稳定的，如果存在常数α, β, η > 0，使得当‖x(0) – x*‖ < η时，有‖x(t) – x*‖ ≤ α‖x(0) – x*‖e^{-βt}。

1.2.3 线性系统的稳定性分析

对于线性系统：

text

dx/dt = Ax, A ∈ ℝ^{n×n}

平衡点（通常是原点）的稳定性完全由矩阵A的特征值决定：

• 如果A的所有特征值都有负实部，则原点是渐近稳定的

• 如果A至少有一个特征值有正实部，则原点是不稳定的

• 如果A的特征值都有非正实部，且零实部特征值对应的Jordan块都是一阶的，则原点是Lyapunov稳定但不渐近稳定

1.2.4 非线性系统的线性化方法

对于非线性系统，在平衡点x*附近进行Taylor展开：

text

f(x) = f(x*) + Df(x*)(x – x*) + O(‖x – x*‖²)

其中Df(x*)是Jacobian矩阵。线性化系统为：

text

d(Δx)/dt = Df(x*)Δx, 其中Δx = x – x*

根据Hartman-Grobman定理，如果Df(x)没有零实部特征值，则非线性系统在x附近的稳定性与线性化系统一致。

1.2.5 Lyapunov直接方法

Lyapunov直接方法避免了求解微分方程的困难，通过构造能量函数（Lyapunov函数）来判断稳定性。

定理（Lyapunov稳定性定理）：对于系统dx/dt = f(x)，f(0) = 0，如果存在定义在原点邻域D上的连续可微函数V: D → ℝ，满足：

则原点是Lyapunov稳定的。如果进一步有dV/dt < 0（对所有x ∈ D{0}），则原点是渐近稳定的。

1.2.6 结构稳定性与分岔理论

结构稳定性研究的是系统在参数变化时，定性行为是否保持不变。当参数通过临界值时，系统稳定性可能发生突变，这种现象称为分岔。

常见分岔类型包括：

• 鞍结分岔（saddle-node bifurcation）

• 跨临界分岔（transcritical bifurcation）

• 叉形分岔（pitchfork bifurcation）

• Hopf分岔（Hopf bifurcation）

1.3 离散系统的稳定性

对于离散动力系统：

text

x_{k+1} = f(x_k), x ∈ ℝⁿ

平衡点x满足x = f(x)。稳定性判据类似：平衡点渐近稳定当且仅当f在x处的Jacobian矩阵的所有特征值的模都小于1。

1.4 时滞系统的稳定性

时滞微分方程形式为：

text

dx/dt = f(x(t), x(t-τ))

时滞可能 destabilize 系统，使原本稳定的系统变得不稳定。分析方法包括特征方程法和Lyapunov-Krasovskii泛函法。

1.5 随机系统的稳定性

对于受随机扰动的系统：

text

dx = f(x)dt + g(x)dW

其中W是Wiener过程。需要定义随机稳定性概念，如均方稳定性、几乎必然稳定性等。

2. 适用场景与典型赛题类型

2.1 适用场景分析

稳定性分析适用于任何涉及"平衡"、"可持续性"、"鲁棒性"概念的问题场景：

2.2 数学建模竞赛中的典型赛题类型

2.2.1 连续动力系统类

示例：美国大学生数学建模竞赛2016年A题"热水浴缸温度模型"涉及热力学系统稳定性分析；2021年D题"音乐的影响力"可以建模为文化传播的动力系统。

特点：

• 问题可用微分方程描述

• 关注长期行为而非瞬时状态

• 需要确定参数阈值或临界条件

2.2.2 离散动力系统类

示例：中国大学生数学建模竞赛2018年B题"智能RGV的动态调度策略"可视为离散事件系统的稳定性问题。

特点：

• 系统状态在离散时间点变化

• 可能涉及迭代过程、递归关系

• 稳定性表现为收敛到固定点或周期轨道

2.2.3 时滞系统类

示例：网络舆情传播、供应链管理、具有反馈延迟的控制系统。

特点：

• 当前状态受过去状态影响

• 时滞可能导致振荡或失稳

• 需要特殊分析方法

2.2.4 随机扰动系统类

示例：金融风险评估、受随机干扰的生态系统、通信网络可靠性。

特点：

• 系统受随机因素影响

• 需要概率意义的稳定性

• 常用Ito随机微分方程描述

2.2.5 多稳定态与切换系统类

示例：气候系统突变、意识状态转换、多模态控制系统。

特点：

• 系统有多个可能的稳定状态

• 状态间可能存在切换

• 吸引域分析至关重要

2.3 稳定性分析在竞赛中的价值体现

3. 具体建模步骤与关键技巧

3.1 稳定性分析的标准流程

步骤1：问题识别与系统界定

• 明确系统中的状态变量、参数和控制输入

• 确定关心的平衡状态或参考轨迹

• 识别可能影响稳定性的扰动类型

步骤2：建立数学模型

• 选择适当的建模框架（连续/离散、确定/随机）

• 基于物理定律、经验关系或数据推导方程

• 验证模型的合理性和一致性

步骤3：寻找平衡点

• 求解f(x) = 0（连续系统）或x = f(x)（离散系统）

• 注意可能存在多个平衡点

• 对复杂的隐式方程使用数值方法

步骤4：线性化分析

• 计算Jacobian矩阵

• 分析特征值分布

• 判断线性化系统的稳定性

步骤5：非线性分析（如需要）

• 当线性化方法不适用时（特征值有零实部）

• 构造Lyapunov函数

• 使用中心流形定理简化分析

步骤6：参数影响分析

• 研究关键参数变化对稳定性的影响

• 绘制稳定性边界（分岔图）

• 确定临界参数值

步骤7：数值验证与模拟

• 使用数值积分验证理论分析

• 模拟不同初始条件和参数下的系统行为

• 可视化相图、时间序列等

步骤8：结果解释与应用

• 将数学结论转化为实际问题解答

• 提出维持稳定或避免失稳的建议

• 讨论模型的局限性和改进方向

3.2 关键技巧与常见陷阱

3.2.1 Lyapunov函数构造技巧

能量类比法：在物理系统中，总能量（动能+势能）往往是天然的Lyapunov函数候选。

变量梯度法：设V(x) = ∫₀ˣ [f(s)]ᵀP ds，其中P为正定矩阵，适当选择P可使dV/dt负定。

平方和形式：尝试V(x) = xᵀPx，这是最常用的二次型Lyapunov函数。

Krasovskii方法：对于系统dx/dt = f(x)，考虑V(x) = f(x)ᵀf(x)。

变量部分法：对于复杂系统，分别构造各子系统的Lyapunov函数，再组合。

3.2.2 处理零实部特征值的技巧

当线性化矩阵有零实部特征值时，需要高阶分析：

3.2.3 数值稳定性分析的注意事项

3.2.4 多尺度系统稳定性分析

对于快慢变量分离的系统：

text

εdx/dt = f(x,y)
dy/dt = g(x,y)

其中ε << 1。可采用奇异摄动理论，分别分析快慢子系统。

3.3 稳定性判据与定理总结

Routh-Hurwitz判据：不计算特征值，直接由多项式系数判断稳定性。

Nyquist判据：频率域稳定性判据，特别适用于控制系统。

Circle判据：处理不确定非线性系统的鲁棒稳定性。

Small-gain定理：互联系统的输入-输出稳定性。

LaSalle不变原理：Lyapunov函数导数半负定时仍可证明渐近稳定性。

4. 常用求解工具/代码示例

4.1 MATLAB/Python工具包介绍

4.1.1 MATLAB相关工具

• eig()：计算矩阵特征值

• lyap()：求解Lyapunov方程

• ode45/ode15s：常微分方程数值求解

• Control System Toolbox：控制系统稳定性分析

• Symbolic Math Toolbox：符号计算

4.1.2 Python生态系统

• NumPy/SciPy：数值计算和线性代数

• SymPy：符号计算

• Matplotlib：结果可视化

• Control Systems Library (python-control)：控制理论工具

4.2 代码示例合集

示例1：线性系统稳定性分析（Python）

python

import numpy as np
import matplotlib.pyplot as plt
from scipy import linalg

def analyze_linear_stability(A):
"""分析线性系统dx/dt=Ax的稳定性"""
eigenvalues, eigenvectors = linalg.eig(A)

print("特征值:", eigenvalues)

# 判断稳定性
max_real = np.max(np.real(eigenvalues))
if max_real < 0:
stability = "渐近稳定"
elif max_real <= 1e-10: # 考虑数值误差
# 检查零实部特征值的代数重数
zero_eig = np.abs(np.real(eigenvalues)) < 1e-10
if np.all(np.imag(eigenvalues[zero_eig]) == 0):
stability = "临界稳定"
else:
stability = "需进一步分析"
else:
stability = "不稳定"

print(f"系统稳定性: {stability}")

# 可视化特征值分布
plt.figure(figsize=(8, 6))
plt.scatter(np.real(eigenvalues), np.imag(eigenvalues),
c='red', s=100, marker='o')
plt.axvline(x=0, color='k', linestyle='–', alpha=0.5)
plt.axhline(y=0, color='k', linestyle='–', alpha=0.5)
plt.xlabel('实部')
plt.ylabel('虚部')
plt.title('特征值分布图')
plt.grid(True, alpha=0.3)
plt.show()

return eigenvalues, stability

# 示例矩阵
A = np.array([[-2, 1],
[0.5, -1]])
eigenvalues, stability = analyze_linear_stability(A)

示例2：非线性系统Lyapunov函数构造（MATLAB）

matlab

% 定义系统
syms x1 x2 real
f1 = -x1 + 2*x1*x2;
f2 = -x2 + x1^2 – x2^2;
f = [f1; f2];

% 尝试二次型Lyapunov函数 V = x'*P*x
P = sym('P', [2, 2]);
assume(P, 'real');
V = [x1, x2] * P * [x1; x2];

% 计算V沿系统轨迹的导数
gradV = gradient(V, [x1, x2]);
V_dot = gradV(1)*f1 + gradV(2)*f2;

% 将V_dot表示为二次型形式
V_dot_quad = collect(V_dot, [x1, x2]);

% 我们希望V正定，V_dot负定
% 可以通过求解矩阵不等式寻找合适的P
% 这里展示符号推导过程
disp('Lyapunov函数:')
pretty(V)
disp('其导数:')
pretty(V_dot)

% 数值求解示例：使用特定P
P_num = eye(2); % 单位矩阵
V_num = [x1, x2] * P_num * [x1; x2];
V_dot_num = simplify(jacobian(V_num, [x1, x2]) * f);

disp('数值Lyapunov函数导数:')
pretty(V_dot_num)

示例3：分岔分析（Python）

python

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

def saddle_node_bifurcation(x, t, r):
"""鞍结分岔标准形式"""
return r + x**2

def bifurcation_analysis():
"""参数分岔分析"""
# 参数范围
r_values = np.linspace(-2, 2, 100)

# 存储平衡点
equilibria = []

for r in r_values:
# 寻找平衡点: r + x^2 = 0
if r < 0:
# 两个平衡点
x1 = np.sqrt(-r)
x2 = -np.sqrt(-r)
equilibria.append((r, x1))
equilibria.append((r, x2))
elif r == 0:
# 一个平衡点（退化）
equilibria.append((r, 0))
else:
# 无实平衡点
pass

# 转换为数组以便绘图
equilibria = np.array(equilibria)

# 绘制分岔图
plt.figure(figsize=(10, 6))

if len(equilibria) > 0:
plt.plot(equilibria[:, 0], equilibria[:, 1],
'b-', linewidth=2, label='平衡点')

# 稳定性分析
# 对于每个平衡点，计算导数: f'(x) = 2x
# f'(x) < 0 稳定，f'(x) > 0 不稳定

# 标记稳定性
if len(equilibria) > 0:
stable_mask = 2 * equilibria[:, 1] < 0
unstable_mask = ~stable_mask

plt.scatter(equilibria[stable_mask, 0], equilibria[stable_mask, 1],
color='green', s=50, label='稳定', zorder=5)
plt.scatter(equilibria[unstable_mask, 0], equilibria[unstable_mask, 1],
color='red', s=50, label='不稳定', zorder=5)

plt.axvline(x=0, color='k', linestyle='–', alpha=0.5, label='分岔点')
plt.xlabel('参数 r')
plt.ylabel('平衡点 x*')
plt.title('鞍结分岔图')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

bifurcation_analysis()

示例4：时滞系统稳定性（MATLAB）

matlab

% 时滞微分方程稳定性分析示例
% 系统: dx/dt = -x(t) + a*x(t-τ)

clear; close all;

% 参数设置
a_values = 0:0.01:2; % 参数a的范围
tau = 1; % 固定时滞

% 存储稳定性结果
stability = zeros(size(a_values));

% 对每个参数值分析稳定性
for i = 1:length(a_values)
a = a_values(i);

% 特征方程: λ + 1 – a*exp(-λτ) = 0
% 使用数值方法寻找根

% 定义特征方程函数
char_eq = @(lambda) lambda + 1 – a*exp(-lambda*tau);

% 在复平面上搜索根
% 简单判据：当a<1时稳定，但严格分析需要计算特征根

% 这里使用简化判据：对于这个特定系统，稳定条件是|a|<1
if abs(a) < 1
stability(i) = 1; % 稳定
else
stability(i) = 0; % 不稳定
end
end

% 绘制稳定性区域
figure;
plot(a_values, stability, 'b-', 'LineWidth', 2);
xlabel('参数 a');
ylabel('稳定性 (1=稳定, 0=不稳定)');
title(sprintf('时滞系统稳定性 vs 参数a (τ=%g)', tau));
grid on;
ylim([-0.1, 1.1]);

% 标记临界点
hold on;
plot([1, 1], [0, 1], 'r–', 'LineWidth', 1.5);
text(1.05, 0.5, '临界点 a=1', 'Color', 'red');

% 数值模拟验证
figure;
subplot(2,1,1);
% 稳定情况模拟
a_stable = 0.5;
[t_stable, x_stable] = dde23(@(t,x,Z) -x + a_stable*Z, tau, …
@(t) 0.5, [0, 20]);
plot(t_stable, x_stable, 'g-', 'LineWidth', 2);
title(sprintf('稳定情况: a=%g', a_stable));
xlabel('时间 t');
ylabel('x(t)');
grid on;

subplot(2,1,2);
% 不稳定情况模拟
a_unstable = 1.5;
[t_unstable, x_unstable] = dde23(@(t,x,Z) -x + a_unstable*Z, tau, …
@(t) 0.5, [0, 20]);
plot(t_unstable, x_unstable, 'r-', 'LineWidth', 2);
title(sprintf('不稳定情况: a=%g', a_unstable));
xlabel('时间 t');
ylabel('x(t)');
grid on;

示例5：随机系统稳定性（Python）

python

import numpy as np
import matplotlib.pyplot as plt

def stochastic_system_simulation():
"""随机微分方程模拟与稳定性分析"""
np.random.seed(42)

# 系统参数
mu = -0.5 # 漂移系数
sigma = 0.3 # 扩散系数
x0 = 1.0 # 初始条件
T = 10.0 # 总时间
dt = 0.01 # 时间步长
n_steps = int(T/dt)

# 数值模拟（Euler-Maruyama方法）
t = np.linspace(0, T, n_steps+1)
x = np.zeros(n_steps+1)
x[0] = x0

for i in range(n_steps):
dW = np.random.normal(0, np.sqrt(dt)) # Wiener增量
x[i+1] = x[i] + mu*x[i]*dt + sigma*x[i]*dW

# 理论分析：几何布朗运动的均方稳定性
# dx = mu*x dt + sigma*x dW
# 均方稳定条件: 2*mu + sigma^2 < 0

stability_condition = 2*mu + sigma**2
if stability_condition < 0:
stability = "均方稳定"
else:
stability = "均方不稳定"

# 绘制结果
plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
plt.plot(t, x, 'b-', linewidth=1.5, label='样本路径')
plt.xlabel('时间 t')
plt.ylabel('状态 x(t)')
plt.title(f'随机系统模拟\\nmu={mu}, sigma={sigma}')
plt.legend()
plt.grid(True, alpha=0.3)

plt.subplot(1, 2, 2)
# 多次模拟观察统计特性
n_simulations = 50
x_all = np.zeros((n_simulations, n_steps+1))

for j in range(n_simulations):
x_temp = np.zeros(n_steps+1)
x_temp[0] = x0
for i in range(n_steps):
dW = np.random.normal(0, np.sqrt(dt))
x_temp[i+1] = x_temp[i] + mu*x_temp[i]*dt + sigma*x_temp[i]*dW
x_all[j] = x_temp
plt.plot(t, x_temp, 'gray', alpha=0.3)

# 均值±标准差
mean_x = np.mean(x_all, axis=0)
std_x = np.std(x_all, axis=0)

plt.plot(t, mean_x, 'r-', linewidth=2, label='均值')
plt.fill_between(t, mean_x – std_x, mean_x + std_x,
alpha=0.3, color='red', label='±1标准差')

plt.xlabel('时间 t')
plt.ylabel('状态 x(t)')
plt.title(f'多次模拟统计特性\\n稳定性: {stability}')
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"稳定性分析:")
print(f" 参数: mu = {mu}, sigma = {sigma}")
print(f" 稳定性条件: 2*mu + sigma^2 = {stability_condition}")
print(f" 系统状态: {stability}")

stochastic_system_simulation()

5. 一个完整的美赛简化案例（含问题、建模、求解与分析）

5.1 问题背景：社交媒体信息传播的动态与控制

问题描述：在社交媒体平台上，信息的传播速度和范围受多种因素影响。假设某平台上有两种竞争性信息（如真实新闻和虚假新闻）同时传播。平台管理者希望了解：

5.2 数学模型建立

5.2.1 基本假设

5.2.2 模型方程

基于竞争性SIR模型，建立如下方程：

text

dS/dt = μN – β₁SI₁/N – β₂SI₂/N – μS
dI₁/dt = β₁SI₁/N – γ₁I₁ – μI₁ + α₁I₂ – δ₁I₁I₂/N + u₁(t)
dI₂/dt = β₂SI₂/N – γ₂I₂ – μI₂ + α₂I₁ – δ₂I₁I₂/N + u₂(t)
dR/dt = γ₁I₁ + γ₂I₂ – μR

其中：

• βᵢ：信息i的传播率

• γᵢ：信息i的失去兴趣率

• μ：用户进入/离开率

• αᵢ：从另一种信息转换到信息i的率

• δᵢ：竞争导致的传播抑制系数

• uᵢ(t)：平台控制输入（如可见度调整）

约束条件：S + I₁ + I₂ + R = N

5.3 稳定性分析

5.3.1 简化与无量纲化

令s = S/N, i₁ = I₁/N, i₂ = I₂/N, r = R/N，系统简化为：

text

ds/dt = μ – β₁si₁ – β₂si₂ – μs
di₁/dt = β₁si₁ – (γ₁+μ)i₁ + α₁i₂ – δ₁i₁i₂ + u₁(t)
di₂/dt = β₂si₂ – (γ₂+μ)i₂ + α₂i₁ – δ₂i₁i₂ + u₂(t)

由于s + i₁ + i₂ + r = 1，可以消去一个变量。

5.3.2 平衡点分析

先考虑无控制情况（u₁=u₂=0）。设平衡点为(s, i₁, i₂*)，满足：

text

0 = μ – β₁s*i₁* – β₂s*i₂* – μs*
0 = β₁s*i₁* – (γ₁+μ)i₁* + α₁i₂* – δ₁i₁*i₂*
0 = β₂s*i₂* – (γ₂+μ)i₂* + α₂i₁* – δ₂i₁*i₂*

平衡点1：信息灭绝点 (s=1, i₁=0, i₂=0)

平衡点2：仅真实信息存在 (s₁, i₁, 0)，其中：

text

s₁* = (γ₁+μ)/β₁
i₁* = μ(β₁ – γ₁ – μ)/[β₁(γ₁+μ)]

存在条件：β₁ > γ₁ + μ（基本再生数R₀₁ > 1）

平衡点3：仅虚假信息存在 (s₂, 0, i₂)

平衡点4：共存平衡点 (s, i₁, i₂*)，需数值求解

5.3.3 局部稳定性分析

计算Jacobian矩阵：

text

J = [ -β₁i₁-β₂i₂-μ -β₁s -β₂s
β₁i₁ β₁s-(γ₁+μ)-δ₁i₂ α₁-δ₁i₁
β₂i₂ α₂-δ₂i₂ β₂s-(γ₂+μ)-δ₂i₁ ]

在信息灭绝点(1,0,0)：

text

J(1,0,0) = [ -μ -β₁ -β₂
0 β₁-(γ₁+μ) α₁
0 α₂ β₂-(γ₂+μ) ]

特征值为：λ₁ = -μ，λ₂ = β₁-(γ₁+μ)，λ₃ = β₂-(γ₂+μ)

因此，当R₀₁ = β₁/(γ₁+μ) < 1且R₀₂ = β₂/(γ₂+μ) < 1时，灭绝点局部渐近稳定。

在仅真实信息平衡点(s₁, i₁, 0)：
稳定性条件更复杂，但可以证明当真实信息的竞争优势足够大时稳定。

5.3.4 基本再生数与阈值现象

定义信息i的基本再生数：

text

R₀ᵢ = βᵢ/(γᵢ+μ)

这表示一个传播者在全易感人群中能产生的新传播者数量。

• 若R₀ᵢ < 1，信息i无法持续传播

• 若R₀ᵢ > 1，信息i可能持续传播

5.4 控制策略设计与稳定性

5.4.1 控制目标

设计控制输入u₁(t), u₂(t)使得：

5.4.2 线性反馈控制设计

在期望平衡点附近线性化，设计状态反馈：

text

u₁(t) = -k₁₁(i₁ – i₁*) – k₁₂(i₂ – i₂*)
u₂(t) = -k₂₁(i₁ – i₁*) – k₂₂(i₂ – i₂*)

通过极点配置或LQR方法确定增益矩阵K。

5.4.3 Lyapunov-based控制设计

构造Lyapunov函数：

text

V(i₁, i₂) = ½(i₁ – i₁*)² + ½(i₂ – i₂*)²

设计控制律使dV/dt负定。

5.5 数值模拟与结果分析

python

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

def information_spread_model(state, t, params, control=False):
"""竞争性信息传播模型"""
s, i1, i2 = state
mu, beta1, beta2, gamma1, gamma2, alpha1, alpha2, delta1, delta2 = params

# 控制输入（简单比例控制）
if control:
i1_target, i2_target = 0.3, 0.1 # 期望平衡点
k1, k2 = 0.5, 0.5 # 控制增益
u1 = -k1 * (i1 – i1_target)
u2 = -k2 * (i2 – i2_target)
else:
u1, u2 = 0, 0

dsdt = mu – beta1*s*i1 – beta2*s*i2 – mu*s
di1dt = beta1*s*i1 – (gamma1+mu)*i1 + alpha1*i2 – delta1*i1*i2 + u1
di2dt = beta2*s*i2 – (gamma2+mu)*i2 + alpha2*i1 – delta2*i1*i2 + u2

# 保持总和为1（近似）
total = s + i1 + i2
if total > 1:
scale = 1/total
dsdt *= scale
di1dt *= scale
di2dt *= scale

return [dsdt, di1dt, di2dt]

def simulate_scenarios():
"""模拟不同场景"""
# 参数设置
mu = 0.01 # 用户更新率
beta1, beta2 = 0.5, 0.6 # 传播率（虚假信息传播更快）
gamma1, gamma2 = 0.1, 0.08 # 失去兴趣率
alpha1, alpha2 = 0.05, 0.03 # 信息转换率
delta1, delta2 = 0.2, 0.1 # 竞争抑制

params = (mu, beta1, beta2, gamma1, gamma2, alpha1, alpha2, delta1, delta2)

# 初始条件：大多数易感，少量传播者
s0, i10, i20 = 0.95, 0.03, 0.02
initial_state = [s0, i10, i20]

# 时间点
t = np.linspace(0, 100, 1000)

# 场景1：无控制
sol_no_control = odeint(information_spread_model, initial_state, t,
args=(params, False))

# 场景2：有控制
sol_with_control = odeint(information_spread_model, initial_state, t,
args=(params, True))

# 计算基本再生数
R0_truth = beta1/(gamma1+mu)
R0_fake = beta2/(gamma2+mu)

# 可视化
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

# 无控制情况
axes[0,0].plot(t, sol_no_control[:,0], 'b-', label='易感者S', linewidth=2)
axes[0,0].plot(t, sol_no_control[:,1], 'g-', label='真实信息传播者', linewidth=2)
axes[0,0].plot(t, sol_no_control[:,2], 'r-', label='虚假信息传播者', linewidth=2)
axes[0,0].set_xlabel('时间')
axes[0,0].set_ylabel('比例')
axes[0,0].set_title('无控制情况')
axes[0,0].legend()
axes[0,0].grid(True, alpha=0.3)

# 有控制情况
axes[0,1].plot(t, sol_with_control[:,0], 'b-', label='易感者S', linewidth=2)
axes[0,1].plot(t, sol_with_control[:,1], 'g-', label='真实信息传播者', linewidth=2)
axes[0,1].plot(t, sol_with_control[:,2], 'r-', label='虚假信息传播者', linewidth=2)
axes[0,1].set_xlabel('时间')
axes[0,1].set_ylabel('比例')
axes[0,1].set_title('有控制情况')
axes[0,1].legend()
axes[0,1].grid(True, alpha=0.3)

# 相图（无控制）
axes[0,2].plot(sol_no_control[:,1], sol_no_control[:,2], 'b-', alpha=0.7)
axes[0,2].scatter(sol_no_control[0,1], sol_no_control[0,2],
color='green', s=100, marker='o', label='起点')
axes[0,2].scatter(sol_no_control[-1,1], sol_no_control[-1,2],
color='red', s=100, marker='s', label='终点')
axes[0,2].set_xlabel('真实信息传播者')
axes[0,2].set_ylabel('虚假信息传播者')
axes[0,2].set_title('相图 (无控制)')
axes[0,2].legend()
axes[0,2].grid(True, alpha=0.3)

# 控制效果对比
axes[1,0].plot(t, sol_no_control[:,1], 'g–', label='真实信息(无控制)', linewidth=2)
axes[1,0].plot(t, sol_with_control[:,1], 'g-', label='真实信息(有控制)', linewidth=2)
axes[1,0].set_xlabel('时间')
axes[1,0].set_ylabel('比例')
axes[1,0].set_title('真实信息传播者对比')
axes[1,0].legend()
axes[1,0].grid(True, alpha=0.3)

axes[1,1].plot(t, sol_no_control[:,2], 'r–', label='虚假信息(无控制)', linewidth=2)
axes[1,1].plot(t, sol_with_control[:,2], 'r-', label='虚假信息(有控制)', linewidth=2)
axes[1,1].set_xlabel('时间')
axes[1,1].set_ylabel('比例')
axes[1,1].set_title('虚假信息传播者对比')
axes[1,1].legend()
axes[1,1].grid(True, alpha=0.3)

# 信息优势比
ratio_no_control = sol_no_control[:,1] / (sol_no_control[:,2] + 1e-10)
ratio_with_control = sol_with_control[:,1] / (sol_with_control[:,2] + 1e-10)

axes[1,2].plot(t, ratio_no_control, 'b–', label='无控制', linewidth=2)
axes[1,2].plot(t, ratio_with_control, 'b-', label='有控制', linewidth=2)
axes[1,2].axhline(y=1, color='r', linestyle='–', alpha=0.5, label='平衡线')
axes[1,2].set_xlabel('时间')
axes[1,2].set_ylabel('真实信息/虚假信息')
axes[1,2].set_title('信息优势比')
axes[1,2].legend()
axes[1,2].grid(True, alpha=0.3)
axes[1,2].set_yscale('log')

plt.suptitle(f'社交媒体信息传播稳定性分析\\nR0(真实)={R0_truth:.2f}, R0(虚假)={R0_fake:.2f}',
fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

# 稳定性分析总结
print("=== 稳定性分析结果 ===")
print(f"基本再生数:")
print(f" 真实信息: R0₁ = {R0_truth:.3f} {'(可持续传播)' if R0_truth > 1 else '(会自然消失)'}")
print(f" 虚假信息: R0₂ = {R0_fake:.3f} {'(可持续传播)' if R0_fake > 1 else '(会自然消失)'}")

print(f"\\n最终状态 (无控制):")
print(f" 真实信息传播者: {sol_no_control[-1,1]:.4f}")
print(f" 虚假信息传播者: {sol_no_control[-1,2]:.4f}")
print(f" 优势比: {ratio_no_control[-1]:.2f}")

print(f"\\n最终状态 (有控制):")
print(f" 真实信息传播者: {sol_with_control[-1,1]:.4f}")
print(f" 虚假信息传播者: {sol_with_control[-1,2]:.4f}")
print(f" 优势比: {ratio_with_control[-1]:.2f}")

print(f"\\n控制效果:")
print(f" 真实信息增加: {((sol_with_control[-1,1]-sol_no_control[-1,1])/sol_no_control[-1,1]*100):.1f}%")
print(f" 虚假信息减少: {((sol_no_control[-1,2]-sol_with_control[-1,2])/sol_no_control[-1,2]*100):.1f}%")

simulate_scenarios()

5.6 案例总结与竞赛应用启示

通过这个简化案例，我们展示了：

竞赛应用建议：

• 在美赛中，类似问题可以扩展到时滞、随机、网络结构等因素

• 稳定性分析能提供深刻的洞察，而不仅仅是数值结果

• 结合控制理论可以设计优化策略

• 可视化结果能有效展示模型行为

6. 该方法的优缺点及改进方向

6.1 稳定性分析方法的优点

6.2 稳定性分析方法的局限性

6.3 当前研究前沿与改进方向

6.3.1 计算Lyapunov函数的新方法

6.3.2 全局稳定性分析

6.3.3 鲁棒稳定性分析

6.3.4 数据驱动的稳定性分析

6.3.5 网络化系统的稳定性

6.4 数学建模竞赛中的应用建议

6.5 稳定性分析的未来展望

随着计算能力的提升和数学理论的发展，稳定性分析正朝着以下方向发展：