本文代码均以Python实现
首先给出按照定义法进行求解的程序，这里选取Henon映射为示例
#-*-coding:utf-8-*-
'''
多变量非线性方程求解
'''
from sympy import *
import numpy as np
np.set_printoptions(suppress=True)
n = 20000#控制迭代次数
def Henon(x,y,n):
for i in range(n):
x1 = 1 - 1.4 * x ** 2 + y
y1 = 0.3 * x
x = x1
y = y1
return x,y
def Jacobian():
count=0
a = 0.123456789
b = 0.123456789
# 使用符号方式求解
x, y = symbols("x,y")
f_mat = Matrix([1 - 1.4 * x ** 2 + y, 0.3 * x])
# 求解雅各比矩阵
jacobi_mat = f_mat.jacobian([x, y])#带变量的雅各比矩阵形式是固定的
a,b=Henon(a,b,5001)#先迭代1000次，消除初始影响.以第1001次的值作为初始值
# 这里为获取初始雅各比矩阵，将第一次放置在循环外
result = jacobi_mat.subs({x: a, y: b}) # 将变量替换为当前迭代值，得到当前的雅各比矩阵（数字）
J = result # 得到初始的雅各比矩阵
a, b = Henon(a, b, 1) # 每次迭代一次获取当前的迭代值
while(count<n-1):
result = jacobi_mat.subs({x: a, y: b}) # 将变量替换为当前迭代值，得到当前的雅各比矩阵(数字)
J = J*result # 计算累积的雅各比矩阵
a, b = Henon(a, b, 1) # 每次迭代一次获取当前的迭代值
count=count+1
return J
def LE_calculate(J):
eig_dic = J.eigenvals()#传入一个累积的雅各比矩阵
eig_list = list(eig_dic)#求累积雅各比矩阵的特征值
eig_1 = eig_list[0]
eig_2 = eig_list[1]
LE1=N(ln(abs(eig_1))/n)
LE2=N(ln(abs(eig_2))/n)
print(LE1)
print(LE2)
if __name__ == '__main__':
J=Jacobian()
LE_calculate(J)
值得注意的是，按照上面这种做法得到的Lyapunov指数通常
不可靠
，因为在计算Jacobian矩阵时，会遭遇
病态数值
的问题。在参考文献【1】【2】中已经简单论证了这一点。
于是，下面给出按照Wolf在1985年论文【3】中提出的
正交法
来计算Lyapunov指数的程序
#-*-coding:utf-8-*-
'''
多变量非线性方程求解
'''
from sympy import *
import numpy as np
# np.set_printoptions(suppress=True)
import copy
n = 5000#控制迭代次数
def Henon(x,y,n):
for i in range(n):
x1 = 1 - 1.4 * x ** 2 + y
y1 = 0.3 * x
x = x1
y = y1
return x,y
def LE_calculate():
sum_Lambda1 = 0
sum_Lambda2 = 0
a = 0.123456789
b = 0.123456789
# 使用符号方式求解
x, y = symbols("x,y")
f_mat = Matrix([1 - 1.4 * x ** 2 + y, 0.3 * x])
# 求解雅各比矩阵
jacobi_mat = f_mat.jacobian([x, y]) # 带变量的雅各比矩阵形式是固定的
a, b = Henon(a, b, 5001) # 先迭代5000次，消除初始影响.以第5001次的值作为初始值
U1 = Matrix([1, 0]) # 初始列向量
U2 = Matrix([0, 1])
for i in range(n):
J = jacobi_mat.subs({x: a, y: b}) # 将变量替换为当前迭代值，得到当前的雅各比矩阵（数字）
column_vector1 = U1#初始列向量为上一次的U1和U2
column_vector2 = U2
vector1 = J * column_vector1 # 初始列向量乘上雅各比矩阵之后得到的向量
vector2 = J * column_vector2
V1 = vector1 # 将vector1复制给V1
U1 = V1 / (V1.norm(2)) # 向量U1等于向量V1除以向量V1的模(2范数)
V2 = vector2 - (vector2.dot(U1)) * U1 # dot为点乘(內积)
U2 = V2 / (V2.norm(2))
Lambda1 = ln(V1.norm(2))
Lambda2 = ln(V2.norm(2))
sum_Lambda1 = sum_Lambda1 + Lambda1
sum_Lambda2 = sum_Lambda2 + Lambda2
a, b = Henon(a,b,1)#进行下一次迭代
LE1=sum_Lambda1/n
LE2=sum_Lambda2/n
print(LE1)
print(LE2)
if __name__ == '__main__':
LE_calculate()
这种方法同样也适用于计算1维离散混沌映射的Lyaounov指数。
下面给出采用正交法计算Logistic映射Lyaouov指数的代码
#-*-coding:utf-8-*-
'''
多变量非线性方程求解
'''
from sympy import *
import numpy as np
# np.set_printoptions(suppress=True)
import copy
n = 5000#控制迭代次数
def Logistic(x,n):
for i in range(n):
y = 4 * x * (1 - x)
x = y
return x
def LE_calculate():
sum_Lambda1 = 0
a = 0.123456789
# 使用符号方式求解
x = symbols("x")
f_mat = Matrix([4 * x * (1 - x)])
# 求解雅各比矩阵
jacobi_mat = f_mat.jacobian([x]) # 带变量的雅各比矩阵形式是固定的
a = Logistic(a, 5001) # 先迭代5000次，消除初始影响.以第5001次的值作为初始值
U1 = Matrix([1]) # 初始列向量
for i in range(n):
J = jacobi_mat.subs(x,a) # 将变量替换为当前迭代值，得到当前的雅各比矩阵（数字）
column_vector1 = U1#初始列向量为上一次的U1
vector1 = J * column_vector1 # 初始列向量乘上雅各比矩阵之后得到的向量
V1 = vector1 # 将vector1复制给V1
U1 = V1 / (V1.norm(2)) # 向量U1等于向量V1除以向量V1的模(2范数)
Lambda1 = ln(V1.norm(2))
sum_Lambda1 = sum_Lambda1 + Lambda1
a= Logistic(a,1)#进行下一次迭代
LE1=sum_Lambda1/n
print(LE1)
if __name__ == '__main__':
LE_calculate()
二维Arnold（猫映射）Lyapunov的计算一（正交法）
#-*-coding:utf-8-*-
'''
多变量非线性方程求解
'''
from sympy import *
import numpy as np
# np.set_printoptions(suppress=True)
import copy
n = 5000#控制迭代次数
def Arnold(x,y,n):
for i in range(n):
x1 = (x+y) % 1
y1 = (x+2*y) %1
x = x1
y = y1
# print(x)
# print(y)
return x,y
def LE_calculate():
sum_Lambda1 = 0
sum_Lambda2 = 0
a = 0.1
b = 0.1
# 使用符号方式求解
x, y = symbols("x,y")
f_mat = Matrix([x+y, x+2*y])
# 求解雅各比矩阵
J = f_mat.jacobian([x, y]) # 带变量的雅各比矩阵形式是固定的
# print(J)
a, b = Arnold(a, b, 5001) # 先迭代5000次，消除初始影响.以第5001次的值作为初始值
U1 = Matrix([1, 0]) # 初始列向量
U2 = Matrix([0, 1])
for i in range(n):
column_vector1 = U1#初始列向量为上一次的U1和U2
column_vector2 = U2
vector1 = J * column_vector1 # 初始列向量乘上雅各比矩阵之后得到的向量
vector2 = J * column_vector2
#这里值得注意的是，需要调用N将表达式计算成数字。否则会因为表达式太复杂，使得程序无法计算（在迭代10次就会出现）
V1 = vector1 # 将vector1复制给V1
U1 = N(V1 / (V1.norm(2))) # 向量U1等于向量V1除以向量V1的模(2范数)
V2 = N(vector2 - (vector2.dot(U1)) * U1) # dot为点乘(內积)
U2 = N(V2 / (V2.norm(2)))
Lambda1 = ln(V1.norm(2))
Lambda2 = ln(V2.norm(2))
sum_Lambda1 = sum_Lambda1 + Lambda1
sum_Lambda2 = sum_Lambda2 + Lambda2
a, b = Arnold(a,b,1)#进行下一次迭代
print(i)
LE1=N(sum_Lambda1/n)
LE2=N(sum_Lambda2/n)
print(LE1)
print(LE2)
if __name__ == '__main__':
LE_calculate()
二维猫映射的LE（特征值法），它的两个LE实际上就是猫映射的Jacobian矩阵的两个特征值取对数。这是因为该映射的Jacobian矩阵是一个常数矩阵
#-*-coding:utf-8-*-#-*-coding:utf-8-*-
'''
多变量非线性方程求解
'''
from sympy import *
import numpy as np
np.set_printoptions(suppress=True)
n = 5000#控制迭代次数
def Arnold(x,y,n):
for i in range(n):
x1 = (x+y)%1
y1 = (x+2*y)%1
x = x1
y = y1
return x,y
def Jacobian():
# 使用符号方式求解
x, y = symbols("x,y")
f_mat = Matrix([x+y, x+2*y])
# 求解雅各比矩阵
result = f_mat.jacobian([x, y])#带变量的雅各比矩阵形式是固定的
J=result
return J
def LE_calculate(J):
eig_dic = J.eigenvals()#传入一个累积的雅各比矩阵
eig_list = list(eig_dic)#求累积雅各比矩阵的特征值
eig_1 = eig_list[0]
eig_2 = eig_list[1]
LE1=N(ln(eig_1))
LE2=N(ln(eig_2))
print(LE1)
print(LE2)
if __name__ == '__main__':
J=Jacobian()
LE_calculate(J)
总结：计算Lyapunov指数有很多种方法，不同的方法计算出来的Lyapunov指数会有差别，但差别不会很大。这是因为这些方法的思想都是线性逼近。这种差别一般控制在0.01~0.1之间。
当差别很大时，考虑两个问题
（1）计算方法错误
（2）迭代次数不够
参考文献：
【1】黄润生
【2】非光滑动力系统Lyapunov指数谱的计算方法-金俐，陆启超
【3】Determing lyapunov exponents from a time series-Alan Wolf
版权声明：本文为u014593570原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。