Symplectic SRK

[ChrisRackauckas]

https://www.cambridge.org/core/journals/communications-in-computational-physics/article/construction-of-symplectic-rungekutta-methods-for-stochastic-hamiltonian-systems/B09EE7BFE177D4B79DCB7C840F599319

[matteoettam09]

Hello, as said here, how can we best work on this? Has someone already started? Do you have some literature?

Thank you.

[ChrisRackauckas]

The literature is in the link above. The implementation would be very similar to the implementation of the SymplecticEuler method of OrdinaryDiffEq.jl. That's done via:

https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/src/perform_step/symplectic_perform_step.jl#L3-L76

It would need to be on the dynamical SDE form, which exists but is not documented yet since there's not specific solvers for it:

https://github.com/SciML/SciMLBase.jl/blob/master/src/problems/sde_problems.jl#L169-L194

How to add it is similar to the ODE case, which is showcased in the devdocs:

http://devdocs.sciml.ai/latest/contributing/adding_algorithms/

[matteoettam09]

Great, thank you, will try to do this.

[bennib22]

Hi @ChrisRackauckas, @matteoettam09 and I took a look at the paper and there's a thing that we don't understand. In the following I'm referring to equation (2.3) on p 240. As in a typical Runge-Kutta solver we need to calcuate the supporting values for each stage. In the paper they are called G_i_(0), G_i_(k) and G_i_(k). (Note: i is the subscript and (0) and (k) are superscripts. ) Suppose we want to calculate the first supporting value G_i(0) of the first stage. This means s= 1 and hence i=1. We are now calculating G_1(0). G_1(0) now depends on G_j_(0) where j will become 1 in the first stage. So we need G_1_(0) in order to calculate G_1_(0). How do we break out of this loop?

[ChrisRackauckas]

I just realized, we do have a method for Hamiltonian systems. You might want to give BAOBB a try.

https://github.com/SciML/StochasticDiffEq.jl/blob/v6.49.1/test/sde/sde_dynamical.jl https://github.com/SciML/StochasticDiffEq.jl/pull/397

We haven't documented it yet, but it follows the same pattern as the dynamical ODE function definition.

[matteoettam09]

Great, thank you, perfect timing, will now take a look.

I think I didn't specify that our diffusion matrix looks like this, so we are dealing with nondiagonal noise:

function σ_3bp(du,u,p,t)
    r = sqrt(u[1]^2 + u[2]^2)
    du[1,1] = 0.0
    du[1,2] = 0.0
    du[2,1] = 0.0
    du[2,2] = 0.0
    du[3,1] =  σ_r*u[1]/r
    du[3,2] = -σ_theta*u[2]
    du[4,1] =  σ_r*u[2]/r
    du[4,2] =  σ_theta*u[1]
end

[ChrisRackauckas]

Oh you'd need to take a look at the paper if it can be done for non-diagonal noise systems. I'm not sure that method extends.

[matteoettam09]

Ok, will do, thank you