# Poincare map DifferentialEquations.jl

**URL:** <https://discourse.julialang.org/t/poincare-map-differentialequations-jl/92308>\
**Category:** Modelling & Simulations\
**Created:** [December 30, 2022, 9:53am UTC](https://discourse.julialang.org/t/poincare-map-differentialequations-jl/92308 "2022-12-30T09:53:53Z")\
**Posts on this page:** 15\
**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:** [December 30, 2022, 9:53am UTC](https://discourse.julialang.org/t/poincare-map-differentialequations-jl/92308/1 "2022-12-30T09:53:53Z")

</div>

Hello, i have system:

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

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

---

<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:** [December 30, 2022, 10:58am UTC](https://discourse.julialang.org/t/poincare-map-differentialequations-jl/92308/2 "2022-12-30T10:58:48Z")

</div>

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

---

<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:** [December 30, 2022, 11:10am UTC](https://discourse.julialang.org/t/poincare-map-differentialequations-jl/92308/3 "2022-12-30T11:10:48Z")

</div>

Can you share example with this ?

---

<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:** [December 30, 2022, 11:16am UTC](https://discourse.julialang.org/t/poincare-map-differentialequations-jl/92308/4 "2022-12-30T11:16:30Z")

</div>

Look at DynamicalSystems.jl’s implementation.

---

<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:** [December 30, 2022, 4:05pm UTC](https://discourse.julialang.org/t/poincare-map-differentialequations-jl/92308/5 "2022-12-30T16:05:54Z")

</div>

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

---

<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:** [December 30, 2022, 4:38pm UTC](https://discourse.julialang.org/t/poincare-map-differentialequations-jl/92308/6 "2022-12-30T16:38:10Z")

</div>

> [@ChrisRackauckas](#):
>
> It’s not the most efficient way but it’s fine.

what is?

---

<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:** [December 30, 2022, 4:42pm UTC](https://discourse.julialang.org/t/poincare-map-differentialequations-jl/92308/7 "2022-12-30T16:42:35Z")

</div>

@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

---

<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:** [December 30, 2022, 4:46pm UTC](https://discourse.julialang.org/t/poincare-map-differentialequations-jl/92308/8 "2022-12-30T16:46:23Z")

</div>

> [@rveltz](#):
>
> what is?

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.

---

<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:** [December 30, 2022, 4:52pm UTC](https://discourse.julialang.org/t/poincare-map-differentialequations-jl/92308/9 "2022-12-30T16:52:42Z")

</div>

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](https://github.com/bifurcationkit/BifurcationKit.jl/blob/master/src/periodicorbit/ShootingDE.jl#L63) based on `VectorContinuousCallback`. I am interested in a better approach if any.

---

<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:** [December 30, 2022, 4:59pm UTC](https://discourse.julialang.org/t/poincare-map-differentialequations-jl/92308/10 "2022-12-30T16:59:29Z")

</div>

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.

> [@rveltz](#):
>
> The return map is coded [here](https://github.com/bifurcationkit/BifurcationKit.jl/blob/master/src/periodicorbit/ShootingDE.jl#L63) based on `VectorContinuousCallback`. I am interested in a better approach if any.

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

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [December 30, 2022, 5:01pm UTC](https://discourse.julialang.org/t/poincare-map-differentialequations-jl/92308/11 "2022-12-30T17:01:22Z")

</div>

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)](https://juliadynamics.github.io/DynamicalSystems.jl/dev/chaos/orbitdiagram/#ChaosTools.poincaresos) it tells you you can define a plane of any orientation by providing the normal vector to the plane.

> [@rveltz](#):
>
> I am surprised this is not provided in DiffEq

Well DynamicalSystems.jl builds on DifferentialEquations.jl so why duplicate effort 😃

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [December 30, 2022, 5:02pm UTC](https://discourse.julialang.org/t/poincare-map-differentialequations-jl/92308/12 "2022-12-30T17:02:27Z")

</div>

> [@Sergey\_Novak](#):
>
> I tried use DynamicalSystems.jl and it didn’t work out

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

---

<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:** [December 30, 2022, 5:17pm UTC](https://discourse.julialang.org/t/poincare-map-differentialequations-jl/92308/13 "2022-12-30T17:17:20Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [December 30, 2022, 11:32pm UTC](https://discourse.julialang.org/t/poincare-map-differentialequations-jl/92308/14 "2022-12-30T23:32:15Z")

</div>

> [@Sergey\_Novak](#):
>
> According to the documentation need define vector D+1 but i have vector of coefficients of plane for D-1.

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.

> [@Sergey\_Novak](#):
>
> i need plot Poincare map on plane, which is built on three points

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](https://math.stackexchange.com/questions/1874920/is-a-hyperplane-defined-by-four-points)

> [@Sergey\_Novak](#):
>
> Therefore, the Poincare mapping causes quite a lot of problems

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](https://en.wikipedia.org/wiki/Hyperplane#Affine_hyperplanes)

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](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)](https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/blob/main/test/poincaremap_tests.jl#L25-L45)

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?

---

<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:** [January 18, 2023, 9:02am UTC](https://discourse.julialang.org/t/poincare-map-differentialequations-jl/92308/15 "2023-01-18T09:02:24Z")

</div>

Thank you