# Second order differential equation blows up numerically, and I have no clue why. Any ideas?

Hi

I’m still in the newbie category so I dare ask this question as a “First Steps” question.

I read this blog:

`````` https://www.johndcook.com/blog/2020/02/08/arenstorf-orbit/
``````

… about an equation that was used to construct orbits for the Apollo program. It looked interesting, and not too complex, so I thought I could implement that in Julia. So I tried, and failed. The program is “numerically unstable” Julia tells me.

Dear Julia community, can you at a glance see what I’m doing wrong? I’ve tried the DPRKN12, DPRKN6 and KahanLi8 solvers, but they seem to be failing in more or less the same way so it’s probably something other than the solver that is wrong, but I have no idea what.

``````
using DifferentialEquations
using OrdinaryDiffEq
using Plots

function plot_arensdorf()
μ   = 0.012277471
μ′  = 1 - μ

function arensdorf_orbit(du, u, z, p, t)
x, y  = z
x′,y′ = u

D(z) = ((x + z )^2  + y^2)^(3/2)
D1 = D(μ)
D2 = D(-μ′)

x′′ = x + 2y′  + μ′*(x + μ)/D1 - μ*(x - μ′)/D2
y′′ = y + 2x′  + μ′*y/D1       - μ*y/D2

du[1] = x′′
du[2] = y′′
end

initial_positions =  [0.994,  0]
initial_velocities = [0.0,   -2.001585106]
tspan = (0.0, 5000.0)
prob = SecondOrderODEProblem(arensdorf_orbit, initial_velocities, initial_positions, tspan)
sol = solve(prob, DPRKN12(), dt=1/100);
plot(sol,vars=(1,2))
end
``````
``````
julia> plot_arensdorf()
┌ Warning: Instability detected. Aborting
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/YIwj5/src/integrator_interface.jl:339
Error showing value of type Plots.Plot{Plots.GRBackend}:
ERROR: InexactError: trunc(Int64, -Inf)
Stacktrace:
[1] trunc at ./float.jl:702 [inlined]
[2] ceil at ./float.jl:365 [inlined]
[3] optimize_ticks_typed(::Float64, ::Float64, ::Bool, ::Array{Tuple{Float64,Float64},1}, ::Int64, ::Int64, ::Int64, ::Float64, ::Float64, ::Float64, ::Float64, ::Bool, ::Nothing) at /Users/rmz/.julia/packages/PlotUtils/6AioP/src/ticks.jl:183
[4] #optimize_ticks#36 at /Users/rmz/.julia/packages/PlotUtils/6AioP/src/ticks.jl:137 [inlined]
[5] (::PlotUtils.var"#kw##optimize_ticks")(::NamedTuple{(:k_min, :k_max),Tuple{Int64,Int64}}, ::typeof(optimize_ticks), ::Float64, ::Float64) at ./none:0
[6] optimal_ticks_and_labels(::Plots.Subplot{Plots.GRBackend}, ::Plots.Axis, ::Nothing) at /Users/rmz/.julia/packages/Plots/qZHsp/src/axes.jl:188
[7] optimal_ticks_and_labels at /Users/rmz/.julia/packages/Plots/qZHsp/src/axes.jl:156 [inlined]
[8] get_ticks(::Plots.Subplot{Plots.GRBackend}, ::Plots.Axis) at /Users/rmz/.julia/packages/Plots/qZHsp/src/axes.jl:266
[9] _update_min_padding!(::Plots.Subplot{Plots.GRBackend}) at /Users/rmz/.julia/packages/Plots/qZHsp/src/backends/gr.jl:902
[10] iterate at ./generator.jl:47 [inlined]
[11] _collect(::Array{AbstractLayout,2}, ::Base.Generator{Array{AbstractLayout,2},typeof(Plots._update_min_padding!)}, ::Base.EltypeUnknown, ::Base.HasShape{2}) at ./array.jl:635
[12] collect_similar(::Array{AbstractLayout,2}, ::Base.Generator{Array{AbstractLayout,2},typeof(Plots._update_min_padding!)}) at ./array.jl:564
[13] map(::Function, ::Array{AbstractLayout,2}) at ./abstractarray.jl:2073
[14] _update_min_padding!(::Plots.GridLayout) at /Users/rmz/.julia/packages/Plots/qZHsp/src/layouts.jl:310
[15] prepare_output(::Plots.Plot{Plots.GRBackend}) at /Users/rmz/.julia/packages/Plots/qZHsp/src/plot.jl:261
[16] display(::Plots.PlotsDisplay, ::Plots.Plot{Plots.GRBackend}) at /Users/rmz/.julia/packages/Plots/qZHsp/src/output.jl:143
[17] display(::Any) at ./multimedia.jl:323
[18] #invokelatest#1 at ./essentials.jl:709 [inlined]
[19] invokelatest at ./essentials.jl:708 [inlined]
[20] print_response(::IO, ::Any, ::Bool, ::Bool, ::Any) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/REPL/src/REPL.jl:156
[21] print_response(::REPL.AbstractREPL, ::Any, ::Bool, ::Bool) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/REPL/src/REPL.jl:141
[22] (::REPL.var"#do_respond#38"{Bool,REPL.var"#48#57"{REPL.LineEditREPL,REPL.REPLHistoryProvider},REPL.LineEditREPL,REPL.LineEdit.Prompt})(::Any, ::Any, ::Any) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/REPL/src/REPL.jl:719
[23] #invokelatest#1 at ./essentials.jl:709 [inlined]
[24] invokelatest at ./essentials.jl:708 [inlined]
[25] run_interface(::REPL.Terminals.TextTerminal, ::REPL.LineEdit.ModalInterface, ::REPL.LineEdit.MIState) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/REPL/src/LineEdit.jl:2306
[26] run_frontend(::REPL.LineEditREPL, ::REPL.REPLBackendRef) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/REPL/src/REPL.jl:1045
[27] run_repl(::REPL.AbstractREPL, ::Any) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/REPL/src/REPL.jl:201
[28] (::Base.var"#770#772"{Bool,Bool,Bool,Bool})(::Module) at ./client.jl:382
[29] #invokelatest#1 at ./essentials.jl:709 [inlined]
[30] invokelatest at ./essentials.jl:708 [inlined]
[31] run_main_repl(::Bool, ::Bool, ::Bool, ::Bool, ::Bool) at ./client.jl:366
[32] exec_options(::Base.JLOptions) at ./client.jl:304
[33] _start() at ./client.jl:460

