 # Second Order DE producing incorrect results

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^2 + u^2)^1.5

du = u
du = u
du = u*(f + 3) + 2u
du = u*f - 2u

end

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

u0 = zeros(4)

for i in 1:index

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

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^2 + u^2)^1.5

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

end

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

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

for i in 1:index

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

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

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^2 + u^2)^1.5

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

end

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

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

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

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^2 + u^2)^1.5

du = u
du = u
du = u*(f + 3) + 2u
du = u*f - 2u

end

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

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

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
1 Like

Thank you!

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