Exploring Iterated Maps

Python 2.3.3 (#1, Aug 23 2004, 20:06:57)
[GCC 3.3 20030304 (Apple Computer, Inc. build 1640)] on darwin
Type "copyright", "credits" or "license()" for more information.

****************************************************************
Personal firewall software may warn about the connection IDLE
makes to its subprocess using this computer's internal loopback
interface. This connection is not visible on any external
interface and no data is sent to or received from the Internet.
****************************************************************

IDLE 1.0.2

Set up the map explorer for the logistic map

Set the parameter to .5 and look at the function and an orbit

>>> m.param = .5
>>> w1 = m.plotComposition(1)
>>> w2 = m.plotOrbit(.5,20)
>>> w1.delete()
>>> w2.delete()

>>> m.param = 1.5
>>> w1 = m.plotComposition(1)
>>> w2 = m.plotOrbit(.001,20)
>>> w2.delete()
>>> w2 = m.plotOrbit(.001,40)
>>> w1.delete()
>>> w2.delete()

Increase the parameter further. This time look at the second composition, and look for a stable period-2 orbit. Look at an orbit. Is the fixed point still stable

>>> m.param = 2.5
>>> w1 = m.plotComposition(1)
>>> w2 = m.plotComposition(2)
>>> w3 = m.plotOrbit(.001,40)
>>> for w in [w1,w2,w3]: w.delete()

Now increase the parameter again. Is there a stable period-2 orbit now? Where does the orbit go if we start it off near the unstable fixed point?

>>> m.param = 3.2
>>> w1 = m.plotComposition(1)
>>> w2 = m.plotComposition(2)
>>> xFixed = (3.2-1)/3.2
>>> xFixed
0.6875
>>> w3 = m.plotOrbit(xFixed + .01,50)
>>> for w in [w1,w2,w3]: w.delete()

The map explorer class has a method that finds periodic orbits and their stability. Here we look at period 2 orbits.

>>> orbits,stability = m.findPeriodicOrbits(2)
>>> orbits
[0.0, 0.51251251251251251, 0.68668668668668664, 0.79879879879879878]
>>> stability
[1.1631508098056809, -0.93703717919219898, 0.18273824089059557, -0.87138324820310187]

The stability parameter given above is the exponential growth rate of perturbations to the orbit. If the value is positive, the orbit is unstable. If the value is negative, the orbit is stable. Note that the list of orbits contains the two unstable fixed points as well as the two points making up the stable period-2 orbit.

Now make an orbit diagram showing the periodic orbits up to period 10, as the parameter of the map is varied over the interval [3,4].

>>> from orbitDiagram import *
>>> w = orbitDiagram(map1D.logistic,3.,4.,100)
Working on 3.0
Working on 3.01
Working on 3.02
Working on 3.03
Working on 3.04
Working on 3.05
(... etc.)
Working on 3.7
Working on 3.71
Working on 3.72
Working on 3.99

Now pick a parameter value in the chaotic zone, and look at how two neighboring orbits behave.

>>> m.param = 3.8
>>> w1 = m.plotOrbit(.5,100)
>>> orbit1 = m.orbit(.5,100)
>>> orbit2 = m.orbit(.5001,100)
>>> w2 = ezplot(orbit1,orbit2)

>>> diff = [abs(orbit1[i]-orbit2[i]) for i in range(len(orbit1))]
>>> w3 = ezplot(diff)
>>> w3.delete()
>>> w3 = ezplot(diff[0:30])

Try an orbit that starts much closer to the first one

>>> orbit3 = m.orbit(.500001,100)
>>> diff = [abs(orbit1[i]-orbit3[i]) for i in range(len(orbit1))]
>>> w3.delete()
>>> w3 = ezplot(diff)
>>> w3.delete()
>>> w3 = ezplot(diff[0:55])

The Lyapunov exponent is the exponential growth rate of deviations between neighboring orbits. You can get it from the slope of the log of the deviation vs. iteration (time).

>>> import math
>>> logdiff = [math.log(d) for d in diff]
>>> w4 = ezplot(logdiff)

>>> dir(m)
['__doc__', '__init__', '__module__', 'compose', 'f', 'findPeriodicOrbits', 'fprime', 'lyap', 'orbit', 'param', 'plotComposition', 'plotOrbit']

The mapExplorer class has a method to compute this slope over a specified number of iterates. The Lyapunov exponent is defined as the limit of the growth rate for large numbers of iterates. Let's try it out.

Over 100 iterates, the estimate of the Lyapunov exponent seems quite dependant on where we start the orbit. The sensitivity becomes less if we take orbits of length 1000, since the orbit then has time to explore more of the space, and sample the full behavior.

>>> m.lyap(.5,1000)
0.41825336668042562
>>> m.lyap(.01,1000)
0.43848255428504934
>>>