Find homoclinic DifferentialEquations.jl, DynamicalSystems.jl

Specific Domains Modelling & Simulations
1

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

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

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.

3

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

2 Likes
4

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

-9.853733057577516 + 0.0im
3.4937394347204993 - 16.249167707360805im
3.4937394347204993 + 16.249167707360805im

Code of system:

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

5

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

6

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

7

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
Unfortunately most references on saddle-focus refer to Shilnikov phenomenon. I couldn’t find any case similar to your.

1 Like
8

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
image1200×900 203 KB

9

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
10

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

11

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.

12

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

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

13

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.

14

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

image
image776×291 13.5 KB

15

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.

16

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

17

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.

18

Thanks, I think I understand now