Hello, I need to repeat the following orbit diagram, where π₯||πππ₯ is maximums from π₯|| = x1 + x2. I know how calculate him from two variable but i donβt know how ploting him, examples in documentation only for one variable.
I using next code, for calculating from two variable use i = [1, 4]
function sigma(x)
return 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, k, el_link = p
x1, y1, z1, x2, y2, z2 = u
du1 = y1 + b * x1 ^ 2 - a * x1 ^3 - z1 + I - k * ( 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 - k * ( 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.;
k = 0.0532;
el_link = 0.0
condition1 = SA[-1.5, 0.0, 0.0, -2.5, 0.0, 0.0]
condition2 = 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, k, el_link]
p = [a, b, c, d, s, xr, r, I, xv, k, el_link]
ds_HR1 = ContinuousDynamicalSystem(HR, condition1, p )
ds_HR2 = ContinuousDynamicalSystem(HR, condition2, p );
ds_HR1
ds = ds_HR1
pvalues = range(0.050, stop = 0.056, length = 5000)
i = 1
plane = (2, 0.0)
tf = 5000.0
tr = 5000.0
p_index = 10
output = produce_orbitdiagram(ds, plane, i, p_index, pvalues,
tfinal = tf, Ttr = tr, printparams = true;
diffeq = (alg = AutoVern9(Rodas5()), abstol = 1e-11, reltol = 1e-11, maxiters = 10000000)
);
fig = Figure()
ax = Axis(fig[1,1]; xlabel = L"k", ylabel = L"x_{2}")
for (j, p) in enumerate(pvalues)
scatter!(ax, fill(p, length(output[j])), output[j];
color = ("black", 0.5), markersize = 1.5
)
end
fig