julia>

``````

Looking at the link, it seems like some of the signs are off here; it should be

``````x′′ = x + 2y′  - μ′*(x + μ)/D1 - μ*(x - μ′)/D2
y′′ = y - 2x′  - μ′*y/D1       - μ*y/D2
``````

With that fix, I don’t get instabilities warnings with `KahanLi8()`, but it looks like one also wants much smaller `tspan` to get something similar to that in the blog post (e.g. `tspan = (0.0, 10.0)`). Going up to 5000, it spirals out massively; not sure if that’s numerical error adding up or if it’s supposed to do that.

4 Likes

Thank you so much . The typo is embarrassing, but your comments were uplifting

3 Likes

I put it into a notebook . https://github.com/la3lma/julia-playground/blob/master/Arensdorf.ipynb

3 Likes

Check with a smaller dt.

Thanks for the hint-- at `dt=1/1000` the magnitude of the drift went down a lot (but is still ~`1e75`), and at `dt=1/2000` it’s still running (with `tspan = (0.0,5000.0)`), but I assume it will be even better. I don’t quite understand though: I thought symplectic integrators were supposed to stay within some bounded distance from the real trajectory for all time, since the discretized system also lives on some nearby symplectic manifold? I’ve read your great stack exchange post and thought I understood their limitations, but it seemed like this is a perfect application (very long time behavior with a desire for moderate accuracy).

Okay, you made me look into it. You really did. Turns out, the reason why Hairer investigated it was because it’s a good toy problem to make integrators go bad. If you look at what Hairer plots you’ll see why:

