wolfram-mathematica - Mathematica 的 NDSolve 正在快速找到一个刚性系统，试图求解大量微分方程

Question

我是一名物理学博士生，解决与圆柱形像素中二元液滴蒸发相关的润滑方程（非线性偏微分方程），以模拟形状演化。所以我将系统分成 N 个 ODE，每个代表 i 点的高度演变 (dh/dt)，写为有限差分，并同时使用 NDSolve 求解。对于单液滴，这非常有效，液滴蒸发干净。

然而，对于二元液滴，我在每个点包括 N 个额外的 ODE，用于组成分数演变 (dX/dt)。这也为 dh/dt 方程增加了一个新项（marangoni 应力）。NDSolve 立即说道：

在 t == 0.00029140763302667777`，步长实际上为零；
怀疑是奇点或僵硬系统。

我在这个 t 值处绘制了液滴高度，它在原点显示了一个狭窄的尖峰（明显爆炸，因此刚度）；绘制成分图在原点处显示了一个狭窄的直线下降（也在爆炸，只是负数。此外，物理学显然不应该允许成分分数低于 0）。

最后，两组方程的项都取决于成分梯度 dX/dr。如果将其设置为零并且我也将蒸发相互作用设置为零（意味着两种液体以相同的速率蒸发并且不会有 X 梯度），那么 X 在任何地方都应该没有变化，它应该减少到单一液体的情况（dX/dt = 0 且 dh/dt 不再取决于 X）。然而，该过程在 X 中引入了一些小的梯度，这会爆炸并导致相同的数值不稳定性。

我的问题是：NDSolve 有什么可能导致这种情况的吗？我已经研究过方程和离散化一百次，并且确信它是正确的。我还查看了 NDSolve 的文档，但没有找到任何对我有帮助的东西。它会在成分梯度中引入一个小的数值误差吗？

我可以在下面发布一个 MRE 代码，但它非常密集并且显然是用数学代码编写的（不能很好地转移到现实世界......）所以我不知道它有多大帮助。无论如何，感谢您阅读本文！

标签： wolfram-mathematica