数值积分约束动力学的简单方法？

Simple way to numerically integrate constrained dynamics?

例如，一个简单的单摆（字符串长度=1）可以描述为

mx'' = λ*2x

我的'' = -mg + λ*2y

x^2 + y^2 - 1 = 0

其中素数表示对w.r.t时间的导数，λ*2x和λ*2y是来自第三个方程的约束力，也就是约束。

当然，在这种情况下，通过选择弦的振幅θ作为唯一坐标，可以不使用λ。在这种情况下，我认为 Python(Scipy) 的 odeint 将是积分方程的最简单的免费方法。

然而，这样的坐标简化并不总是那么容易。

有没有一种简单的方法来积分像这样的约束 ODE（没有 Mathematica/Maple/Matlab）？性能和准确性对我来说都不重要，我只是不会快速检查结果。

（我不确定这个问题是否适合在这里或物理论坛，但对我来说，物理论坛似乎更多的是理论，而不是关于程序和数值的东西。）

非常感谢您的帮助！！

这种类型的系统称为微分代数方程。 netlib 中有标准求解器，例如 DASSL。

必须小心使用涵盖DAE的（微分）索引的方法。通常，求解器涵盖 index-1 或 index-2 系统。如果指数较高，则必须先使用符号指数降低方法才能应用数值方法。

摆的指数为3，这意味着系统的方程必须微分两次才能微分3次才能提取出一个普通的微分方程组。或者 2 个微分减少到一个显式索引 1 系统。