原始论文

M.T. Rosenstein, J.J. Collins, and C.J. De Luca. A practical method for calculating largest Lyapunov exponents from small data sets. Physica D, 65:117-134, 1993.

下载地址：https://www.physionet.org/content/lyapunov/1.0.0/

python 相关代码

NonLinear Time Series Analysis（nolitsa）
NOnLinear measures for Dynamical Systems (nolds)

混沌系统的常见指标

区分确定性混沌系统与噪声已成为许多不同领域的重要问题。

对于实验产生的时间序列，可以计算这些混沌系统的指标：

相关维数（D2D_2D2），
Kolmogorov 熵
Lyapunov 特征指数。

相关维度是对系统复杂程度的估计，熵和特征指数是对混沌程度的估计。

最大李亚普诺夫指数的含义

LLE 描述了相空间中相近的两点（初始间距为CCC）随时间推移指数分离的速率：
d(t)=Ceλ1td(t) = Ce^{\lambda_1 t} d(t)=Ceλ1t其中d(t)d(t)d(t)表示分离距离，CCC表示初始间距，λ1\lambda_1λ1 为最大李氏指数。

算法流程图

python 代码模块

最近邻

import numpy as npfrom scipy import stats
from scipy.spatial import cKDTree as KDTree
from scipy.spatial import distancedef neighbors(y, metric='chebyshev', window=0, maxnum=None):"""Find nearest neighbors of all points in the given array.Finds the nearest neighbors of all points in the given array usingSciPy's KDTree search.Parameters----------y : ndarrayN-dimensional array containing time-delayed vectors.metric : string, optional (default = 'chebyshev')Metric to use for distance computation.  Must be one of"cityblock" (aka the Manhattan metric), "chebyshev" (aka themaximum norm metric), or "euclidean".window : int, optional (default = 0)Minimum temporal separation (Theiler window) that should existbetween near neighbors.  This is crucial while computingLyapunov exponents and the correlation dimension.maxnum : int, optional (default = None (optimum))Maximum number of near neighbors that should be found for eachpoint.  In rare cases, when there are no neighbors that are at anonzero distance, this will have to be increased (i.e., beyond2 * window + 3).Returns-------index : arrayArray containing indices of near neighbors.dist : arrayArray containing near neighbor distances."""if metric == 'cityblock':p = 1elif metric == 'euclidean':p = 2elif metric == 'chebyshev':p = np.infelse:raise ValueError('Unknown metric.  Should be one of "cityblock", ''"euclidean", or "chebyshev".')tree = KDTree(y)n = len(y)if not maxnum:maxnum = (window + 1) + 1 + (window + 1)else:maxnum = max(1, maxnum)if maxnum >= n:raise ValueError('maxnum is bigger than array length.')dists = np.empty(n)indices = np.empty(n, dtype=int)for i, x in enumerate(y):for k in range(2, maxnum + 2):dist, index = tree.query(x, k=k, p=p)valid = (np.abs(index - i) > window) & (dist > 0)if np.count_nonzero(valid):dists[i] = dist[valid][0]indices[i] = index[valid][0]breakif k == (maxnum + 1):raise Exception('Could not find any near neighbor with a ''nonzero distance.  Try increasing the ''value of maxnum.')return np.squeeze(indices), np.squeeze(dists)

maximum Lyapunov exponent

def mle(y, maxt=500, window=10, metric='euclidean', maxnum=None):"""Estimate the maximum Lyapunov exponent.Estimates the maximum Lyapunov exponent (MLE) from amulti-dimensional series using the algorithm described byRosenstein et al. (1993).Parameters----------y : ndarrayMulti-dimensional real input array containing points in thephase space.maxt : int, optional (default = 500)Maximum time (iterations) up to which the average divergenceshould be computed.window : int, optional (default = 10)Minimum temporal separation (Theiler window) that should existbetween near neighbors (see Notes).maxnum : int, optional (default = None (optimum))Maximum number of near neighbors that should be found for eachpoint.  In rare cases, when there are no neighbors that are at anonzero distance, this will have to be increased (i.e., beyond2 * window + 3).Returns-------d : arrayAverage divergence for each time up to maxt.Notes-----This function does not directly estimate the MLE.  The MLE should beestimated by linearly fitting the average divergence (i.e., theaverage of the logarithms of near-neighbor distances) with time.It is also important to choose an appropriate Theiler window so thatthe near neighbors do not lie on the same trajectory, in which casethe estimated MLE will always be close to zero."""index, dist = utils.neighbors(y, metric=metric, window=window,maxnum=maxnum)m = len(y)maxt = min(m - window - 1, maxt)d = np.empty(maxt)d[0] = np.mean(np.log(dist))for t in range(1, maxt):t1 = np.arange(t, m)t2 = index[:-t] + t# Sometimes the nearest point would be farther than (m - maxt)# in time.  Such trajectories needs to be omitted.valid = t2 < mt1, t2 = t1[valid], t2[valid]d[t] = np.mean(np.log(utils.dist(y[t1], y[t2], metric=metric)))return d

RANSAC 拟合曲线

需要先安装 sklearn 库

