我正在解决一个非线性 ODE 系统，其形式为

$Q_{nm} \ddot{y}_m + S_{nkl}\dot{y}_k\dot{y}_l +V_n = 0$

其中假设爱因斯坦求和，是因（复）变量， 、和都是因变量中的低阶多项式，，其中是的模型。这些方程与初始条件和一起完全说明了问题。$y_i$$Q_{nm}$$S_{nkl}$$V_n$$n =1,...,N$$N$$resolution$$y_i(0)$$\dot{y}_i(0)$

该方案使用（Fortran 接口）SUNDIALS 解决。特别是，由于非平凡系数，使用隐式微分代数求解器 (IDA) 来求解方程。$Q$

我已经成功地实现了这些方程，并且测试用例产生了一致的结果。计算时间受我们模型的分辨率$N$

现在，无论我多么精确地为某些物理相关场景指定初始条件，我发现我的系统中有一个虚假的增长，这导致方程变得僵硬，最终 Q 不再可逆并且积分失败。（暂时）缓解了这个问题，但是无论我可以合理地使多高，仍然会出现不稳定性。所以我想看看是否有人对平滑操作有任何建议，或者可能添加某种类型的粘度，以减轻这种不稳定性，同时保持结果物理相关。$N$$N$

作为一个例子，我把 IC 对应于已知解的方程。物理空间解与系数为的傅立叶级数的总和有关。在时间 t=0，我看到我的系统的数值解与我期望的结果之间的差异是好的（见图 1；注意这是针对的！），$y_i(t)$$N=16$

并且是 order，这是我在积分中的公差的顺序。请注意，解决方案是 order。现在，在时间 t=0.1，我所期望的和计算的解之间的差异不再那么微不足道，即我们看到差异大了 5 个数量级！这以这个指数速率增加，直到系统最终变得僵硬并停止积分。我不确定这是否是由于物理或数字引起的某种不稳定，但无论哪种情况，我都想摆脱它。$10^{-8}$$10^{-1}$

天真地，似乎集成步骤之间的某种平滑可能是有用的，但这有很大的缺点。IDA 使用内部雅可比行列式，因此在集成步骤之间冲动地更改解决方案不是一个好主意。

可能最吸引人的技术是向系统添加某种，以抑制或消除这些日益增长的不稳定性。解的不稳定部分增长最快，因此在这些方程中添加一些可以找出这一点的东西是理想的。我不清楚如何做到这一点，因此任何关于类似问题的建议或参考将不胜感激！$viscosity$

谢谢，

缺口