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

rmz · February 10, 2020, 6:07pm

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>

ericphanson · February 10, 2020, 6:57pm

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.

rmz · February 10, 2020, 7:49pm

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

rmz · February 10, 2020, 9:20pm

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

ChrisRackauckas · February 10, 2020, 11:43pm

Check with a smaller dt.

ericphanson · February 11, 2020, 2:54am

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

ChrisRackauckas · February 11, 2020, 5:57am

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

Tamas_Papp · February 11, 2020, 8:03am

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

ericphanson · February 11, 2020, 12:02pm

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.

ericphanson · February 11, 2020, 12:15pm

Thanks very much for the investigation

Tamas_Papp · February 11, 2020, 12:21pm

thanks, I missed that one.

rmz · February 11, 2020, 12:22pm

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

Topic		Replies	Views
My second order ODE explodes Modelling & Simulations	1	75	August 29, 2024
Numerical instability of partial differential equations General Usage diffeq	6	964	May 17, 2020
Vern9 from differentialequations.jl and ode89 from matlab Modelling & Simulations	15	453	May 25, 2024
Instability detected in solving ODE with DifferentialEquations.jl Numerics differentialequation	7	405	June 9, 2023
Solving ODEs arising from spherical and cylindrical systems? General Usage diffeq	8	1744	April 25, 2019

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

Related topics