def poly_fit(x, y, degree, fit="RANSAC"):# check if we can use RANSACif fit == "RANSAC":try:# ignore ImportWarnings in sklearnwith warnings.catch_warnings():warnings.simplefilter("ignore", ImportWarning)import sklearn.linear_model as sklinimport sklearn.preprocessing as skpreexcept ImportError:warnings.warn("fitting mode 'RANSAC' requires the package sklearn, using"+ " 'poly' instead",RuntimeWarning)fit = "poly"if fit == "poly":return np.polyfit(x, y, degree)elif fit == "RANSAC":model = sklin.RANSACRegressor(sklin.LinearRegression(fit_intercept=False))xdat = np.asarray(x)if len(xdat.shape) == 1:# interpret 1d-array as list of len(x) samples instead of# one sample of length len(x)xdat = xdat.reshape(-1, 1)polydat = skpre.PolynomialFeatures(degree).fit_transform(xdat)try:model.fit(polydat, y)coef = model.estimator_.coef_[::-1]except ValueError:warnings.warn("RANSAC did not reach consensus, "+ "using numpy's polyfit",RuntimeWarning)coef = np.polyfit(x, y, degree)return coefelse:raise ValueError("invalid fitting mode ({})".format(fit))

例子：计算洛伦兹系统的最大李雅普诺夫指数

import warnings
from nolitsa import data, lyapunov
import numpy as np
import matplotlib.pyplot as pltdt = 0.01
x0 = [0.62225717, -0.08232857, 30.60845379]
x = data.lorenz(length=4000, sample=dt, x0=x0,sigma=16.0, beta=4.0, rho=45.92)[1]
plt.plot(range(len(x)),x)
plt.show()# Choose appropriate Theiler window.
meanperiod = 30
maxt = 250
d = lyapunov.mle(x, maxt=maxt, window=meanperiod)
t = np.arange(maxt) *dt
coefs = poly_fit(t, d, 1)
print('LLE = ', coefs[0])plt.title('Maximum Lyapunov exponent for the Lorenz system')
plt.xlabel(r'Time $t$')
plt.ylabel(r'Average divergence $\langle d_i(t) \rangle$')
plt.plot(t, d, label='divergence')
plt.plot(t, t * 1.50, '--', label='slope=1.5')
plt.plot(t, coefs[1] +coefs[0]* t, '--', label='RANSAC')
plt.legend()
plt.show()

求最大李雅普诺夫指数（Largest Lyapunov Exponents，LLE）的 Rosenstein 算法相关推荐

Matlab求解混沌系统最大李雅普诺夫指数
李雅普诺夫指数(Lyapunov)是一个较为典型的判断一个系统是否具有混沌特性以及混沌的程度分析方法. 李雅普诺夫指数:在相空间中初始时无限接近的两个轨道,随着时间的不断推移按指数收敛或发散的平均变化 ...
matlab混沌指数的计算,matlab实现混沌系统最大李雅普诺夫指数
李雅普诺夫指数(Lyapunov)是一个较为典型的判断一个系统是否具有混沌特性以及混沌的程度分析方法. 李指数:在相空间中初始时无限接近的两个轨道,随着时间的不断推移按指数收敛或发散的平均变化率,它可 ...
Matlab编程之混沌系统李雅普诺夫指数分析
简介李雅普诺夫指数是衡量混沌系统的一个重要参数,下面截图是对其具体解释. 代码实现: clc; clear; global kk;e=0ina1=0;final2=10;for kk=ina1:1: ...
李雅普诺夫指数对应的特征方向
洛伦兹系统的李雅普诺夫指数的符号为 (+,0,−)(+,0,-)(+,0,−) Any continuous time-dependent dynamical system without a fix ...
C语言使用1到9求出所有k个数字的所有组合的算法(附完整源码)
C语言使用1到9求出所有k个数字的所有组合的算法 C语言使用1到9求出所有k个数字的所有组合的算法完整源码(定义,实现,main函数测试) C语言使用1到9求出所有k个数字的所有组合的算法完整源码(定 ...
【机器人】基于指数积的机械臂正运动学算法
基于指数积的机械臂正运动学算法 1.前言 2.指数积公式建立过程 3.PoE实例 4.PoE与DH对比 1.前言在学习了刚体运动的指数坐标表示和运动旋量后,我又对使用指数积法(PoE)对机械臂进行正 ...
分段线性映射PWLCM的李雅普诺夫指数Lyapunov的matlab实现
clc;clear; p1=0:0.0001:0.5; p2=0.5:0.0001:1; p=0:0.0001:1;for i=1:length(p1)y1(i)=log(abs(1/p1(i))); ...
matlab求莫兰指数程序,python计算莫兰指数(Moran's I)并绘制地区热力图——以中国各省pm2.5为例...
[TOC] 程序简介计算省的pm2.5平均值作为观测矩阵,省会的距离的倒数作为空间权重矩阵,计算全局莫兰指数为0.49,显著性检验p值为3.75>1.96,得出中国地区的pm2.5存在空间正相 ...
统计学怎么求加权指数_指数统计学试卷
统计指数一.单项选择题 1. 按照个体价格指数和报告期销售额计算的价格总指数是( ) A. 综合指数 B. 平均指标指数 C. 加权算术平均数指数 D. 加权调和平均数指数 2. 在由3个指数所组成 ...
学习记录563@求模下指数幂的快速算法（求模指幂快速算法）
令a,x,n 为正整数且 a < n.公钥密码体系常需要求模下指数署 axa^xax mod n,如果先求y=axa^xax再求y mod n则所需时间太多,y也太大,因为axa^xax mod ...

求最大李雅普诺夫指数（Largest Lyapunov Exponents，LLE）的 Rosenstein 算法

文章目录