使用Python绘制一个非常漂亮的李雅普诺夫指数，及其函数值

import numpy as np
import matplotlib.pyplot as plt
# show plots in notebook
# % matplotlib inline

result = []
lambdas = []
maps = []

# define range of r
rvalues = np.arange(0, 2, 0.01)

# loop through r
for r in rvalues:
    x = 0.1
    result = []
    # iterate system 100 times
    for t in range(100):
        x = x + r - x**2
        # calculate log of the absolute of the derivative
        # 计算李雅普诺夫指数最关键的地方是，计算导数；
        result.append(np.log(abs(1 - 2*x)))
    # take average
    lambdas.append(np.mean(result))
    # for the map ignore first 100 iterations as transient time and iterate anew
    for t in range(20):
        x = x + r - x**2
        maps.append(x)

fig = plt.figure(figsize=(10,7))
ax1 = fig.add_subplot(1,1,1)

xticks = np.linspace(0, 2, 4000)
# zero line
zero = [0]*4000
ax1.plot(xticks, zero, 'g-')
# plot map
ax1.plot(xticks, maps, 'r.',alpha = 0.3, label = 'Map')
ax1.set_xlabel('r')
# plot lyapunov
ax1.plot(rvalues, lambdas, 'b-', linewidth = 3, label = 'Lyapunov exponent')
ax1.grid('on')
ax1.set_xlabel('r')
ax1.legend(loc='best')
ax1.set_title('Map of x(t+1) = x(t) + r - x(t)^2 versus Lyapunov exponent')
plt.show()

代码解析：

1、xticks的长度是4000,200*20；

2、一共for循环了200个不同的r值（即模型参数）；

3、每个r值计算出来一个李雅普诺夫指数；

4、每个r值计算出来20个不同的幅值，（如果振动比较小，这个20个值，很接近），如果进入混沌态，这20个值差异就很大，这说明此时整个系统对模型参数的微小变化非常敏感；

5、