向前欧拉法求解刚性常微分方程的数值稳定性边界实测偏差问题

1. 稳定区间推导验证 #

你的稳定区间推导过程完全正确：

• 该线性方程组的系数矩阵特征值为λ₁=-1、λ₂=-100，匹配向前欧拉法的绝对稳定判据|1 + hλ| ≤ 1。
• 分别代入两个特征值计算允许的步长范围：
  • 对λ=-1：|1 - h| ≤ 1 → 允许步长范围0 ≤ h ≤ 2
  • 对λ=-100：|1 - 100h| ≤1 → 允许步长范围0 ≤ h ≤ 1/50 = 0.02s
• 取两个范围的交集，全局稳定的步长上限确实是0.02s，推导逻辑没有错误。

2. 临界步长下数值解匹配度差的核心原因 #

需要明确绝对稳定的边界仅代表数值解不会无界发散，不代表能复现真实解的变化规律：

• 当h取临界值0.02s时，对应快衰减模式（y的方程，时间常数0.01s）的增长因子为1 + h*(-100) = -1，此时y的数值迭代关系为y_{n+1} = -y_n，会出现永久振荡、无衰减的行为，和真实解y(t)=0.1e^{-100t}的快速指数衰减特性完全不符：真实y在t=0.02s时已经衰减到初始值的13.5%，但数值解会在+0.1、-0.1之间来回跳变。
• 该快衰减模式会通过耦合项进入x的方程，导致x的数值解也被振荡的y值污染，自然和解析解偏差极大。
• 当h取1e-3时，对应快模式的增长因子为1 - 100*1e-3 = 0.9，和真实的指数衰减因子e^{-100*1e-3}≈0.905非常接近，能正确复现衰减行为，结果自然和解析解高度吻合。
• 对于刚性方程组，向前欧拉法的实际可用步长远小于稳定上限，只有当步长小到能正确模拟快模式的衰减速率时，才能得到精度合格的结果，不能卡着稳定边界选择步长。

3. 现有测试代码的存在问题 #

你的代码存在多处不完整/错误的地方，也会影响结果展示：

• 缺少必要的依赖导入：没有导入numpy，直接调用np.zeros会触发运行报错；
• 引用了未定义的变量：绘图用的时间序列t1、解析解数组s、第二个子图对象ax1都没有定义，代码本身无法正常运行；
• 绘图逻辑不完整：只创建了1行2列布局的第一个子图，第二个子图完全没有实现，还重复覆盖了ylim的设置；
• 提前定义的Euler函数没有在循环中实际调用，存在冗余代码。

核心解析解计算的参考实现片段：

import numpy as np
import matplotlib.pyplot as plt

tstop = 10.0
# 生成解析解参考曲线
t_analytic = np.linspace(0, tstop, 1000)
y_analytic = 0.1 * np.exp(-100 * t_analytic)
x_analytic = 0.1 * (98/99) * np.exp(-t_analytic) + 0.1/99 * np.exp(-100 * t_analytic)

内容的提问来源于stack exchange，提问作者user11607936