# Find homoclinic DifferentialEquations.jl, DynamicalSystems.jl

**URL:** <https://discourse.julialang.org/t/find-homoclinic-differentialequations-jl-dynamicalsystems-jl/95800>\
**Category:** Modelling & Simulations\
**Created:** [March 9, 2023, 2:21pm UTC](https://discourse.julialang.org/t/find-homoclinic-differentialequations-jl-dynamicalsystems-jl/95800 "2023-03-09T14:21:45Z")\
**Posts on this page:** 18\
**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:** [March 9, 2023, 2:21pm UTC](https://discourse.julialang.org/t/find-homoclinic-differentialequations-jl-dynamicalsystems-jl/95800/1 "2023-03-09T14:21:45Z")

</div>

Hello, I’m trying to find a homoclinic orbit. To do this, I calculate a very long trajectory, and then calculate the minimum distance between all points of the trajectory and a fixed point. How can I optimize this?  
Fragment code for this:  
UPD: Correct calculate distance

```julia
tr = trajectory(ds, t, Δt = tstep;diffeq = integ_set)
trange = range(0.0, t, step = tstep)
fp, eigs, stable = fixedpoints(ds, box);
distance = zeros(length(tr))
for i in range(1, length(distance), step = 1)
    distance[i] = sqrt( (tr[i][1] - fp[1][1])^2 + (tr[i][2] - fp[1][2])^2 + (tr[i][3] - fp[1][3])^2)
end
minimum(distance)

```

---

<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:** [March 9, 2023, 7:26pm UTC](https://discourse.julialang.org/t/find-homoclinic-differentialequations-jl-dynamicalsystems-jl/95800/2 "2023-03-09T19:26:11Z")

</div>

You gave insufficient information on the vector field that defines your ODE. Is it a planar vector field or one in dimension d\>2.  
What is the dimension of the unstable and the stable eigenspace associated to the linearization of the vector field at its equilibrium?

In 2d, if p is a saddle equilibrium point, and u is the unit vector in the unstable direction, then the initial state (condition) should be p+epsilon\*u, for epsilon\>0, very small. In a “long” interval of integration, theoretically your orbit should return to the equilibrium point along the homoclinic orbit. But due to errors it can diverge from it. You have to perform a careful management of round off errors.  
Lately have been developed methods for computation of homoclinic orbits based on interval arithmetic. But DynamicalSystems.jl hasn’t implemented yet functions for computing these invariant sets.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [March 9, 2023, 7:31pm UTC](https://discourse.julialang.org/t/find-homoclinic-differentialequations-jl-dynamicalsystems-jl/95800/3 "2023-03-09T19:31:16Z")

</div>

You asked your question a bit too early… I released a new package for this [HclinicBifurcationKit.jl](https://github.com/bifurcationkit/HclinicBifurcationKit.jl) but it needs an update (a tag) of BifurcationKit.  
In any case, what you are asking for is not an easy matter.

---

<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:** [March 10, 2023, 5:30am UTC](https://discourse.julialang.org/t/find-homoclinic-differentialequations-jl-dynamicalsystems-jl/95800/4 "2023-03-10T05:30:50Z")

</div>

The system has a single state of equilibrium, it is saddle focus,its eigenvalues:

```julia
-9.853733057577516 + 0.0im
3.4937394347204993 - 16.249167707360805im
3.4937394347204993 + 16.249167707360805im

```

Code of system:

```julia
function TM(u, p, t)
    U(y) = p[8] + p[9] / ( 1.0 + exp( -50.0 * (y - p[7]) ) )
    σ(x) = 1.0 / ( 1.0 + exp( -20.0 * (x-p[6]) ) )
    du1 = (-u[1] + p[1] * log( 1.0 + exp( (p[5] * U(u[3]) * u[2] * u[1] + p[11] ) / (p[1]) ) ) ) / p[2]
    du2 = (1.0 - u[2])/p[3] - U(u[3])*u[2]*u[1]
    du3 = (-u[3])/p[4] + p[10] * σ(u[2])
    return SA[du1, du2, du3]
end;

```

I was advised to check the presence of a homoclinic by searching for the minimum distance between the equilibrium state and the trajectory points, and if it is quite small, then it is a homoclinic.

I was looking for such a distance for another system where there is a homoclinic and there the distance was quite small

---

<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:** [March 10, 2023, 5:48am UTC](https://discourse.julialang.org/t/find-homoclinic-differentialequations-jl-dynamicalsystems-jl/95800/5 "2023-03-10T05:48:31Z")

</div>

Okay, then I’ll be waiting for your package. Do you know the simplest methods of checking for homoclinic?

---

<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:** [March 10, 2023, 7:30am UTC](https://discourse.julialang.org/t/find-homoclinic-differentialequations-jl-dynamicalsystems-jl/95800/6 "2023-03-10T07:30:14Z")

</div>

I found your page and read about your work, I noticed that you work with models of neurons. I deal with them too, and that’s why I wanted to ask you for advice on how best to create activity mode maps. That is, how best to determine the burst and spike modes. I have already built such maps, but my algorithm is not universal

---

<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:** [March 10, 2023, 7:41am UTC](https://discourse.julialang.org/t/find-homoclinic-differentialequations-jl-dynamicalsystems-jl/95800/7 "2023-03-10T07:41:46Z")

</div>

Meanwhile you can try to get the homoclini orbit as follows: Since your saddle focus is opposite to the case of Shilnikov type saddle focus, illustrated in the attached image, you shouldd take as initial condition, `p+epsilon*vs`, `epsilon>0`  
(vs is the unit vector in the stable direction) and integrate backward in time. If after a long time of integration the orbit stays close to the equilibrium point, you can conclude that your computed orbit is an approximation of the homoclinic orbit. Since you don’t know in which region is included the homoclinic orbit with respect to the unstable plane, you should also try to compute the orbit starting at `p-epsilon*vs`, and integrate again backward in time.  
 ![silnik](https://global.discourse-cdn.com/julialang/original/3X/6/d/6d1e0d5c7a73ef8941583e6b6e85fe589989eb25.jpeg)  
Unfortunately most references on saddle-focus refer to Shilnikov phenomenon. I couldn’t find any case similar to your.

---

<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:** [March 10, 2023, 11:12am UTC](https://discourse.julialang.org/t/find-homoclinic-differentialequations-jl-dynamicalsystems-jl/95800/8 "2023-03-10T11:12:53Z")

</div>

Thanks, I try in backward + eps and - eps for all directions, and the trajectory goes to infinity. For direct time, when I use + eps for example for y, the trajectory is attracted to another attractor. As a result, I have two attractors, chaotic and regular.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/2/1/2103cc2b8980239aa2831b9babde775c9276e8c1.png)

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [March 10, 2023, 11:55am UTC](https://discourse.julialang.org/t/find-homoclinic-differentialequations-jl-dynamicalsystems-jl/95800/9 "2023-03-10T11:55:22Z")

</div>

do you have a time series that looks like a homoclinic orbit? Is there a Bogdanov-Takens biufrcation nearby? If not my package will not help you.

Indeed, in this case, one has to implement the homotopy method:

- E.J. Doedel, et. al. 1994. On locating connecting orbits
- E.J. Doedel, et. al. 1997. Successive continuation for locating connecting orbits

---

<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:** [March 10, 2023, 12:14pm UTC](https://discourse.julialang.org/t/find-homoclinic-differentialequations-jl-dynamicalsystems-jl/95800/10 "2023-03-10T12:14:29Z")

</div>

No, I have a stable fixed point, then a Hopf bifurcation, as a result of which a limit cycle is born, after which a Period-doubling cascade lead into a chaos. And I’m trying to find a homoclinic orbit in chaos

---

<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:** [March 10, 2023, 12:35pm UTC](https://discourse.julialang.org/t/find-homoclinic-differentialequations-jl-dynamicalsystems-jl/95800/11 "2023-03-10T12:35:59Z")

</div>

You said you calculated in all directions, but it’s important to specify which direction, the stable direction or one of the unstable directions ( the unstable space is 2d). As long as the real eigenvalue is simple and negative, it follows that the stable manifold is one-dimensional. If vs is a unit eigenvector in the space corresponding to \lambda\in\mathbb{R}, \lambda\<0, the trajectory in the forward direction tends to the equilibrium point, while in the backward direction along the stable manifold which eventually continues as a homoclinic orbit.

---

<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:** [March 10, 2023, 12:45pm UTC](https://discourse.julialang.org/t/find-homoclinic-differentialequations-jl-dynamicalsystems-jl/95800/12 "2023-03-10T12:45:44Z")

</div>

With a high probability I did everything wrong, so I will attach the code

```julia
t = 5000.0; tstep = 0.001; Tt = 500.0
trange = range(0.0, t, step = tstep)
integ_set = (alg = RK4(), adaptive = false, dt = tstep);

const τ = 0.013; const τD = 0.080; const τy = 3.3; const J = 3.07; const β = 0.300
const xthr = 0.75; const ythr = 0.4
const α = 1.58; const ΔU0 = 0.305;

p = SA[α, τ, τD, τy, J, xthr, ythr, 0.3, ΔU0, β, -1.599305]

ds = ContinuousDynamicalSystem(TM, SA[0.7935673221978937, -0.0011744045682187255, -4.56582255633564], p)
tr = trajectory(ds, t, Δt = tstep; diffeq = integ_set)
fp, eigs, stable = fixedpoints(ds, box);

eps = 1e-6
Efp, xfp, yfp = fp[1]
point = SA[Efp-eps, xfp, yfp]

ds1 = ContinuousDynamicalSystem(TM, point, p)
tr1 = trajectory(ds1, -t, Δt = -tstep; diffeq = integ_set);

```

I tried to subtract eps alternately from each coordinate of the equilibrium state

---

<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:** [March 10, 2023, 3:27pm UTC](https://discourse.julialang.org/t/find-homoclinic-differentialequations-jl-dynamicalsystems-jl/95800/13 "2023-03-10T15:27:10Z")

</div>

Once you have the fixed point (equilibrium point), `fixp`, and \lambda\<0 as the real eigenvalue of the Jacobian evaluated at `fixp`, calculate the eigenspace corresponding to this eigenvalue. Let v\_s be a vector in this eigenspace (||v\_s||=1).  
The backward trajectory starting at fixp+\epsilon \*v\_s, computed carefully, is the stable manifold associated to your saddle focus.  
It seems that in your code the stable direction, v\_s, is not involved.

---

<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:** [March 10, 2023, 4:10pm UTC](https://discourse.julialang.org/t/find-homoclinic-differentialequations-jl-dynamicalsystems-jl/95800/14 "2023-03-10T16:10:28Z")

</div>

What did you mean by “eigenspace” is eigenvectors?  
If I understood correctly, did you mean this?

 ![image](https://global.discourse-cdn.com/julialang/original/3X/f/b/fb7e87e31f8ff1ed838bd58242e65f6ffc1e750e.png)

---

<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:** [March 10, 2023, 6:00pm UTC](https://discourse.julialang.org/t/find-homoclinic-differentialequations-jl-dynamicalsystems-jl/95800/15 "2023-03-10T18:00:17Z")

</div>

If \lambda is an eigenvalue of a n\times n matrix, J, (here J is your Jacobian at the equilibriun point), the corresponding eigenspace is S\_\lambda=\{v\in \mathbb{R}^n\:|\: Jv=\lambda v\}, i.e. is the set of corresponding eigenvectors. Since in your case \lambda \<0 is a simple eigenvalue (i.e. a simple root of the characteristic polynomial), the subspace (eigenspace) S\_\lambda has dimension one.  
In your data posted above, v\_s has the coordinates in the first column of vectors. It defines the basis in the 1-dimensional subspace S\_\lambda.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [March 12, 2023, 5:54pm UTC](https://discourse.julialang.org/t/find-homoclinic-differentialequations-jl-dynamicalsystems-jl/95800/16 "2023-03-12T17:54:21Z")

</div>

On a side note, I consider indexing parameters with an array as bad practice. So easy to get it wrong, you should use named tuples

---

<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:** [March 13, 2023, 5:45am UTC](https://discourse.julialang.org/t/find-homoclinic-differentialequations-jl-dynamicalsystems-jl/95800/17 "2023-03-13T05:45:38Z")

</div>

Previously, I parsed parameters and phase variables in function, but then I noticed that if this is not do the time during which the system integrates is reduced. I have not tried using named tuples, I will have to see how this affects the integration time.

The thing is that earlier I calculate a two-parameter map of the spectrum of Lyapunov exponents, map of modes and orbit diagram, and all this took a decent time. Therefore, had to optimize the function as much as possible to the detriment of code readability

UPD: I still don’t know how to write fast and good code in Julia, and the algorithms I writed and used to inherit initial conditions and determine modes are not very good from an optimization point of view.

---

<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:** [March 13, 2023, 5:46am UTC](https://discourse.julialang.org/t/find-homoclinic-differentialequations-jl-dynamicalsystems-jl/95800/18 "2023-03-13T05:46:24Z")

</div>

Thanks, I think I understand now