He was showcasing how quickly things can drift. Indeed, with a simple Tsit5 we get a good solution. Note how I’m plotting in position space.

``````using OrdinaryDiffEq
using Plots

const μ   = 0.012277471
const μ′  = 1 - μ

function arensdorf_orbit(du, u, z, p, t)
x, y  = z
x′,y′ = u

D(z) = ((x + z )^2  + y^2)^(3/2)
D1 = D(μ)
D2 = D(-μ′)

x′′ = x + 2y′  - μ′*(x + μ)/D1 - μ*(x - μ′)/D2
y′′ = y - 2x′  - μ′*y/D1       - μ*y/D2

du[1] = x′′
du[2] = y′′
end

initial_positions =  [0.994,  0]
initial_velocities = [0.0,   -2.00158510637908252240537862224]
N = 1
tend = 17.0652165601579625588917206249
tspan = (0.0, N*tend)
prob = SecondOrderODEProblem(arensdorf_orbit, initial_velocities, initial_positions, tspan)
sol = solve(prob, Tsit5())
plot(sol,vars=(3,4))
scatter!([sol[3,1]],[sol[4,1]])
scatter!([sol[3,end]],[sol[4,end]])
``````

Bingo, it’s periodic… right? No, it has slightly drifted. And oh boy, a slight drift really throws it off for the second period:

``````N = 2
tend = 17.0652165601579625588917206249
tspan = (0.0, N*tend)
prob = SecondOrderODEProblem(arensdorf_orbit, initial_velocities, initial_positions, tspan)
sol = solve(prob, Tsit5())
plot(sol,vars=(3,4))
scatter!([sol[3,1]],[sol[4,1]])
scatter!([sol[3,end]],[sol[4,end]])
savefig("plot.png")
``````

But we’re good with a highly accurate solution:

``````N = 2
tend = 17.0652165601579625588917206249
tspan = (0.0, N*tend)
prob = SecondOrderODEProblem(arensdorf_orbit, initial_velocities, initial_positions, tspan)
sol = solve(prob, Vern9(), abstol=1e-12, reltol=1e-12)
plot(sol,vars=(3,4))
scatter!([sol[3,1]],[sol[4,1]])
scatter!([sol[3,end]],[sol[4,end]])

savefig("plot.png")
``````

However, by 5 periods we get:

``````N = 5
tend = 17.0652165601579625588917206249
tspan = (0.0, N*tend)
prob = SecondOrderODEProblem(arensdorf_orbit, initial_velocities, initial_positions, tspan)
sol = solve(prob, Vern9(), abstol=1e-12, reltol=1e-12)
plot(sol,vars=(3,4))
scatter!([sol[3,1]],[sol[4,1]])
scatter!([sol[3,end]],[sol[4,end]])

savefig("plot.png")
``````

so then in the 6th:

``````N = 6
tend = 17.0652165601579625588917206249
tspan = (0.0, N*tend)
prob = SecondOrderODEProblem(arensdorf_orbit, initial_velocities, initial_positions, tspan)
sol = solve(prob, Vern9(), abstol=1e-12, reltol=1e-12)
plot(sol,vars=(3,4))
scatter!([sol[3,1]],[sol[4,1]])
scatter!([sol[3,end]],[sol[4,end]])

savefig("plot.png")
``````

And boom, notice how when it’s slightly different it can end up exponentially diverging in the trajectory. This is a sign that the long term behavior is just caused by sensitive dependence to initial conditions, and not even symplectic properties can save you from chaos.

Though it looks like DPRKN6/DPRKN12 have a bug https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/issues/1030

10 Likes

Sorry if I missed something, but are any of the integrators in the code above symplectic?

The `KahanLi8()` I was using earlier in the thread is symplectic, and I think Chris was responding to my confusion about why that wasn’t enough to give good long time behavior.

1 Like

Thanks very much for the investigation

1 Like

thanks, I missed that one.

1 Like

I’ll second that. I’ve had a smile on my face all day because of this Thank you so much.

2 Likes