一种看待它的方式是说你有兴趣解决

$$y'(t) = f(y(t)), \qquad g(y(t)) = 0,$$ 并且您只近似地强制执行约束，从而解决 $$y'(t) = f(y(t)), \qquad g(y(t)) = O(\epsilon).$$

如果约束表现良好（比如说可微分，尽管这可能有点多），那么（通过隐函数定理）$g(y(t)) = O(\epsilon)$暗示$y(t)$“接近”真正的解决方案$y(t) = \tilde y(t) + \delta y(t)$和$g(\tilde y(t)) = 0$并且错误大约是

$$|\delta y(t)| = O(|\nabla g^{-1}|\epsilon).$$ （比如说，在最小二乘的意义上，可能奇异雅可比行列式的逆。）

因此，如果两者都近似$y$和真实的$\tilde y$满足相同的 ODE，则

$$(y-\tilde y)' = f(y) - f(\tilde y) \sim (\nabla f)(y-\tilde y).$$ 反之，求解近似系统就像求解扰动系统 $$y' = f(y) + (\nabla f)\delta y$$ 和$\delta y$是一个顺序上的不可知函数$|\nabla g^{-1}|\epsilon$.

所以重要的不是你使用的数值积分器，而是你的 ODE 对小扰动的敏感性。换句话说，如果积分器正确地求解了正确的方程，那么只有 ODE 的属性才是重要的。你关心的是错误是否$|y-\tilde y|$在较长时间内保持较小。

如果 ODE 是稳定的，即从相似初始条件开始的解在某种意义上保持接近，则误差的影响将是$O(|\nabla f||\nabla g^{-1}|\epsilon)$，所以希望足够小。如果 ODE 是向后稳定的，即，对于任何近似解，您都可以找到一个“附近”的初始条件，它会产生真正的解，那么从这个意义上说，误差会足够小。

在这两种情况下，问题在于 ODE 对扰动的敏感性，而不是数值方法。所以我想说检查任何近似解决方案就足够了$y$是否确实有“附近”的真解$\tilde y$, 这可以通过分析约束函数来完成。当然，您可以提出非常糟糕的约束。

通常，ODE 的微扰理论是一个众所周知的主题。

我知道的大多数 DAE 求解器都没有利用辛结构，但其中一个完全有可能做到这一点。

通常，如果您使用 DAE 求解器，则您使用的是隐式方法，或者可能是半隐式方法；代数变量是通过非线性求解确定的，我不知道您是如何解决的。我不知道您如何强制执行约束。如果您通过非线性求解修改位置，则半隐式 DAE 求解器可能具有竞争力；如果修改位置比非线性求解更便宜，那么 DAE 求解器可能更昂贵。完全隐式 DAE 求解器可能比您当前的方法更慢；然而，DAE 求解器背后有坚实的收敛理论（参见 Hairer 和 Wanner 的经典著作；Brenan 和 Petzold 的书，以及 Roswitha Marz 的书），这将使数值误差分析更加稳固。

也存在其他满足约束的方法。您可以尝试查看分子动力学社区如何在其机械系统中强制执行约束，例如通过拉格朗日乘法或专门的算法来强制执行诸如分子中原子的键长和键角之类的事情。听起来您正在做的是解决投影动态系统，并且您要小心该投影操作。我记得我的一位顾问讨论了投影如何影响稳定性，例如由于浮点运算导致的约束流形的扰动可能会产生不稳定的解决方案，因为基础 ODE 系统（投影之前）在关闭时不稳定多方面的方向。您还可以查看微分变分不等式。PETSc 中提供了变分不等式求解器和 DAE 求解器；您可以在 SUNDIALS 中找到隐式 DAE 求解器。

使用机械系统要注意的一件事是 DAE 的索引。（指数有很多定义；我只是假设“指数”的意思是“微分指数”，尽管我隐约意识到使用该定义存在缺陷。）机械系统可能有多种 DAE 公式，还有一些其中是高指数（实际上，这意味着指数大于 2）。大多数 DAE 求解器倾向于求解 index-1 或半显式 Hessenberg index-2 DAE；如果您的系统不是这种形式，您将需要一个专门的求解器（有一些专门用于机械系统的 index-3 求解器），或者您需要通过坐标更改等方式重新制定您的问题。