Poincare map DifferentialEquations.jl

Specific Domains Modelling & Simulations
1

Hello, i have system:

function model(u, p ,t)
    
    E, x, u_, y = u
    τ, α, τ_D, J, U0, ΔU0, τ_y, β, xthr, ythr, I0 = p
    
    U(y, U0, ΔU0, ythr) = U0 + ΔU0 / ( 1 + exp( -50 * ( y - ythr ) ) )
    σ(x, xthr)= 1 / (1 + exp(-20 * (x - xthr)));
    
    du1 = (-E + α * log( 1 + exp( ( J * u_ * x * E + I0 ) / α ) )) / τ
    du2 = ( 1 - x ) / τ_D - u_ * x * E
    du3 = (U(y, U0, ΔU0, ythr) - u_) / τ_F + U(y, U0, ΔU0, ythr) * ( 1 - u_ ) * E
    du4 = -y /  τ_y + β * σ(x, xthr)
    
    return SVector(du1, du2, du3, du4)
    
end;

and i need plot Poincare map on plane, which is built on three points, which are the equilibrium states of the system. I have coordinates equilibrium points and i plot plane as follow

function plane_v2(x, y, fp1, fp2, fp3)
    """
    E 1
    x 2
    u 3
    y 4
    """
    x0 = fp1[2]
    x1 = fp2[2]
    x2 = fp3[2]
    
    y0 = fp1[3]
    y1 = fp2[3]
    y2 = fp3[3]
    
    z0 = fp1[1]
    z1 = fp2[1]
    z2 = fp3[1]
    
    a = ( y1 - y0 ) * ( z2 - z0 ) - ( y2 - y0 ) * ( z1 - z0 )
    b = ( z1-z0 ) * ( x2  -x0 ) - ( x1 - x0 ) * (z2 - z0)
    c = ( x1 - x0 ) * ( y2 - y0 ) - ( x2 - x0 ) * ( y1 - y0 )
    
    z =  z0 .+ ( -( x .- x0 ) * a - (y .- y0) * b ) / c
    return z
    
end

x_range = range( 0.0, 1.0, length = 1000 )
u_range = range( 0.0, 1.0, length = 1000 )
E_range = plane_v2( x_range, u_range, first_eq, second_eq, third_eq  )

how do I properly use callbacks to get a poincare mapping? I tried use DynamicalSystems.jl and it didn’t work out

2

Just make a continuouscallback at the roots that you want. It’s not the most efficient way but it’s fine.

3

Can you share example with this ?

4

Look at DynamicalSystems.jl’s implementation.

5

In code of poincare map DynamicalSystems.jl i can not find callback

6

what is?

7

@Datseris Hello, can you share link on code from github with callback for Poincare map. I know coefficients from equation of plane for 3d projection and i want try use this for plot Poincare map. That is, I have a plane in a 3-dimensional projection and I want to use it to plota Poincare mapping

8

A discrete callback with NonlinearSolve rootfinding. It’s similar but not the same because ContinuousCallbacks will pull back to the time of the of the root. If you you cross multiple roots in an interval, the ContinuousCallback will pull back to the first one, then do a step, rootfind, pullback, step, rootfind, etc. until it’s gone through all of the roots. However, this is unnecessary if your callback doesn’t modify the ODE at all: you might as well do a big steps, find all of the roots, and then take the next big step. Since in theory a special DiscreteCallback is better, though in practice the ContinuousCallback way is much simpler.

9

I am surprised this is not provided in DiffEq, this is such common functionality that “everybody” wants and if you know better than anybody how to do it right, that’s a great benefit for all of us… ;). Is it documented?

I am saying this because in BifurcationKit, the Poincaré shooting is a little brittle. The return map is coded here based on VectorContinuousCallback. I am interested in a better approach if any.

10

We could add one to the DiffEqCallbacks library, it would fit in as one of those premade callbacks. But @Datseris 's DynamicalSystems.jl has been a good Poincure map world for so long that I haven’t cared to implement this version.

Definitely the discrete form would do better since VectorContinuousCallback doesn’t handle simultaneous events, among other things (I may fix that rather soon though).

11

Poincare map in DynamicalSystems.jl is the most efficient version as far as I know. We have compared with several variants and several callback versions. That was 3 years ago though, maybe things changed. But I think not, because our Poincare code is not a general code, it can only work for poincare planes, utilizing NLD knowledge. I have to admit we have not tested a variant where the find of the root happens via Optim or other optimization packages, but we use the standard Roots.jl interface. Howeever we search very close to the root due to system knowledge. Also we have benchmarked that simply using step! is faster than using callbacks.

You can define a custom plane @Sergey_Novak , why don’t you have a look at the documentation of the package? Orbit Diagrams & PSOS · DynamicalSystems.jl (juliadynamics.github.io) it tells you you can define a plane of any orientation by providing the normal vector to the plane.

Well DynamicalSystems.jl builds on DifferentialEquations.jl so why duplicate effort :smiley:

12

also please understand that this isn’t helpful. Say exactly what didn’t work and exactly what you expected to happen…

13

According to the documentation need define vector D+1 but i have vector of coefficients of plane for D-1.

The plane itself in 3d is constructed in such a way that all equilibrium states lie on it, and there are only three of them. That is, I received the plane by three points and by their three coordinates. As I understand it, in the 4th dimension I need a 4 points, and i I can’t get another dot.

Therefore, the Poincare mapping causes quite a lot of problems

14

I can’t understand this sentence unfortunately. In any case you indeed need D+1 real numbers to uniquely define a hyperplane in D dimensions. See the wikipedia page on hyperplanes.

Your dynamical system is 4D. So the state space is 4D. You can’t define a hyper plane in 4 dimensions using only three points, you need at least four. See geometry - Is a hyperplane defined by four points? - Mathematics Stack Exchange

It isn’t the poincare map that is causing the problems. It is your lack of knowledge and/or understanding of hyperplanes. Have a look at the wikipedia page on hyper planes. There you will find the equation that I quote in the documentation string of poincaresos , that says you need to provide a vector of D+1 real numbers to define a hyperplane in D-dimensional space. Here it is: Hyperplane - Wikipedia

In fact, in the very first published tutorial of DynamicalSystems.jl , published 4 years ago, I show exactly what you want to do here, making a Poincare Section on a hyperplane connecting the (unstable) fixed ponts of a dynamical system. Here you are: https://www.youtube.com/live/13hqE_1a158?feature=share&t=2277

and here is the code that does the same thing from the test suite. The comment in the code block literally says # create hyperplane spanned by these two points. DynamicalSystemsBase.jl/poincaremap_tests.jl at main · JuliaDynamics/DynamicalSystemsBase.jl (github.com)

Do you see how gis_plane(μ) = [cross(Np(μ), Nm(μ))..., 0] is a vector of four numbers even though the Gissinger system I use it for is 3-dimensional?

2 Likes
15

Thank you