N维离散混沌映射Lyapunov指数的计算

本文代码均以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