Basins of attractros for 6D system DynamicalSystems.jl

Sergey_Novak · August 15, 2022, 2:21pm

I trye found basins for attractors, but i have some problems.
First problem is out of memmory when i use big length for
xg1 = yg1 = zg1 = xg2 = yg2 = zg2 = range(-6.0, 6.0; length = 10)
If i using for example length = 100 i have error out of memory. Can i fix this or i need using poincare map?

Second problem is projection basins and attractors on 3D. I think I can do a projection of attractors on 3D, but I have a strange result
Image:

For this i using this code:

xg1 = yg1 = zg1 = xg2 = yg2 = zg2 = range(-6.0, 6.0; length = 10)
mapper = AttractorsViaRecurrences(ds_HR, (xg1, yg1, zg1, xg2, yg2, zg2))
basins, attractors = basins_of_attraction(mapper; show_progress = true)


function scatter_attractors!(ax, attractors)
    
    for k ∈ keys(attractors)
        x1, y1, z1, x2, y2, z2  = columns(attractors[k])
        scatter!(ax, x1, y1, z1;
            color = Cycled(k),
            strokewidth = 3, strokecolor = :white
        )
    end
end

fig = Figure()

ax = fig[1,1] = Axis3(fig; title = "found attractors")

scatter_attractors!(ax, attractors)

fig

Also i don’t understand how need modification variable basins for projection on 3D.

Full code:

function sigma(x)
    return @fastmath 1.0 / ( 1.0 + exp( -10.0 * ( x  - ( - 0.25 ) ) ) )
end

function HR(u, p, t)
        
    a, b, c, d, s, xr, r,  I, vs, k1, k2, el_link  = p
    x1, y1, z1, x2, y2, z2 = u
    
    du1 = y1 + b * x1 ^ 2 - a * x1 ^3 - z1 + I - k1 * ( x1 - vs ) * sigma(x2) + el_link * ( x2 - x1 )
    du2 = c - d * x1 ^2 - y1
    du3 = r * ( s * ( x1 - xr ) - z1 )
    
    du4 = y2 + b * x2 ^ 2 - a * x2 ^3 - z2 + I - k2 * ( x2 - vs ) * sigma(x1) + el_link * ( x1 - x2 )
    du5 = c - d * x2 ^2 - y2
    du6 = r * ( s * ( x2 - xr ) - z2 )
    return SVector(du1, du2, du3,
                    du4, du5, du6)
end

a = 1.; b = 3.; c = 1.; d = 5.;
xr = -1.6; r = 0.01; s = 5.; I = 4.; xv = 2.;
k1= -0.17; k2 = -0.17;
k = 1.0
condition = SA[-1.5, 0.0, 0.0, -2.5, 0.0, 0.0]
# condition = SA[-0.5, -0.5, -0.5, -0.5, 0-0.5, -0.5]

p = SA[a, b, c, d,
        s, xr, r, I, xv, k1, k2, k]
p = [a, b, c, d, s, xr, r, I, xv, k1, k2, k]

ds_HR = ContinuousDynamicalSystem(HR,  condition, p )

xg1 = yg1 = zg1 = xg2 = yg2 = zg2 = range(-6.0, 6.0; length = 10)
mapper = AttractorsViaRecurrences(ds_HR, (xg1, yg1, zg1, xg2, yg2, zg2))
basins, attractors = basins_of_attraction(mapper; show_progress = true)

function scatter_attractors!(ax, attractors)
    
    for k ∈ keys(attractors)
        x1, y1, z1, x2, y2, z2  = columns(attractors[k])
        scatter!(ax, x1, y1, z1;
            color = Cycled(k),
            strokewidth = 3, strokecolor = :white
        )
    end
end

fig = Figure()

ax = fig[1,1] = Axis3(fig; title = "found attractors")

scatter_attractors!(ax, attractors)

fig

Sergey_Novak · August 15, 2022, 2:29pm

If i using poincare map i can fix fisrt problem, but can i fix this without poincare map?
Also i have this progress bar when using length = 50

Datseris · August 15, 2022, 2:58pm

what’s the strange thing?

Datseris · August 15, 2022, 2:59pm

You have to wait a couple of weeks for our new sparse algorithm to be published!

Sergey_Novak · August 15, 2022, 3:04pm

I don’t understand why i have so many dots. I have 3 attractors

Sergey_Novak · August 15, 2022, 3:07pm

Can i reshape basins from 6D to 3D or 2D and then plot him?

Datseris · August 15, 2022, 3:14pm

Yeah and your plot has 3 colors, because you have 3 attractors. Attractors may have “many points”, i.e., not be fixed points. I’m a bit confused by the question.

You can project them if you want, after computing them, but it makes no sense to “reshape” the basins into 3D. What sense does it make to reshape a sphere into a circle?

Given your questions, perhaps you should consider going through a couple of textbooks of nonlinear dynamics to sharpen your understanding in these topics. Unfortunately, it is not possible to get this understanding only by running code

Sergey_Novak · August 15, 2022, 3:25pm

