Mathematica中NDSolve奇点/刚性系统错误的解决问询

问题概述 #

使用Wolfram Mathematica的NDSolve求解含变量分母的非线性微分方程组，结合Manipulate可视化结果时，触发NDSolve::ndsz警告：

NDSolve::ndsz: At t == ... step size is effectively zero; singularity or stiff system suspected

虽能生成图像，但警告提示数值求解过程中出现异常，需确认问题原因并修复以获得可靠解。

原因分析 #

1. 分母趋近于零：方程组中多处出现分母S[t] + 0.1*Ih[t]，当该值趋近于0时，分式项会趋向无穷大，导致数值求解器无法继续推进步长，触发奇点警告。
2. 刚性系统特性：方程组是非线性的，部分参数或初始值下，变量的变化速率差异极大（比如某变量快速衰减，另一变量缓慢变化），标准显式求解器难以适应这种剧烈的步长需求。
3. 初始值合理性：若初始值S0和Ih0过小，会直接让初始时刻的分母接近0，加速触发数值异常。

解决方案 #

1. 避免分母为零 #

给分母添加极小的偏移量，既不影响解的精度，又能彻底避免除以零的情况。例如将S[t] + 0.1*Ih[t]修改为S[t] + 0.1*Ih[t] + 1e-8。

2. 使用刚性系统专用求解器 #

在NDSolve中指定针对刚性系统的求解方法，推荐两种：

• Method -> "StiffnessSwitching"：自动检测刚性，在显式和隐式求解器间切换
• Method -> "BDF"：隐式向后差分格式，专门处理刚性微分方程

3. 调整初始值范围 #

限制S0和Ih0的下限，确保初始时刻分母不会过小。例如将初始值控件的下限从0改为0.01，避免初始状态就接近奇点。

4. 优化求解终止条件 #

添加WhenEvent，当变量达到无物理意义的阈值时自动停止求解，避免无效的数值计算。比如当S[t] < 1e-6时终止求解。

修正后的代码示例 #

Manipulate[
 Module[{plt1, plt2, plt3, sol, S0 = SS0, Ih0 = IhIh0, Y0 = YY0, denom},
  denom[t_] := S[t] + 0.1*Ih[t] + 1e-8; (* 添加偏移量避免分母为0 *)
  sol = NDSolve[{
      S'[t] == 3.5*S[t]*(1 - (S[t] + Ih[t])/50) - 1*S[t]*Ih[t] - (0.34*(S[t])^2*Y[t])/denom[t],
      Ih'[t] == 1*S[t]*Ih[t] - (2.2*(Ih[t])^2*Y[t])/denom[t] - 0.02*Ih[t],
      Y'[t] == (0.9*0.34*(S[t])^2*Y[t])/denom[t] - (0.02*2.2*(Ih[t])^2*Y[t])/denom[t] - 0.022*Y[t],
      S[0] == S0, Ih[0] == Ih0, Y[0] == Y0,
      WhenEvent[S[t] < 1e-6, "StopIntegration"] (* 自动终止无效求解 *)
    }, {S, Ih, Y}, {t, 0, 500},
    Method -> "StiffnessSwitching" (* 使用刚性切换求解器 *)
  ];
  plt1 = Plot[Evaluate[S[t] /. sol], {t, 0, 500}, PlotRange -> All, AspectRatio -> 1, PlotStyle -> {Red, Thick}, AxesLabel -> {"t", "S"}];
  plt2 = Plot[Evaluate[Ih[t] /. sol], {t, 0, 500}, PlotRange -> All, AspectRatio -> 1, PlotStyle -> {Green, Thick}, AxesLabel -> {"t", "Ih"}];
  plt3 = Plot[Evaluate[Y[t] /. sol], {t, 0, 500}, PlotRange -> All, AspectRatio -> 1, PlotStyle -> {Blue, Thick}, AxesLabel -> {"t", "Y"}];
  Show[plt1, plt2, plt3, ImageSize -> {300, 300}]
 ],
 Style["微分方程组 :", Bold],
 Style["S'= rS(1-(S+I)/K)-λSI-(pS² Y)/(S+αI)", Bold],
 Style["I'= λSI -(cI² Y)/(S+αI) - γI", Bold],
 Style["Y'= (δpS² Y)/(S+αI) - (ηcI² Y)/(S+αI) - dY", Bold],
 Delimiter,
 Style["参数", Bold, 10],
 Delimiter,
 Style["初始条件", Bold, 10],
 {{SS0, 4, "S0"}, 0.01, 10, .01, ImageSize -> Small, Appearance -> "Labeled"},
 {{IhIh0, 1, "Ih0"}, 0.01, 10, .01, ImageSize -> Small, Appearance -> "Labeled"},
 {{YY0, 6, "Y0"}, 0, 10, .01, ImageSize -> Small, Appearance -> "Labeled"},
 ControlPlacement -> Left, SynchronousUpdating -> False
]

注：代码中还将ParametricPlot替换为更高效的Plot[Evaluate[...]]，减少冗余计算。

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