# Second Order DE producing incorrect results

**URL:** <https://discourse.julialang.org/t/second-order-de-producing-incorrect-results/59380>\
**Category:** Numerics\
**Tags:** diffeq, ode\
**Created:** [April 16, 2021, 1:52am UTC](https://discourse.julialang.org/t/second-order-de-producing-incorrect-results/59380 "2021-04-16T01:52:31Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![daviiinsky](https://avatars.discourse-cdn.com/v4/letter/d/db5fbb/32.png) [@daviiinsky](https://discourse.julialang.org/u/daviiinsky)\
**Post date:** [April 16, 2021, 1:52am UTC](https://discourse.julialang.org/t/second-order-de-producing-incorrect-results/59380/1 "2021-04-16T01:52:31Z")

</div>

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](https://en.wikipedia.org/wiki/Hill_differential_equation). 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!

```nohighlight
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()

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 16, 2021, 9:17am UTC](https://discourse.julialang.org/t/second-order-de-producing-incorrect-results/59380/2 "2021-04-16T09:17:01Z")

</div>

If you simplify the problem:

```julia
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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 16, 2021, 9:33am UTC](https://discourse.julialang.org/t/second-order-de-producing-incorrect-results/59380/3 "2021-04-16T09:33:19Z")

</div>

I made an issue: [DPRKN6 gives substantially different answer on a second order ODE than Tsit5 · Issue #1372 · SciML/OrdinaryDiffEq.jl · GitHub](https://github.com/SciML/OrdinaryDiffEq.jl/issues/1372)

---

<div class="post-metadata">

**Author:** ![daviiinsky](https://avatars.discourse-cdn.com/v4/letter/d/db5fbb/32.png) [@daviiinsky](https://discourse.julialang.org/u/daviiinsky)\
**Post date:** [April 16, 2021, 10:38pm UTC](https://discourse.julialang.org/t/second-order-de-producing-incorrect-results/59380/4 "2021-04-16T22:38:10Z")

</div>

Thank you!

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