# based on Sage interact examples. Marshall Hampton. Creative Commons by SA
# Staffan Liljegren simplified and added more 1d maps and made the Cython code
# more flexible to pass maps. Plots the iterated map and the bifurcation diagram for 1d
# chaotic maps with flip/periodic doubling bifurcations. Added "time series" plots
|
|
First we define some functions in Cython, so people can subclass Function and add new 1D maps. The orbit function iterates the function f N times and returns the value. It is used for plotting bifurcation and can also be used if you want to; skip the transient part. trajectory is a function that iterates N times and keeps the itterates in a list which is returned- This is used for plotting time series of (x_n,n): Both are used in more advanced analysis for for calculating Liapunov value and Feigenbaums invariant \delta_i
%auto
%cython
cdef class Function:
cpdef double evaluate(self, double x,double k) except *:
return 0
# use values 0<=k<=4, so map x <= 1
cdef class LogisticMap(Function):
cpdef double evaluate(self, double x, double k) except *:
return (k*x*(1-x))
# 0<=k<=2 (try k= -0.25 then it becomes a saddle node bifurcation)
cdef class SquareMap(Function):
cpdef double evaluate(self, double x, double k) except *:
return (k - x**2)
# -4.3<k<0
cdef class DoubleUndoubleMap(Function):
cpdef double evaluate(self, double x, double k) except *:
return (k + 11.5*x/(1+x**2))
cpdef double orbit(Function f, double k,long N,double x0):
cdef double x = x0
cdef long i
cdef double s, dx
if f is None:
raise ValueError("function cannot be None")
for i in range(N):
x = f.evaluate(x,k)
return x
cpdef trajectory(Function f,double k,long N, double x0):
cdef double x = x0
xvals = []
cdef long i
if f is None:
raise ValueError("function cannot be None")
for i in range(N):
x = f.evaluate(x,k)
xvals.append(x)
return xvals
#logistics map equilibrium and the fixed points
var('x,k')
lmap = k*x*(1-x)
fps = solve(x == lmap,x)
eqs = solve(0 == lmap,x)
fps,eqs,derivative(lmap,x)
([x == (k - 1)/k, x == 0], [x == 0, x == 1], -(x - 1)*k - k*x) ([x == (k - 1)/k, x == 0], [x == 0, x == 1], -(x - 1)*k - k*x) |
Traceback (click to the left of this block for traceback) ... SyntaxError: invalid syntax Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "_sage_input_29.py", line 10, in <module>
exec compile(u"print _support_.syseval(python, u'(orbit(LogisticMap(),k,2,0) -2', __SAGE_TMP_DIR__)" + '\n', '', 'single')
File "", line 1, in <module>
File "/sagenb/sage_install/sage-4.7/devel/rkirov-flask/sagenb/misc/support.py", line 473, in syseval
return system.eval(cmd, sage_globals, locals = sage_globals)
File "/sagenb/sage_install/sage-4.7/local/lib/python2.6/site-packages/sage/misc/python.py", line 53, in eval
eval(compile(s, '', 'exec'), globals, globals)
File "", line 3
(orbit(LogisticMap(),k,2,0) -2
^
SyntaxError: invalid syntax
|
t = trajectory(LogisticMap(),3.0,100,.1)
#t[1::] # get every other item, starting with the second
list_plot(zip(t,t[1::]))
|
|
|
|
# uncomment next if you want to see how much time it takes
#%time
# leave
lp = [list_plot(trajectory(LogisticMap(),4.0,60,1e-9*random()),plotjoined =True,thickness=0.4) for i in range(9)]
graphics_array(lp,3,3).show()
|
lp = [list_plot(trajectory(LogisticMap(),k,200,0.5)) for k in srange(2.0,4.0,0.1)]
a = animate(lp, figsize=[4,4],xmin=0.0,xmax=100.,ymin=0.0, ymax=1.0)
a.show(delay=200)
|
#just for one value of k sofar - correct? no...
def liapunov(f,k):
x0= random()
transient = orbit(f,k,200,x0)
lp = orbit(f,3.0,20000,transient)
lyap = ln(abs(k-2*k*lp))
return lyap
|
|
def cobweb(f, k, start, iterations = 10, xmin = 0.0, xmax = 1.0,skip_transient=False):
# plot the function we are given in blue
basic_plot = plot(f.evaluate(start,k),xmin = xmin, xmax = xmax ,rgbcolor = (0,0,1),thickness= 0.6,axes_labels=['$x_n$','$x_{n+1}$'])
# add diagonal which shows fixed points where f crosses
basic_plot = basic_plot + line([(xmin,xmin),(xmax,xmax)],thickness=0.3)
# now we show iterations as red cobweb diagrams
if skip_transient:
start = orbit(f,k,100,start)
iter_list = [[start,start]]
else:
iter_list = [[start,0]]
current = start
for i in range(iterations):
iter_list.append([current,f.evaluate(current,k)])
current = f.evaluate(current,k)
iter_list.append([current,current])
cobweb = line(iter_list, rgbcolor = (1,0,0),thickness=0.4)
return basic_plot + cobweb
|
|
cobweb(LogisticMap(), 3.449,0.1,200,skip_transient = True)
|
|
liapunov(LogisticMap(),3.8) # > 0 so aperiodic or chaotic
0.22631731805389263 0.22631731805389263 |
def bifurcation_plot(f,k_min = 3.0, k_max = 4.0):
dk = (k_max - k_min)/1000.0
xpts = []
x = 0.5 #random()
for k in srange(k_min,k_max,dk):
# skip transient
x = orbit(f,k,200,x)
# plot trajectories for this k
ks = trajectory(f,k,100,x)
xpts = xpts + [[k,q] for q in ks]
return points(xpts, size = 1)
|
|
bifurcation_plot(DoubleUndoubleMap(),-4.0,-0.8)
|
|
# here we make a simple sage Sage "app"
@interact
def chaosxplorer(k = (0.0,4.0,.05),map =["Logistic","Square","Sine","UnOne"],\
start_value = (0,1,0.1), iteration = (1,100,1),skip_transient=False):
if map:
tom = eval(map + 'Map()')
cw = cobweb(tom, k,start_value,iterations = iteration,skip_transient=skip_transient)
bp = bifurcation_plot(tom,k-1,k)
g = graphics_array([cw, bp])
show(g, axes=True, figsize=(9,4),gridlines=False)
html("A %s map with k=%.1f"%(map, k))
Click to the left again to hide and once more to show the dynamic interactive window |
# leave and add fp's
#cobweb(0*x*(1-x),0) + cobweb(1*x*(1-x),0) + cobweb(2*x*(1-x),0) + cobweb(3*x*(1-x),0) + cobweb(4*x*(1-x),0)
|
|
bp1 = bifurcation_plot(LogisticMap(),2.9, 4.0)
bp1.show(axes=True, figsize=(9,9))
|
|
bpa = animate([bifurcation_plot(LogisticMap(),r-0.005,r) for r in srange(3.0,4.0,0.05)], figsize=[4,4],xmin=1.0,xmax=4,ymin=0.0, ymax=1.0)
bpa.show(delay=100)
|
bp3p = animate([bifurcation_plot(LogisticMap(),r-0.001,r) for r in srange(3.8,3.9,0.002)], figsize=[4,4],xmin=3.8,xmax=4.,ymin=0.0, ymax=1.0)
bp3p.show(delay=300)
|
# keep
#a = animate([cobweb(LogisticMap(),r*x*(1-x),0) for r in srange(0,4.,0.1)],xmin=0, xmax=1, figsize=[3,3])
|
|
#a.show()
|
|
# based on a blog post by william stein
# staffan liljegren added roessler and made it interactive
# next step is to make this interactive and add some more common attractors
@interact
def chaosxplorer2d(k = (0.0,4,.05),flow=["lorentz","roessler","atmosphere","ocean"], start_value = (0,1,0.1), iteration = (1,100,1)):
def lorenz(t,y,params):
return [params[0]*(y[1]-y[0]),y[0]*(params[1]-y[2])- y[1],y[0]*y[1]-params[2]*y[2]]
def lorenz_jac(t,y,params):
return [ [-param[0],params[0],0],[(params[1]-y[2]),-1,-y[0]],[y[1],y[0],-params[2]],[0,0,0]]
def roessler(t,y,params):
return [-(y[2] + y[1]),y[0] + params[0]*y[1],params[1]+ (y[0]-params[2])*y[2]]
def roessler_jac(t,y,params):
return None
roessler_par = [0.15,0.2,10.0]
T=ode_solver()
#T.algorithm="bsimp"
T.function=roessler
#T.jacobian=lorenz_jac
T.ode_solve(y_0=[.5,.5,.5],t_span=[0,155],params=roessler_par,num_points=20000)
l=[T.solution[i][1] for i in range(len(T.solution))]
line3d(l,thickness=0.3, figsize=5).show()
Click to the left again to hide and once more to show the dynamic interactive window |
hj = jacobian(henon, (x,y))
Traceback (click to the left of this block for traceback) ... in Python, and has the wrong precedence. Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "_sage_input_27.py", line 10, in <module>
exec compile(u'print _support_.syseval(python, u"var(\'x,y,a,b\')\\nhenon= matrix([[1+y-a*x^2],[b*x]]) \\nhj = jacobian(henon, (x,y))", __SAGE_TMP_DIR__)
File "", line 1, in <module>
File "/sagenb/sage_install/sage-4.7/devel/rkirov-flask/sagenb/misc/support.py", line 473, in syseval
return system.eval(cmd, sage_globals, locals = sage_globals)
File "/sagenb/sage_install/sage-4.7/local/lib/python2.6/site-packages/sage/misc/python.py", line 53, in eval
eval(compile(s, '', 'exec'), globals, globals)
File "", line 2, in <module>
File "element.pyx", line 699, in sage.structure.element.Element.__xor__ (sage/structure/element.c:5309)
RuntimeError: Use ** for exponentiation, not '^', which means xor
in Python, and has the wrong precedence.
|
det(hj)
Traceback (click to the left of this block for traceback) ... NameError: name 'hj' is not defined Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "_sage_input_26.py", line 10, in <module>
exec compile(u"print _support_.syseval(python, u'det(hj)', __SAGE_TMP_DIR__)" + '\n', '', 'single')
File "", line 1, in <module>
File "/sagenb/sage_install/sage-4.7/devel/rkirov-flask/sagenb/misc/support.py", line 473, in syseval
return system.eval(cmd, sage_globals, locals = sage_globals)
File "/sagenb/sage_install/sage-4.7/local/lib/python2.6/site-packages/sage/misc/python.py", line 56, in eval
eval(z, globals)
File "", line 1, in <module>
NameError: name 'hj' is not defined
|
|
|