Projecting u during propagation

I am solving a set of equations of motion that I know is only approximate and will lead to unphysical results. After each propagation step, I want to project u back onto a physical subspace.

I was wondering if something like this is possible using Tsit5 of the DifferentialEquations module. I have heard that one should not change u manually during propagation and wanted to make sure that what I am doing is not leading to unexpected results. For example, I am unsure whether Tsit5() uses the projected values of u to calculate the next update or whether it uses the values before the projection.

My code looks something like this:

function always(u, t, integrator)
    return true
end

function projection_func!(slv::Solve.Instance)
    return function (integrator)
        # modify u here
        project!(
            integrator.u, 
            max_steps = slv.n_project_max, 
        )
    end
end

function solve_projected(slv, uâ‚€, p)
    cb = DiscreteCallback(always, projection_func!(slv))

    # Do time evolution
    prob = ODEProblem(eoms!, uâ‚€, slv.tspan, p)
    sol = solve(
        prob,
        Tsit5(),
        callback = cb,
        dtmax = slv.dtmax,
        maxiters = slv.maxiters,
        reltol = slv.reltol,
        abstol = slv.abstol,
        saveat = slv.saveat,
    )
    return sol
end

Check ManifoldProjection in DiffEqCallbacks.jl.

Or for a lot more detail on this topic:

1 Like

Thanks for the helpful ressources! The manifold projection is definitely something I will consider.

But coming back to my original question: Will my approach with Tsit5() work in principle? Or is what I am doing “evil” in the sense that the solver does not expect it?

Nothing wrong there, just set save_positions correctly and know that the projection will mean that you only can get linear interpolations

1 Like