In [3]:
#=
Compute first few periodic orbits of the Rossler system.
=#
#Call python's plotting routine
using PyPlot
#call ODE package; ODE.jl is similar to matlab functions, but sundials is an alternative
using ODE
# call interpolation package for return maps
using Interpolations
In [10]:
function rosslermap(t, r)
# Constants
a = 0.1
b = 0.1
c = 14
# store r as 3-tuple
(x, y, z) = r
# The Rossler equations
rossler1 = -y - z
rossler2 = x + a * y
rossler3 = b + z * (x - c)
#return image of rossler map at r as vector
rossler = [rossler1; rossler2; rossler3]
return rossler
end
Out[10]:
In [11]:
function rosslerintegrate(initialcondition, endtime, timestep)
timespan = collect(0 : timestep: endtime)
timespan, pos = ode45(rosslermap, initialcondition, timespan)
x = map(v -> v[1], pos)
y = map(v -> v[2], pos)
z = map(v -> v[3], pos)
return (timespan, x, y, z)
end
Out[11]:
In [12]:
rosslerintegrate([1.0,1.0,1.0], 10.0, 0.01)
Out[12]:
In [13]:
function rosslerplot(initialcondition, endtime, timestep)
clf()
timespan = collect( 0 : timestep : endtime)
(timespan, x, y, z) = rosslerintegrate(initialcondition, endtime, timestep)
plot3D(x, y, z);
xlabel("x")
ylabel("y")
zlabel("z")
axis("tight")
end
Out[13]:
In [14]:
function rosslerreturn(deta, dt, tf)
#construct poincare section as matrix, with each point on the section a column vector
eta = collect(0: deta: 1.0)
poincare = zeros(3, length(eta))
poincare[2,:] = -22.0 + 14.0 * eta
time = collect(0 : dt : tf)
# compute trajectory starting at each point on poincare section,
# check each pair of successive points on trajectory to find intersections
# with poincare section
#
feta = zeros(2, round(Integer, tf/dt))
v = 1
for j in collect(1 : 1 : length(eta)-1)
#integrate starting at each point on poincare section
r0 = poincare[:,j]
~, x, y, ~ = rosslerintegrate(r0, tf, dt)
for i in collect(1 : 1 :length(y)-1)
if ((x[i] < 0) && (x[i+1] > 0))
feta[1, v] = y[1]
feta[2, v] = (y[i] + y[i+1])/2
v += 1
end
end
end
fetatrunc = zeros(2, findlast(feta[2,:]))
for j in collect(1 : 1 : findlast(feta[2,:]))
fetatrunc[1, j] = (feta[1, j] + 22.0) / 14.0
fetatrunc[2, j] = (feta[2, j] + 22.0) / 14.0
end
return fetatrunc
end
Out[14]:
In [9]:
function rosslerreturnplot(deta, dt, tf)
y = rosslerreturn(deta, dt, tf)
clf()
scatter(y[1,:], y[1,:])
hold(true)
scatter(y[1,:], y[2,:])
end
Out[9]:
In [16]:
function rosslerperiodicplot(deta, dt, tf1, tf2)
# find fixed points of return map
# use enough eta values to allow accurate interpolation
y = rosslerreturn(deta, dt, tf1)
itp = interpolate(y, BSpline(Quadratic(Reflect())),
OnCell())
# build diagonal; intersections of diagonal and return map curve
# are fixed points
ydiag = [y[1,:]; y[1,:]]
# compute residual of return map and diagonal
A = abs( y[2,:] - ydiag[2,:]).^2/(ydiag[2,:]).^2
# find minimum; this is a good guess for the fixed point
# of the return map, which corresponds to the intersection
# of a periodic orbit of the rossler system with the Poincare section
fp, j_min = findmin(A)
# convert back to standard spatial coordinates
fp = -22.0 + 14.0 * fp
# integrate starting at fixed point
# this is periodic orbit
r0 = [0; fp; 0]
clf()
rosslerplot(r0, tf2, dt)
end
Out[16]:
In [11]:
#= scratch
function rosslerjacobian(r)
# store r as 3-tuple
(x, y, z) = r
# compute Jacobian of rossler map at r=(x,y,z)
Df = [ 0 -1 -1;
1 a 0;
z 0 x-c; ]
return Df
end
=#
#function rossler_return(eta)
#zero = [0; 0; 0]
#alpha = c / 2.0 / a
#bet = 4 * a * b / c /c
#disc = sqrt(1 - bet)
#gammaplus = alpha * (1 + disc)
#gammaminus = alpha * (1 - disc)
#xstarhat = [ a; -1; 1]
#xstarplus = gammaplus * xstarhat
#xstarminus = gammaminus * xstarhat
#compute jacobian at equilibria
#Dfzero = rossler_jacobian(zero)
#Dfplus = rossler_jacobian(xstarplus)
#Dfminus = rossler_jacobian(xstarminus)
#compute eigenvalues and eigenvectors of jacobian at each
#eigvalszero, eigvecszero = eigfact(Dfzero)
#eigvalsplus, eigvecsplus = eigfact(Dfplus)
#eigvalsminus, eigvecsminus = eigfact(Dfminus)
In [17]:
@time rosslerperiodicplot(0.01,0.01,10.0)
In [ ]: