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

**URL:** <https://discourse.julialang.org/t/second-order-differential-equation-blows-up-numerically-and-i-have-no-clue-why-any-ideas/34427>\
**Category:** New to Julia\
**Created:** [February 10, 2020, 6:07pm UTC](https://discourse.julialang.org/t/second-order-differential-equation-blows-up-numerically-and-i-have-no-clue-why-any-ideas/34427 "2020-02-10T18:07:55Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![rmz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rmz/32/9861_2.png) [@rmz](https://discourse.julialang.org/u/rmz)\
**Post date:** [February 10, 2020, 6:07pm UTC](https://discourse.julialang.org/t/second-order-differential-equation-blows-up-numerically-and-i-have-no-clue-why-any-ideas/34427/1 "2020-02-10T18:07:55Z")

</div>

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.

```julia

  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

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> 

```

---

<div class="post-metadata">

**Author:** ![ericphanson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ericphanson/32/215186_2.png) [@ericphanson](https://discourse.julialang.org/u/ericphanson)\
**Post date:** [February 10, 2020, 6:57pm UTC](https://discourse.julialang.org/t/second-order-differential-equation-blows-up-numerically-and-i-have-no-clue-why-any-ideas/34427/2 "2020-02-10T18:57:57Z")

</div>

> [@rmz](#):
>
> ```julia
> x′′ = x + 2y′ + μ′*(x + μ)/D1 - μ*(x - μ′)/D2
> y′′ = y + 2x′ + μ′*y/D1 - μ*y/D2
> 
> ```

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

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

---

<div class="post-metadata">

**Author:** ![rmz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rmz/32/9861_2.png) [@rmz](https://discourse.julialang.org/u/rmz)\
**Post date:** [February 10, 2020, 7:49pm UTC](https://discourse.julialang.org/t/second-order-differential-equation-blows-up-numerically-and-i-have-no-clue-why-any-ideas/34427/3 "2020-02-10T19:49:47Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![rmz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rmz/32/9861_2.png) [@rmz](https://discourse.julialang.org/u/rmz)\
**Post date:** [February 10, 2020, 9:20pm UTC](https://discourse.julialang.org/t/second-order-differential-equation-blows-up-numerically-and-i-have-no-clue-why-any-ideas/34427/4 "2020-02-10T21:20:20Z")

</div>

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

---

<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:** [February 10, 2020, 11:43pm UTC](https://discourse.julialang.org/t/second-order-differential-equation-blows-up-numerically-and-i-have-no-clue-why-any-ideas/34427/5 "2020-02-10T23:43:03Z")

</div>

> [@ericphanson](#):
>
> 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.

Check with a smaller dt.

---

<div class="post-metadata">

**Author:** ![ericphanson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ericphanson/32/215186_2.png) [@ericphanson](https://discourse.julialang.org/u/ericphanson)\
**Post date:** [February 11, 2020, 2:54am UTC](https://discourse.julialang.org/t/second-order-differential-equation-blows-up-numerically-and-i-have-no-clue-why-any-ideas/34427/6 "2020-02-11T02:54:11Z")

</div>

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](https://scicomp.stackexchange.com/a/29154) 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).

---

<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:** [February 11, 2020, 5:57am UTC](https://discourse.julialang.org/t/second-order-differential-equation-blows-up-numerically-and-i-have-no-clue-why-any-ideas/34427/7 "2020-02-11T05:57:44Z")

</div>

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](https://global.discourse-cdn.com/julialang/original/3X/3/e/3e5713561cc7dd9fe9a9efeb02044f7d8ed92e5b.png)

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.

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/1/3/1365474f4fa046e1aa2d455fd4d9c12df74a49ce.png)

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

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/b/d/bdd5b47f40958f9df7a5cfd7efecdf4e4e8802e0.png)

But we’re good with a highly accurate solution:

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/7/9/79fafb49a9cf51a2e233d68502b1f23242ca90d1.png)

However, by 5 periods we get:

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/8/e/8eb1b3386d6f0de012fb678194f71aed8b2fc6c9.png)

so then in the 6th:

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/4/0/40ca56e115335a63d7e44a17f94883fa4570aaa9.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 [DPRKNs diverge on an example · Issue #1030 · SciML/OrdinaryDiffEq.jl · GitHub](https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/issues/1030)

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [February 11, 2020, 8:03am UTC](https://discourse.julialang.org/t/second-order-differential-equation-blows-up-numerically-and-i-have-no-clue-why-any-ideas/34427/8 "2020-02-11T08:03:45Z")

</div>

> [@ChrisRackauckas](#):
>
> not even symplectic properties can save you from chaos

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

---

<div class="post-metadata">

**Author:** ![ericphanson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ericphanson/32/215186_2.png) [@ericphanson](https://discourse.julialang.org/u/ericphanson)\
**Post date:** [February 11, 2020, 12:02pm UTC](https://discourse.julialang.org/t/second-order-differential-equation-blows-up-numerically-and-i-have-no-clue-why-any-ideas/34427/9 "2020-02-11T12:02:07Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ericphanson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ericphanson/32/215186_2.png) [@ericphanson](https://discourse.julialang.org/u/ericphanson)\
**Post date:** [February 11, 2020, 12:15pm UTC](https://discourse.julialang.org/t/second-order-differential-equation-blows-up-numerically-and-i-have-no-clue-why-any-ideas/34427/10 "2020-02-11T12:15:27Z")

</div>

Thanks very much for the investigation 🙂

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [February 11, 2020, 12:21pm UTC](https://discourse.julialang.org/t/second-order-differential-equation-blows-up-numerically-and-i-have-no-clue-why-any-ideas/34427/11 "2020-02-11T12:21:20Z")

</div>

thanks, I missed that one.

---

<div class="post-metadata">

**Author:** ![rmz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rmz/32/9861_2.png) [@rmz](https://discourse.julialang.org/u/rmz)\
**Post date:** [February 11, 2020, 12:22pm UTC](https://discourse.julialang.org/t/second-order-differential-equation-blows-up-numerically-and-i-have-no-clue-why-any-ideas/34427/12 "2020-02-11T12:22:35Z")

</div>

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