# Homoclinic loop accuracy DifferentialEquations.jl

**URL:** <https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456>\
**Category:** Modelling & Simulations\
**Created:** [April 14, 2023, 7:10am UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456 "2023-04-14T07:10:59Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Sergey\_Novak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sergey_novak/32/37716_2.png) [@Sergey\_Novak](https://discourse.julialang.org/u/Sergey_Novak)\
**Post date:** [April 14, 2023, 7:10am UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/1 "2023-04-14T07:10:59Z")

</div>

Hello!  
I try get homoclinic loop to saddle focus for system of Rossler. With the given parameters, the homoclinic definitely exists, according to some articles. To get homoclinic loop, I integrate in backward time. When I do this I get an instability message. Why I get such a message I understand. I suppose the problem is with accuracy. In the same articles, people managed to get a homoclinic loop using the 4-order Runge Kutta. I tried to use Verner9, Feigin12. I used big float and it didn’t help. I also made the time step quite small. What could be the problem? Is it all about accuracy, or am I doing something wrong even before the integration?

Code:

```julia
using StaticArrays, DifferentialEquations, DynamicalSystems

x, y, z = -10..10, -10..10, -10..10
box = x × y × z

using ForwardDiff
using LinearAlgebra

@inbounds function Res(u, p, t)
    du1 = -u[2]-u[3]
    du2 = u[1]+p[1]*u[2]
    du3 = p[2]*u[1]-p[3]*u[3]+u[1]*u[3]
    return SVector(du1, du2, du3)
end
@inbounds function jac_res(u, p, t)
    SMatrix{3, 3}(0.0, 1.0, p[2]+u[3],
                -1.0, p[1], 0.0,
                -1.0, 0.0, -p[3]+u[1])
end

time_set = (t = 500, Ttr = 250, tstep = 0.001)
integ_set = (alg = Vern9(), adaptive = false, dt = time_set.tstep)

u0 = SA[0.13060105886751308, 0.4782615699018184, 0.7435866656901606]
b = 0.3; c = 4.89; a = 0.35;
p = [a, b, c]
ds = CoupledODEs(Res, u0, p, diffeq = integ_set)

fp, ei, _ = fixedpoints(ds, box, jac_res)

Jacobian = jac_res(fp[1], p, 0)
eivecs = eigvecs(Jacobian)
vec_stable = real(eivecs[:, 1])
ϵ = 1e-23
shift = SA[0.0, 0.0, 0.0] + vec_stable*ϵ

Jacobi(u0, p, t) = ForwardDiff.jacobian((x) -> Res(x, p, 0), u0)

prob = ODEProblem(Res, shift, (0.0, -1000.0),p)

sol = solve(prob,alg = RK4(), adaptive = false, dt = 0.001);

```

---

<div class="post-metadata">

**Author:** ![cmarcotte](https://avatars.discourse-cdn.com/v4/letter/c/a3d4f5/32.png) [@cmarcotte](https://discourse.julialang.org/u/cmarcotte)\
**Post date:** [April 14, 2023, 8:52am UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/2 "2023-04-14T08:52:44Z")

</div>

A couple things:  
You didn’t define your `a,b,c` so we can’t run your code as you have (though I just used `a=0.1,b=0.1,c=14.0`).

Additionally, if you plot [my] solution, the divergence begins in earnest around `t=-16.0`, just before which the z-component becomes negative; this is consistent for adaptive stepping and fixed (though the precise time changes). This seems like a problem with solving the system in backward time – \dot{z} \< 0 if z=0 and p\_2 u\_1 \> 0, which is eminently possible. In fact, using `Roots.jl`, we can find the precise times that z=0, `t0 = find_roots((t)->sol(t)[3], tspan)`, and then check whether the (negative) flow has a negative component along z.

 ![Solution of Backward-Time Rossler from perturbed equilibrium](https://global.discourse-cdn.com/julialang/original/3X/8/a/8ac543e7f77d2b0ce5fc0cf7f924b9fa3fbf9c23.png)

---

<div class="post-metadata">

**Author:** ![Sergey\_Novak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sergey_novak/32/37716_2.png) [@Sergey\_Novak](https://discourse.julialang.org/u/Sergey_Novak)\
**Post date:** [April 14, 2023, 9:03am UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/3 "2023-04-14T09:03:47Z")

</div>

Sorry, I added parameter values  
I don’t quite understand what you mean  
I think the problem is the accuracy. If I did everything correctly, then the starting point corresponds to the displacement along the stable manifold from the saddle focus.

---

<div class="post-metadata">

**Author:** ![empet](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/empet/32/221303_2.png) [@empet](https://discourse.julialang.org/u/empet)\
**Post date:** [April 14, 2023, 9:34am UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/4 "2023-04-14T09:34:12Z")

</div>

I recommend reading this recent paper by @ranocha: [https://arxiv.org/abs/2304.02365](https://arxiv.org/abs/2304.02365), to understand why you are surprised by the computed homoclinic orbit.

---

<div class="post-metadata">

**Author:** ![empet](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/empet/32/221303_2.png) [@empet](https://discourse.julialang.org/u/empet)\
**Post date:** [April 14, 2023, 9:56am UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/5 "2023-04-14T09:56:30Z")

</div>

Running your code I noticed that you computed the eigenvalues and eigenvectors of Jacobian matrix evaluated at u0, but u0 is not an equilibrium point for your system. Or you are searching for a homoclinic orbit to a particular equilibrium.  
To check it I defined the vector field, having as components the RHS of your system, ie.:

```julia
F(u::SVector{3,Float64}; p = [0.35, 0.3, 4.89])=(-u[2]-u[3], u[1]+p[1]*u[2], p[2]*u[1]-p[3]*u[3]+u[1]*u[3])

```

F evaluated at u0 gives:

```julia
F(u0)
(-1.221848235591979, 0.2979926083331495, -3.4998452716657327)

```

not (0,0,0).

---

<div class="post-metadata">

**Author:** ![Sergey\_Novak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sergey_novak/32/37716_2.png) [@Sergey\_Novak](https://discourse.julialang.org/u/Sergey_Novak)\
**Post date:** [April 14, 2023, 10:30am UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/6 "2023-04-14T10:30:37Z")

</div>

Why `u0`?  
Eigenvalues and eigenvectors calculate for `fp[1]`  
Fragment code:

```julia
fp, ei, _ = fixedpoints(ds, box, jac_res)

Jacobian = jac_res(fp[1], p, 0)
eivecs = eigvecs(Jacobian)
vec_stable = real(eivecs[:, 1])
ϵ = 1e-23
shift = SA[0.0, 0.0, 0.0] + vec_stable*ϵ

```

Thank you for link on article

---

<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:** [April 14, 2023, 10:42am UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/7 "2023-04-14T10:42:56Z")

</div>

> [@Sergey\_Novak](#):
>
> I suppose the problem is with accuracy.

Well I mean, you turned off the accuracy handling. Don’t use `adaptive=false` if you care about solving to tolerance or high accuracy. Set the `abstol` and `reltol` to the desired local tolerance and let the integrator solve to tolerance.

---

<div class="post-metadata">

**Author:** ![Sergey\_Novak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sergey_novak/32/37716_2.png) [@Sergey\_Novak](https://discourse.julialang.org/u/Sergey_Novak)\
**Post date:** [April 14, 2023, 10:58am UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/8 "2023-04-14T10:58:28Z")

</div>

With next tol i get message about unstable solution

```julia
sol = solve(prob,alg = RK4(), abstol=1e-13, reltol=1e-13, maxiters = 10000000)

```

```julia
dt(-1.1368683772161603e-13) <= dtmin(1.1368683772161603e-13) at t=-13.637725296644767. Aborting. There is either an error in your model specification or the true solution is unstable.

```

---

<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:** [April 14, 2023, 11:03am UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/9 "2023-04-14T11:03:14Z")

</div>

That means RK4 cannot hit that accuracy within the chosen dt, which isn’t too surprising because it’s low order. What about Vern9?

---

<div class="post-metadata">

**Author:** ![Sergey\_Novak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sergey_novak/32/37716_2.png) [@Sergey\_Novak](https://discourse.julialang.org/u/Sergey_Novak)\
**Post date:** [April 14, 2023, 11:05am UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/10 "2023-04-14T11:05:09Z")

</div>

I tried the Vern9. It also displays a message about an unstable solution

---

<div class="post-metadata">

**Author:** ![cmarcotte](https://avatars.discourse-cdn.com/v4/letter/c/a3d4f5/32.png) [@cmarcotte](https://discourse.julialang.org/u/cmarcotte)\
**Post date:** [April 14, 2023, 11:08am UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/11 "2023-04-14T11:08:33Z")

</div>

One issue is it’s not vanilla Rössler:  
\dot{x} = -y-z, \dot{y} = x + a y, and \dot{z} = b + z(x-c), but instead  
\dot{x} = -y -z, \dot{y} = x + a y, and \dot{z} = b x + z(x-c).

Changing it to Rossler original edition, you immediately blow up (even with `Feagin14()` and `abstol=reltol=1e-13`. The issue is that with your initial condition at (approximately) z=0 immediately encounters \dot{z} \< 0 (in backward time), and after that, \dot{z} \propto z, so the solution blows up. I suspect the same happens with your system – it steps over z=0 because the p\_2u\_1 term is negative at that time for the backwards flow. If you use a `ContinuousCallback` to ensure that z\geq 0, then perhaps it will be more reliable. I’m not convinced this isn’t a mathematical issue, at heart, rather than a numerical one. You’re essentially hoping that the attractor is stable and bounded under the reverse flow, but thats not obvious to me.

So my recommendations are:

1. change the flow to Rössler, if that’s what you actually want,
2. Try a different initial condition with low tolerances and high order and adaptivity to see if the flow is actually bounded in backward time,
3. implement a `ContinuousCallback` to keep z \geq 0 (or z \> 0, probably).

---

<div class="post-metadata">

**Author:** ![Sergey\_Novak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sergey_novak/32/37716_2.png) [@Sergey\_Novak](https://discourse.julialang.org/u/Sergey_Novak)\
**Post date:** [April 14, 2023, 11:24am UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/12 "2023-04-14T11:24:00Z")

</div>

The system for which I actually want to get a loop is different. I took the Rössler system to get a loop in a simpler system. How i can correct use Callback? More precisely, what should he do after checking for the condition?

---

<div class="post-metadata">

**Author:** ![empet](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/empet/32/221303_2.png) [@empet](https://discourse.julialang.org/u/empet)\
**Post date:** [April 14, 2023, 11:33am UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/13 "2023-04-14T11:33:22Z")

</div>

You’re right, I don’t know how I read in a hurry on the REPL 🙂

---

<div class="post-metadata">

**Author:** ![cmarcotte](https://avatars.discourse-cdn.com/v4/letter/c/a3d4f5/32.png) [@cmarcotte](https://discourse.julialang.org/u/cmarcotte)\
**Post date:** [April 14, 2023, 11:43am UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/14 "2023-04-14T11:43:11Z")

</div>

For the callback I can’t tell you what _should_ be done – you have the ability to modify the state values to prevent z from becoming negative, but what particular form that takes is up to you and how you wish to perturb\* your system.

\*state modification is not guaranteed to be ‘perturbative’ in the mathematical sense.

---

<div class="post-metadata">

**Author:** ![Sergey\_Novak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sergey_novak/32/37716_2.png) [@Sergey\_Novak](https://discourse.julialang.org/u/Sergey_Novak)\
**Post date:** [April 28, 2023, 12:30pm UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/15 "2023-04-28T12:30:48Z")

</div>

I am experimenting with the Feagin 12. Is it possible to use type of data float64 for this method?  
With next code i don’t get message with warnings about maxiter.

```julia
prob_for = ODEProblem(TM, shift, (0.0, 500), p)
sol_for = solve(prob_for, alg = Feagin12(), abstol = 1e-18, reltol = 1e-18);

```

I tried use Feagin14 and i get message about Instability, when using abstol = 1e-16, reltol = 1e-16

---

<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:** [April 30, 2023, 2:23pm UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/16 "2023-04-30T14:23:14Z")

</div>

> [@Sergey\_Novak](#):
>
> Is it possible to use type of data float64 for this method?

Yes, as always it’s just done by how you choose your initial values.

> [@Sergey\_Novak](#):
>
> `reltol = 1e-18`

```julia
julia> eps(Float64)
2.220446049250313e-16

```

so a relative tolerance that low is generally impossible with Float64. Anything in Float64 range I would recommend Vern9 anyways.

---

<div class="post-metadata">

**Author:** ![Sergey\_Novak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sergey_novak/32/37716_2.png) [@Sergey\_Novak](https://discourse.julialang.org/u/Sergey_Novak)\
**Post date:** [July 18, 2023, 10:03am UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/17 "2023-07-18T10:03:46Z")

</div>

I decided to try implicit methods to solve stiff problems. With them, the maxiters option is not required, and there is no need to calculate the trajectory for a long time to get an instability message.

The radau method does not work and the following error appears. Did I forget any options for this method?

```julia
sol = solve(prob, ODEInterfaceDiffEq.radau())

```

```julia
MethodError: no method matching radau(::ODEInterfaceDiffEq.var"#12#18"{ODEProblem{SVector{3, Float64}, Tuple{Float64, Float64}, false, SVector{11, Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(TM), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tuple{Int64}}, ::Float64, ::Float64, ::SVector{3, Float64}, ::ODEInterface.OptionsODE)

Closest candidates are:
  radau(::Any, ::Real, ::Real, !Matched::Vector, ::ODEInterface.AbstractOptionsODE)
   @ ODEInterface C:\Users\Alex\.julia\packages\ODEInterface\RwRLn\src\Radau.jl:742

Stacktrace:
  [1] __solve(prob::ODEProblem{SVector{3, Float64}, Tuple{Float64, Float64}, false, SVector{11, Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(TM), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::radau{Nothing}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Vector{Any}; saveat::Vector{Float64}, verbose::Bool, save_everystep::Bool, save_on::Bool, save_start::Bool, timeseries_errors::Bool, dense_errors::Bool, callback::Nothing, alias_u0::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ODEInterfaceDiffEq C:\Users\Alex\.julia\packages\ODEInterfaceDiffEq\lGV1B\src\solve.jl:144
  [2] __solve(prob::ODEProblem{SVector{3, Float64}, Tuple{Float64, Float64}, false, SVector{11, Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(TM), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::radau{Nothing}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Vector{Any})
    @ ODEInterfaceDiffEq C:\Users\Alex\.julia\packages\ODEInterfaceDiffEq\lGV1B\src\solve.jl:1
  [3] __solve
    @ C:\Users\Alex\.julia\packages\ODEInterfaceDiffEq\lGV1B\src\solve.jl:1 [inlined]
  [4] #solve_call#22
    @ C:\Users\Alex\.julia\packages\DiffEqBase\G15op\src\solve.jl:511 [inlined]
  [5] solve_call
    @ C:\Users\Alex\.julia\packages\DiffEqBase\G15op\src\solve.jl:481 [inlined]
  [6] #solve_up#30
    @ C:\Users\Alex\.julia\packages\DiffEqBase\G15op\src\solve.jl:972 [inlined]
  [7] solve_up
    @ C:\Users\Alex\.julia\packages\DiffEqBase\G15op\src\solve.jl:945 [inlined]
  [8] #solve#28
    @ C:\Users\Alex\.julia\packages\DiffEqBase\G15op\src\solve.jl:882 [inlined]
  [9] solve(prob::ODEProblem{SVector{3, Float64}, Tuple{Float64, Float64}, false, SVector{11, Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(TM), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::radau{Nothing})
    @ DiffEqBase C:\Users\Alex\.julia\packages\DiffEqBase\G15op\src\solve.jl:872
 [10] top-level scope
    @ c:\Users\Alex\Desktop\dynamical-systems\Tsodyks Markram\

```

---

<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:** [July 18, 2023, 10:52am UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/18 "2023-07-18T10:52:29Z")

</div>

Make an mwe

---

<div class="post-metadata">

**Author:** ![Sergey\_Novak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sergey_novak/32/37716_2.png) [@Sergey\_Novak](https://discourse.julialang.org/u/Sergey_Novak)\
**Post date:** [July 18, 2023, 11:04am UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/19 "2023-07-18T11:04:35Z")

</div>

What is you mean by mwe?

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [July 18, 2023, 11:47am UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/20 "2023-07-18T11:47:21Z")

</div>

MWE == Minimum Working Example, i.e. a small, self-contained piece of code that exhibits the problem. You should be able to copy-paste a MWE into the REPL, run it, and trigger the error.

[Next page](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456.md?page=2)