Fix CG3

[mateuszbaran]

Initially implemented in #1 , it seems to not work and was cut from that PR:

"""
    CG3

A Crouch-Grossmann algorithm of second order for problems in the
[`ExplicitManifoldODEProblemType`](@ref) formulation. See tableau 6.1 of [^OwrenMarthinsen1999]:

    0     | 0
    3/4   | 3/4      0
    17/24 | 119/216  17/108  0
    ------------------------------
          | 13/51    -2/3    24/17

[^OwrenMarthinsen1999]:
    > B. Owren and A. Marthinsen, “Runge-Kutta Methods Adapted to Manifolds and Based on
    > Rigid Frames,” BIT Numerical Mathematics, vol. 39, no. 1, pp. 116–142, Mar. 1999,
    > doi: 10.1023/A:1022325426017.

"""
struct CG3{TM<:AbstractManifold,TR<:AbstractRetractionMethod} <: OrdinaryDiffEqAlgorithm
    manifold::TM
    retraction_method::TR
end

alg_order(::CG3) = 3

"""
    CG3Cache

Cache for [`CG3`](@ref).
"""
struct CG3Cache <: OrdinaryDiffEqMutableCache end

function alg_cache(
    alg::CG3,
    u,
    rate_prototype,
    uEltypeNoUnits,
    uBottomEltypeNoUnits,
    tTypeNoUnits,
    uprev,
    uprev2,
    f,
    t,
    dt,
    reltol,
    p,
    calck,
    ::Val{true},
)
    return CG3Cache()
end

function perform_step!(integrator, cache::CG3Cache, repeat_step = false)
    @unpack t, dt, uprev, u, f, p, alg = integrator
    M = alg.manifold

    k1 = f(u, p, t)
    c2h = (3 // 4) * dt
    c3h = (17 // 24) * dt
    a21h = (3 // 4) * dt
    a31h = (119 // 216) * dt
    a32h = (17 // 108) * dt
    b1 = (13 // 51) * dt
    b2 = (-2 // 3) * dt
    b3 = (24 // 17) * dt
    k2u = retract(M, u, k1 * a21h, alg.retraction_method)
    k2 = f(k2u, p, t + c2h)
    k1tk2u = f.f.operator_vector_transport(M, u, k1, k2u, p, t, t + c2h)
    k3u = retract(M, k2u, a31h * k1tk2u + a32h * k2)
    k3 = f(k3u, p, t + c3h)

    k2tu = f.f.operator_vector_transport(M, k2u, k2, u, p, t + c2h, t)
    k3tu = f.f.operator_vector_transport(M, k3u, k3, u, p, t + c3h, t)
    retract!(M, u, u, b1 * k1 + b2 * k2tu + b3 * k3tu, alg.retraction_method)

    return integrator.destats.nf += 3
end

function initialize!(integrator, cache::CG3Cache)
    integrator.kshortsize = 2
    integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)

    integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
    integrator.destats.nf += 1

    integrator.fsallast = zero.(integrator.fsalfirst)
    integrator.k[1] = integrator.fsalfirst
    integrator.k[2] = integrator.fsallast
    return nothing
end

[kellertuer]

I did not yet have time to investigate, but one thing that I had was that exp/log was sometimes less accurate compared to the ProjectionRetraction (for very small steps).

[mateuszbaran]

It's not this I think, it's also unstable on R^n and for relatively large steps (but small enough that CG2 works well). My guess is that the sequence of operations is wrong somewhere but I don't see where.

[kellertuer]

Ok, I will also take a look when I find time.

[mateuszbaran]

Thanks!

[kellertuer]

The first short question, where I am not sure in CG3 is: You are implementing an explicit variant of Eq. (2.3), right? Is that algorithm using retractions or the Lie group action? I am not sure there (since exp(X)p is not that explicitly mentioned what that is).

[mateuszbaran]

That's a relatively complex thing but yes, Eq. (2.3) is I think the most clear presentation of explicit Runge-Kutta frozen coefficient methods I've seen and what I was trying to implement here. There are no group actions as far as I can tell, these are flows of the frozen coefficients fields (as explained above Eq (2.1)), and these circles are, I think, compositions of these flows. Each e^{h b_r F^r_p} looks like a tangent moved from some point to p, and I've interpreted this composition as a linear combination in the tangent space (since all of these tangent vectors are at the same point, but maybe I'm wrong).

[kellertuer]

But why is there then a p at the end of the y_k+1 line and why was it not directly written as a sum? I think understanding exactly that might be the solution here (for CG2 there is just one e^... so there it is clear and sound). I thought that e^(hb_1F_p^1)p would be exp_p(hb_1F_p^1) and one would continue from there (though that would mean that hb_2F_p^2 is from another tangent space and now that I write that it seems strange, too?).

[kellertuer]

Here https://arc.aiaa.org/doi/pdf/10.2514/1.G004578, Eq (29) it reads like it is several exp_p(Xi) which are then combined with the Lie group action (\circ) ?

[mateuszbaran]

No, it doesn't look like Lie group action because above Eq. (2.3) from [^OwrenMarthinsen1999] there is no mention of Lie groups. It's a composition of flows of integral curves (see above Eq. (2.1)) and I think p means that we start from p. But apparently I don't understand how these flows are supposed to work in practice.

[mateuszbaran]

Fixed in https://github.com/JuliaManifolds/ManifoldDiffEq.jl/pull/1/commits/7450d69438c9e2c3e6d41b4609ccdb90e7a1ca72 .