4

I am trying to using the DynamicalSystems julia package to plot an orbit diagram of a discrete dynamical system governed by the following equation:

y_{n+1} = 1- by_n^2

I'm able to get it to plot however I can't seem to get the initial condition to work, I am trying to use the initial condition of y_0 = -0.68.

Here is my code:

using DynamicalSystems
using PyPlot


function eom(dx, x, p, n)
    dx[1] = 1 - p[1] * x[1]^2
end

ds = DiscreteDynamicalSystem(eom, [-0.68], [2.0, 0.0])

i = 1
pvalues = 0.0:0.001:2.0
n = 2000
Ttr = 2000
p_index = 1
output = orbitdiagram(ds, i, p_index, pvalues; n = n,
                      Ttr = Ttr, u0 = get_state(ds))

L = length(pvalues)
x = Vector{Float64}(undef, n*L)
y = copy(x)
for j in 1:L
    x[(1 + (j-1)*n):j*n] .= pvalues[j]
    y[(1 + (j-1)*n):j*n] .= output[j]
end

figure()
PyPlot.title("total points: $(L*n)")
plot(x, y, ls = "None", ms = 0.5, color = "black", marker = "o", alpha = 0.05)
xlim(pvalues[1], pvalues[end]);
xlabel("\$b\$"); ylabel("\$y\$")
tight_layout()

And here is the figure it generates: enter image description here

If possible, could someone also tell me if there is a way to get the bifurcation points using DynamicalSystems.

Cavenfish
  • 366
  • 2
  • 13

1 Answers1

2

Sorry, I was not thinking straight. It already is taking the initial condition into account but the bifurcation diagram is independent of the initial condition. In terms of finding the bifurcation values, i wrote a code to check when y_{n+2} = y_{n} but y_{n+1} != y_{n} Which gives the first bifurcation value.

Cavenfish
  • 366
  • 2
  • 13