SecondOrderODEProblem gives wrong result

[IBUzPE9]

I solved the differential equation in two ways. Should the solution2 be the same as the solution1?

using DifferentialEquations

println("start");
g = 9.81; # gravity acceleration, m/s2
Kd = 0.0006;
V0x = 470.0;
V0y = 170.0;
Tspan = (0.0, 50.0)

# move equations 2nd order 
# x''(t) == -Kd * x'(t) * sqrt( x'(t)^2 + y'(t)^2 )
# y''(t) == -Kd * y'(t) * sqrt( x'(t)^2 + y'(t)^2 ) - g
# x'(0) == V0x
# y'(0) == V0y
# x(0) == 0
# y(0) == 0

function move_eq2(ddu,du,u,p,t)
    local vmod = sqrt(du[1]^2 + du[2]^2);

    ddu[1] = -Kd * du[1] * vmod; # x''
    ddu[2] = -Kd * du[2] * vmod - g; # y''
end 

# move equations 1st order 
# vx'(t) == -Kd * vx(t) * sqrt( vx(t)^2 + vy(t)^2 )
# vy'(t) == -Kd * vy(t) * sqrt( vx(t)^2 + vy(t)^2 ) - g
# x'(t) == vx(t)
# y'(t) == vy(t)
# vx(0) == V0x
# vy(0) == V0y
# x(0) == 0
# y(0) == 0

function move_eq1(du,u,p,t)
    local vmod = sqrt(u[1]^2 + u[2]^2);

    du[1] = -Kd * u[1] * vmod; # vx'
    du[2] = -Kd * u[2] * vmod - g; # vy'
    du[3] = u[1]; # x'
    du[4] = u[2]; # y'
end 

prob1 = ODEProblem(move_eq1,[V0x, V0y, 0.0, 0.0], Tspan);
prob2 = SecondOrderODEProblem(move_eq2,[V0x, V0y], [0.0, 0.0], Tspan);

solution1 = solve(prob1);
solution2 = solve(prob2);

using Plots
#plotly();
gr();

lab = ["vx" "vy" "x" "y"];
plot(
    plot(solution1, title="ODEProblem", label=lab), 
    plot(solution2, title="SecondOrderODEProblem", label=lab), 
    layout = (2, 1), legend=:outertopright, size=(1920,1080))
savefig("soode-bug.png");
println("done");

[ChrisRackauckas]

What if you reduce the tolerances?

[IBUzPE9]

Something like this?

solution2 = solve(prob2, abstol=1e-10, reltol=1e-10);

Looks better, but not the same. The x still far from ODEProblem result.

[IBUzPE9]

With abstol=1e-20, reltol=1e-20 the solutions are almost the same. Is this normal? With the default tolerances the SecondOrderODEProblem solution looks quite strange.

[ChrisRackauckas]

Which algorithms is it comparison? sol.alg. Indeed it looks a little strange. But the issue may come down to the fact that on a linear equation the error estimators can be more loose. It's odd that it's that loose though...

[HenryLangner]

DPRKN6 was used by default here, however it is designed for 2nd order ODEs which are not dependent on the first derivative. This could have lead to the wrong results. More information can be found here:

https://github.com/SciML/OrdinaryDiffEq.jl/issues/1938

The default solver for SecondOrderODEProblems has been changed since so this should not be an issue anymore.