Second Order DE producing incorrect results

Specific Domains Numerics
,
1

Hello!

I was working on the following problem where I’m integrating a second order differential equation. Initially I tried putting it into the 2nd Order Solver in DifferentialEquations.jl, where it did produce a plot set, but when comparing to colleagues they were not correct.

Later, I changed the 2nd Order ODE System to a First Order ODE System, which this one does produce the correct results. I’m able to work with the First Order System, but now I’m more curious as to why my implementation of the 2nd Order version is producing different results. For context, these are the Hill equations. I tried different integrators for the solve() , and some did change the plots, but none reproduced the plots the first order system created.

I provided the code block beneath that should show the two different cases.

Any understanding on why the results are different would be appreciated!

using DifferentialEquations
using Plots

function main()

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

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

    # system_ode(b, t0, tspan)
    second_ode(b, t0, tspan)

end

function system_ode(b, t0, tspan)

    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)

    for i in 1:index

        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)
        orbits[i] = plot(sol[1,:], sol[2,:], label="b = "*string(b[i]), framestyle=:origin, xlim=(-2,2), ylim=(-2,2))

    end

    plot(orbits..., layout=index)
end

function second_ode(b, t0, tspan)

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

        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)

    for i in 1:index

        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)
        orbits[i] = plot(sol[1,:], sol[2,:], label="b = "*string(b[i]), framestyle=:origin, xlim=(-2,2), ylim=(-2,2))

    end

    plot(orbits..., layout=index)
end


main()
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
3

I made an issue: DPRKN6 gives substantially different answer on a second order ODE than Tsit5 · Issue #1372 · SciML/OrdinaryDiffEq.jl · GitHub

1 Like
4

Thank you!

And I should’ve checked the column output more diligently. I was racking my head over why the output looked off.