I’m a little confused why the attractors had so many points, now I understand thanks.
I’m trying to see basins and attractors in 3D and 2D, I just don’t understand how to do it.

Sergey_Novak · August 15, 2022, 3:52pm

Also I found a small typo in the documentation code.
In sectiond 3D in code for ploting 3d Image using variable bsn but it not declared. I think the need replace bsn on basins
Full code from section:

using GLMakie
fig = Figure()
display(fig)
ax = fig[1,1] = Axis3(fig; title = "found attractors")
cmap = cgrad(:dense, 6; categorical = true)

for i in keys(attractors)
    tr = attractors[i]
    markersize = length(attractors[i]) > 10 ? 2000 : 6000
    marker = length(attractors[i]) > 10 ? :circle : :rect
    scatter!(ax, columns(tr)...; markersize, marker, transparency = true, color = cmap[i])
    j = findfirst(isequal(i), bsn)
    x = xg[j[1]]
    y = yg[j[2]]
    z = zg[j[3]]
    tr = trajectory(ds, 100, SVector(x,y,z); Ttr = 100)
    lines!(ax, columns(tr)...; linewidth = 1.0, color = cmap[i])
end

a = range(0, 2π; length = 200) .+ π/4

record(fig, "cyclical_attractors.mp4", 1:length(a)) do i
    ax.azimuth = a[i]
end

Sergey_Novak · August 15, 2022, 5:55pm

I found one of your articles. By my question I meant how to get an image similar to (h) with dots of attractor for example from (c).

Also i needed images how this

Now it has become a little clearer to me how to get them, sorry for spaming

Datseris · August 15, 2022, 7:15pm

Well you can do this if your system is 3D or 2D. Otherwise you have to project the basins and/or the attractors into a lower dimensional space. This is not reshaping, nor is it an invertible operation. You need to know the system to know when this is possible and when it retains some accuracy. In many systems the attractors can’t be separated by projection.

Sounds like a great opportunity for a pull request!

(h) is just a slice of 6D basins. Like plotting basins[1, 1, 1, 1, :, :]. There are no attractors plotted there, because you simply can’t project the attractors into 2D space for this system.

No worries, but what I said above is for your benefit, to spend some time with a textbook or two. It will raise your intuition levels which is the most important thing. Running code is easy, understanding it is the hard part.

Sergey_Novak · August 16, 2022, 7:16am

Thanks, thanks to the code from the article, it became clearer how to do this.
But some code snippets worked a little incorrectly. These fragments code are related to the color of the attractors and they were also found in the documentation.

For example this fragment code from section 3D basins:

fig = Figure()

ax = fig[1,1] = Axis3(fig; title = "found attractors")
cmap = cgrad(:dense, 6; categorical = true)

for i in keys(attractors)
    
    tr = attractors[i]
    markersize = length(attractors[i]) > 10 ? 2000 : 6000
    marker = length(attractors[i]) > 10 ? :circle : :rect
    
    scatter!(ax, columns(tr)...; markersize, marker, transparency = true, color = cmap[i])
    j = findfirst(isequal(i), bsn)
    x = xg[j[1]]
    y = yg[j[2]]
    z = zg[j[3]]
    tr = trajectory(ds, 100, SVector(x,y,z); Ttr = 100)
    lines!(ax, columns(tr)...; linewidth = 1.0, color = cmap[i])
end

I fix this deleting line code markersize = length(attractors[i]) > 10 ? 2000 : 6000

Also i have problem with this visualization on my system. I modification this code

fig = Figure()
ax = Axis3(fig[1,1])
cmap = cgrad(:dense, 8; categorical = true)

for i in keys(attractors)
    
    tr = attractors[i]
    
    marker = length(attractors[i]) > 10 ? :circle : :rect
    
    idxs = [1,2,3]
    
    x,y,z = columns(tr)[idxs]
    
    scatter!(ax, x,y,z; marker, transparency = true, color = cmap[i])
    
    
    tr = trajectory(ds_HR, 5000, tr[1], Ttr = 1000;
                    diffeq = (alg = AutoVern9(Rodas5()), abstol = 1e-11, reltol = 1e-11,  maxiters = 10000000) )
    
    x,y,z = columns(tr)[idxs]
    
    lines!(ax, x,y,z; linewidth = 1.0, color = cmap[i])

end

fig

and get this result

I don’t see on this images trajectory from anothers attractors, for example for blue points. This may be due to the fact that they cannot be projected onto a given subspace? Or this error with visualization?

If i using more time and Ttr, for example 8000 and 4000 i have this images:

fragment of code for search basins and attractors

xg1 = yg1 = zg1 =  xg2 = yg2 = zg2 = range(-5.0, 5.0; length = 15)

mapper = AttractorsViaRecurrences(ds_HR, (xg1, yg1, zg1, xg2, yg2, zg2),
                                    diffeq = (alg = AutoVern9(Rodas5()),
                                                abstol = 1e-11, reltol = 1e-11, dense = false,  maxiters = 10000000))
basins, attractors = basins_of_attraction(mapper)
attractors