# logLyapunov.py

'''
# https://en.wikipedia.org/wiki/Logistic_map
# https://en.wikipedia.org/wiki/Lyapunov_exponent
# https://ipython-books.github.io/121-plotting-the-bifurcation-diagram-of-a-chaotic-dynamical-system/

# Nonlinear Dynamics and Chaos, Ch 10.5; Example 10.5.3

The main significance of this figure is that one can easily distinguish the regions which are chaotic lambda_L > 0 from the regions which tend to a fixed point or limit cycle (DL < 0). 

You see several points (the first is at D = 0.75) where the Lyapuov exponent hits 0 and then goes negative again. These are the period doubling bifurcations. Precisely at the period doubling point the system is at the limit of chaos, but then becomes non-chaotic when the period doubles.  However, at the end of the period doubling regime, at D about 0.8922, DL crosses the axis and the system enters a chaotic regime

'''

import math
import matplotlib.pyplot as plt
from frange import *


NUM_RS = 10000
NUM_ITERS = 500

def logistic(r, x):
	return r * x * (1 - x)

def logistics(r):
  xs = [0.3]
  lyapunov = 0
  for i in range(NUM_ITERS-1):
    xs.append( logistic(r,xs[i]))
    lyapunov += math.log(abs(r - 2*r*xs[-1]))
        # sum of log( | f'(x) | )
  return xs[NUM_ITERS-100:], lyapunov/NUM_ITERS
     #  last 100 iterations


rs = linspace(2.0, 4.0, NUM_RS)  # R between 2 and 4
xVals = []
lyapunovs = []
rVals = []
for r in rs:
  xs, lya = logistics(r)
  xVals.append(xs)
  lyapunovs.append(lya)
  rVals.append([r] * 100) 

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 9),
                 sharex=True)

# display the bifurcation diagram
ax1.scatter(rVals, xVals, s=0.1)
ax1.set_xlabel('R')
ax1.set_ylabel('x')
ax1.set_xlim(2.5, 4)
ax1.set_title("Bifurcation Diagram")

# display the Lyapunov exponent
ax2.plot(rs, lyapunovs, '.k', ms=1)
ax2.axhline(0, color='k', lw=.5, alpha=.5)
ax2.set_xlim(2.5, 4)
ax2.set_ylim(-2, 1)
ax2.set_title("Lyapunov Exponent")
plt.tight_layout()

plt.show()