Relaxation Runge-Kutta

[ChrisRackauckas]

https://arxiv.org/abs/2001.04826

Seems like a good GSoC.

@ranocha

[ranocha]

Thanks for pinging me, @ChrisRackauckas.

The paper above investigates the setting of Hamiltonian problems. Just for reference: The earlier papers

• https://arxiv.org/abs/1905.09847
• https://arxiv.org/abs/1905.09129

describe the general procedure (and contain the basic proofs) and other applications. We are still working on extending and generalizing these results. That's why I haven't made a PR yet. But it definitely seems to be worthwhile to implement such a relatively general approach as used for conservative/Hamiltonian systems.

If a GSoC student is interested in this, I am of course willing to help and assist in mentoring, if you like. Since he is interested in open source projects (and relaxation Runge-Kutta methods, of course), I'll also ping @ketch to keep him up to date.

[ketch]

Yes, implementing those methods in JuliaDiffEq would be very nice. From a software point of view, they require some unusual things that would break most APIs:

• one or two user-defined functions need to be passed in and used inside the time step algorithm
• the step length dt that is actually taken is not the step length requested at the beginning of the step

[ranocha]

I think a good first step would be to implement relaxation methods for conservative systems (Hamiltonians, conserved entropies, ...) as callbacks, similar to https://docs.juliadiffeq.org/latest/features/callback_library/#Manifold-Conservation-and-Projection-1. Then, we can pass all necessary information to such a callback and can also update the current time of the integrator. This should be a rather general approach.

For dissipative systems, I don't see a good way to reuse the existing code, since we need such modifications as described by @ketch above.

[Biswajitghosh98]

@ChrisRackauckas @ranocha I was working on a low storage SSPRK algorithm mentioned in https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/issues/493 I'll try to submit the pull request as soon as possible for the same ( https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/compare/master...Biswajitghosh98:SSPRK53_H ) Can I start working on it alongside, and look up possible solutions for the bottlenecks mentioned by @ketch ?

[ranocha]

As mentioned above, it could be possible to start with an implementation of relaxation (for conservative systems at first) as callback, along the lines of

function relaxation!(integrator)
    told = integrator.tprev
    uold = integrator.uprev
    tnew = integrator.t
    unew = integrator.u

    # solve for the relaxation parameter γ

    @. unew = uold + γ * (unew - uold) # for inplace ops
    DiffEqBase.set_u!(integrator, unew)

    # This is the part that can be a bit nasty.
    # This check is a kind of a hack...
    # Without it, there can be errors because we do not hit
    # the final time anymore.
    if !(tnew ≈ top(integrator.opts.tstops))
        tγ = told + γ * (tnew - told)
        DiffEqBase.set_t!(integrator, tγ)
    end

    nothing
end

relaxation = DiscreteCallback((u,t,integrator) -> true, relaxation!, save_positions=(false,true))

[Biswajitghosh98]

@ranocha Thanks for helping me out. I'll start working on it, and will keep you updated .

[ranocha]

Another reference with generalizations of the theory and application to multistep methods: https://arxiv.org/abs/2003.03012