Bug in DiffEqPhysics ?

[KlausC]

What did I wrong?

Trying to solve one of the simplest dynamical systems using the Hamiltonian method with DiffEqPhysics failed.

We have H(p, q) = p^2/20 (a inert body of 10 kg) p0 = 1.0 (initial momentum = 1 kg * m / s) q0 = 0.0 (initial speed = 0 m/s)

equations of motion: dq/dt = ∂H/∂p = p / 10 dp/dt = - ∂H/∂q = 0

solutions: p(t) = 1.0 q(t) = 0.1 * t

That should be reproducible using the package, but surprisingly delivers two wrong solutions with 2 different solver algorithms.

julia> using DifferentialEquations
julia> H(p, q, param) = p^2/20
H (generic function with 1 method)
julia> p0, q0 = 1.0, 0.0;

julia> prob = HamiltonianProblem(H, p0, q0, (0., 1.))
ODEProblem with uType ArrayPartition{Float64,Tuple{Float64,Float64}} and tType Float64. In-place: false
timespan: (0.0, 1.0)
u0: 1.00.0

julia> sol = solve(prob, SofSpa10(), dt=0.001);
julia> sol(1.0)
(1.0, 0.13545807513475394)

julia> sol = solve(prob, dt=0.001);
julia> sol(1.0)
(1.0, 1.0)

[ChrisRackauckas]

Are the inputs flipped?

[KlausC]

I think no, just trying with flipped (p, q) in the definition of H gave nothing reasonable:

julia> H(q, p, param) = p^2/20
H (generic function with 1 method)

julia> prob = HamiltonianProblem(H, p0, q0, (0., 1.))
ODEProblem with uType ArrayPartition{Float64,Tuple{Float64,Float64}} and tType Float64. In-place: false
timespan: (0.0, 1.0)
u0: 1.00.0

julia> sol = solve(prob, dt=0.001)
retcode: Success
Interpolation: specialized 6th order interpolation
t: 7-element Array{Float64,1}:
 0.0
 0.001
 0.011
 0.07787167758977546
 0.26173425803838424
 0.5841996582033573
 1.0
u: 7-element Array{ArrayPartition{Float64,Tuple{Float64,Float64}},1}:
 (1.0, 0.0)
 (0.9999999500000004, 0.0009999999833333335)
 (0.9999939500061004, 0.010999977816680087)
 (0.9996968154128554, 0.07786380759988994)
 (0.9965767138414221, 0.2614355260679268)
 (0.9829840154332482, 0.5808823065288845)
 (0.9504152802553116, 0.9834164685913813)

julia> sol(1.0)
(0.9504152802553114, 0.9834164685913813)

julia> sol = solve(prob, SofSpa10(), dt=0.001);

julia> sol(1.0)
(0.9980287849980085, 0.03937199652442939)

[ChrisRackauckas]

Try another integrator? Maybe there's an issue there. Is prob.f what you'd expect?

[KlausC]

Is prob.f what you'd expect?

That looks fine.

julia> prob.f
(::DynamicalODEFunction{false,ODEFunction{false,DiffEqPhysics.var"#3#7"{typeof(H),Float64},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},ODEFunction{false,DiffEqPhysics.var"#5#9"{typeof(H),Float64},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing}) (generic function with 2 methods)

julia> prob.f.f1(p0, q0, nothing, 0.0)
0.0

julia> prob.f.f2(p0, q0, nothing, 0.0)
0.1

[KlausC]

Try another integrator?

That one worked as expected:

julia> sol = solve(prob, Tsit5(), dt=0.001);

julia> sol(1.0)
(1.0, 0.09999999999999959)

[KlausC]

Another one again fails:

julia> sol = solve(prob, KahanLi6(), dt=0.001);

julia> sol(1.0)
(1.0, 0.27647264980326564)

[ChrisRackauckas]

Are you sure the definition of f isn't flipped? https://github.com/SciML/DiffEqPhysics.jl/pull/58 might've been wrong.

[KlausC]

What do you mean? I will x-check #58 - there was no intention to change any logic there.

[KlausC]

No changes except (v,x,p) was renamed to (p,q,param) consistently. I will repeat my experiments with the older release...

[KlausC]

I use now the pre#58 releases

julia> VERSION
v"1.5.3"
(tmp) pkg> st
Status `/tmp/Project.toml`
  [055956cb] DiffEqPhysics v3.7.0
  [0c46a032] DifferentialEquations v6.16.0

and there is no difference in the results.

[KlausC]

When I go back more than a year, using the default algo now produces a correct result, but the other algorithms produce exactly the same false results (except of Tsit5, which remains correct).

julia> VERSION
v"1.3.1"

(tmp) pkg> st
    Status `/tmp/Project.toml`
  [055956cb] DiffEqPhysics v3.2.0
  [0c46a032] DifferentialEquations v6.0.0

[ChrisRackauckas]

I think somewhere along the lines the definition flipped the order of the equations and it's now incorrect in this library. Somehow our tests didn't catch that...

[KlausC]

I will take me the time to investigate and keep you informed.

[KlausC]

It seems, that some of the solvers in "symplectic_perform_step.jl" silently assume, that f.f2(p, q) == p for the very first calculation in perform_step! (for example lines 1133 and 1275). The corresponding initialize! functions evaluate f.f2 and leave it in a variable, which is not used by perform_step!. That is the case for KahanLi8 and SofSpa10, and also the constant variants, but is correct for SymplecticEuler. It will be necessary to check and fix each individual algorithm.