I am studying the transition of numerical methods from the classical setting to the manifold case in the context of geometric numerical integration. Precisely, I am studying the paper "GEOMETRIC INTEGRATION OF ORDINARY DIFFERENTIAL EQUATIONS ON MANIFOLDS" (https://link.springer.com/article/10.1023/A:1021989212020), and I am trying to reproduce the numerical examples. My struggle is with algorithm 3.2, the Symmetric projection method. Can anyone provide a Python script of the implementation of such method for the Rigid body example presented in the paper, considering the classical trapezoidal method as the one-step method in the second bullet point?
Per http://www.unige.ch/~hairer/preprints/symproj.pdf
It is important to take the same vector mu in (2.1) and (2.3).
where the equations 2.1, 2.2 and 2.3 describe the method as cited.
Which means that the method is implicit in mu. One can choose how exact the approximation for mu has to be, if it always has to be computed via Newton or if some fixed-point iteration is sufficient, if the Jacobian G needs to be computed in every step, etc.
If you are using a general-purpose non-linear solver, then it is better to implement an all-at-once approach, and not alternate the search directions.
Instead assemble state y1 and perturbation parameters mu into one input vector for the function serving as non-linear system in fsolve and return the residual of the RK step, here again the trapezoidal method, and the residuals of the first integrals, with the exact value being set to zero.
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]
This is more pseudo-code than directly executable. The passing of parameters can be made more explicit, using the args keyword parameter of fsolve.