对称投影法（流形上的几何积分）

Question

我正在研究几何数值积分背景下数值方法从经典设置到流形情况的转变。准确地说，我正在研究论文“流形上常微分方程的几何积分”（https://link.springer.com/article/10.1023/A:1021989212020），并且我正在尝试重现数值示例。 我的困难在于算法 3.2，即对称投影方法。 任何人都可以为论文中提出的刚体示例提供这种方法的实现的 Python 脚本，将经典梯形方法视为第二个要点中的一步方法吗？

.

Answer 1

现在是。它产生了一个保留规范的解决方案，但 的行为并不符合预期（如 http://www.unige.ch/~hairer/preprints/symproj.pdf 中的图 4.1 所示）。

Answer 2

每 http://www.unige.ch/~hairer/preprints/symproj.pdf

重要的是要采取 (2.1) 和 (2.3) 中的向量 mu 相同。

其中方程 2.1、2.2 和 2.3 描述了所引用的方法。

这意味着该方法隐含在mu中。人们可以选择 mu 的近似值必须有多精确，如果它总是必须通过牛顿计算，或者如果某些定点迭代就足够了，如果雅可比行列式 G 需要在每一步中计算，等等。

---

如果您使用通用非线性求解器，那么最好实现一次性方法，而不是交替搜索方向。

而是将状态

y1

和扰动参数

mu

组装成

fsolve

中作为非线性系统的函数的一个输入向量，并返回 RK 步骤的残差，这里又是梯形方法，以及第一个积分，精确值设置为零。

def residuals(u):
    y1, mu = u[:sdim], u[sdim,:]
    y0mu = y0 + G(y0).T @ mu
    y1mu = y1 - G(y1).T @ mu
    f0mu = f(y0mu)
    f1mu = f(y1mu)
    restrap = (y1mu-y0mu)/dt - 0.5*(f0mu+f1mu)
    resfirst = g(y1)
    return [*restrap, *resfirst]

这更多是伪代码而不是直接可执行的代码。使用

args

的

fsolve

关键字参数可以使参数的传递更加明确。