# LogisticMapLCE.py:
#	Estimate the Lyapunov Characteristic Exponent for the Logistic Map.
#   Plot it as a function of the map's control parameter.

# Import modules
# SciPy for the arange function
from scipy import arange
# Plotting: Use this for Linux/Unix/OS X
from pylab import *
# OR
# Plotting: Use this for Windows
#from matplotlib.matlab import *

# Define the Logistic map's function, r in [0,4]
def LogisticMap(r,x):
	return r * x * (1.0 - x)

# Define the Logistic map's derivative
def LogisticMapDeriv(r,x):
	return r * (1.0 - 2.0 * x)

# Setup parameter range
plow  = 2.95
phigh = 4.0

# Setup the plot
# It's numbered 1 and is 8 x 6 inches.
figure(1,(8,6))
# Stuff parameter range into a string via the string formating commands.
TitleString = 'Logistic map Lyapunov Characteristic Exponent for r in [%g,%g]' % (plow,phigh)
title(TitleString)
# Label axes
xlabel('Control parameter r')
ylabel('$\lambda$')

# Set the initial condition used at each parameter value
ic = 0.2
# Establish the arrays to hold the LCE estimates at each parameter value
psweep = [ ] # The parameter values
lce    = [ ] # The LCE values
# The iterates we throw away
nTransients = 100
# The iterations over which we estimate the LCE (more is better)
nIterates = 2500
# This sets how dense the parameter sweep will be
nSteps = 100.0
# Sweep the control parameter over the desired range
for p in arange(plow,phigh,(phigh-plow)/nSteps):
	# Set the initial condition to the reference value
	state = ic
	# Throw away the transient iterations
	for i in range(nTransients):
		state = LogisticMap(p,state)
	# Now estimate the LCE over nIterates
	# We start an estimate at a new parameter, must clear previous LCE
	l = 0.0
	for i in range(nIterates):
		state = LogisticMap(p,state)
		l = l + log(abs(LogisticMapDeriv(p,state)))
	# Convert to the per-iteration average
	l = l / float(nIterates)
	lce.append(l)
	psweep.append(p)

# Plot the list of (p,lce) pairs as lines
plot(psweep, lce, 'r-')
# Since the boundary between stability and chaos occurs at LCE = 0, we draw
axhline(0.0,0.0,1.0,color = 'k')

# Use this to save figure as a bitmap png file
#savefig('LogisticMapLCE', dpi=600)

# Display plot in window
show()