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

New to Julia
1

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>
2

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
3

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

3 Likes
4

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

3 Likes
5

Check with a smaller dt.

6

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).

7

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:

Capture

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]])

plot

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")

plot

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")

plot

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")

plot

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")

plot

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 DPRKNs diverge on an example · Issue #1030 · SciML/OrdinaryDiffEq.jl · GitHub

10 Likes
8

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

9

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
10

Thanks very much for the investigation :slight_smile:

1 Like
11

thanks, I missed that one.

1 Like
12

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

2 Likes