### 李雅普诺夫指数的Python实现 李雅普诺夫指数是衡量动力系统稳定性的关键指标之一。为了计算该指数，可以通过数值方法估计两个初始接近轨迹随时间分离的速度。以下是基于参考文献[^1]和[^3]中的理论框架设计的一个简单Python实现。 #### 实现思路 1. 定义系统的动态方程 $f(x)$。 2. 给定初值点 $x_0$ 和其邻近点 $x_1$。 3. 迭代更新这两个点的位置并记录它们的距离变化。 4. 使用自然对数函数 $\ln(\cdot)$ 对距离增长率取平均值得到最终的李雅普诺夫指数。 下面是具体的代码示例： ```python import numpy as np def lyapunov_exponent(f, df_dx, x0, dt=0.01, T=100, delta=1e-8): """ 计算离散动力系统的最大李雅普诺夫指数 参数: f (function): 动力系统的映射函数 f(x) df_dx (function): 映射函数关于 x 的导数 df/dx x0 (float): 初始条件 dt (float): 时间步长，默认为 0.01 T (int): 总迭代次数，默认为 100 delta (float): 邻近距离的初始扰动量，默认为 1e-8 返回: float: 最大李雅普诺夫指数 λ """ sum_log = 0 x = x0 dx = delta for _ in range(T): # 更新当前状态 x_new = f(x) # 更新微小扰动的增长情况 dx = df_dx(x) * dx # 归一化扰动向量 norm_factor = abs(dx) dx /= norm_factor # 累加 log 值 sum_log += np.log(norm_factor) # 准备下一次迭代 x = x_new # 计算李雅普诺夫指数 lambda_exp = sum_log / (T * dt) return lambda_exp # 测试洛伦兹吸引子的最大李雅普诺夫指数 if __name__ == "__main__": sigma = 10 rho = 28 beta = 8/3 def lorenz_x_dot(state): """Lorenz 方程的第一个分量""" x, y, z = state return sigma * (y - x) def jacobian_lx(state): """Jacobian 关于第一个变量的偏导数""" _, y, _ = state return -sigma + sigma * (y / abs(y)) if y != 0 else -sigma initial_state = [1.0, 1.01, 0] result = lyapunov_exponent( f=lambda s: lorenz_x_dot(s), df_dx=jacobian_lx, x0=np.array(initial_state[:2]), dt=0.01, T=10000, delta=1e-5 ) print(f"Estimated Largest Lyapunov Exponent: {result}") ``` 上述代码实现了针对连续动力系统的最大李雅普诺夫指数估算功能，并以经典的洛伦兹吸引子为例进行了测试。注意此程序仅适用于简单的低维系统；对于高维混沌系统可能需要更复杂的处理方式。 --- ###