Second Order DE producing incorrect results

Specific Domains Numerics
,
2

If you simplify the problem:

using DifferentialEquations
using Plots

b = range(1, 3, length=6)

t0 = -100.
t1 = 150.
tspan = (t0, t1)
i = 1

function hill!(ddu, du, u, p, t)

    @show u

    f = -3/(u[1]^2 + u[2]^2)^1.5

    ddu[1] = u[1]*(f + 3) + 2du[2]
    ddu[2] = u[2]*f - 2du[1]

end

index = size(b, 1)
orbits = Array{Any}(nothing, index)

u0 = zeros(2)
du0 = zeros(2)

u0[1] = b[i] - 8/(3b[i]^2 * (t0^2 + 4/9)^0.5)
u0[2] = -1.5*u0[1]*t0
du0[1] = 0
du0[2] = -1.5*u0[1]

prob = SecondOrderODEProblem(hill!, du0, u0, tspan, abstol=1e-8, reltol=1e-8)

sol = solve(prob, Tsit5())
plot(sol)

function hill!(du, u, p, t)

    f = -3/(u[1]^2 + u[2]^2)^1.5

    du[1] = u[3]
    du[2] = u[4]
    du[3] = u[1]*(f + 3) + 2u[4]
    du[4] = u[2]*f - 2u[3]

end

index = size(b, 1)
orbits = Array{Any}(nothing, index)

u0 = zeros(4)
u0[1] = b[i] - 8/(3b[i]^2 * (t0^2 + 4/9)^0.5)
u0[2] = -1.5*u0[1]*t0
u0[3] = 0
u0[4] = -1.5*u0[1]

prob = ODEProblem(hill!, u0, tspan, abstol=1e-8, reltol=1e-8)

sol = solve(prob, Tsit5())
plot(sol)

You’ll see you plotted the wrong variables and the standard solvers work fine. Something is up with the default DPRKN6 that is being chosen there: something seems flipped. Can you open an issue?

2 Likes