Lyapunov exponents

系统动力学 Lyapunov 混沌
发布日期:   2022-11-11
更新日期:   2024-06-27
文章字数:   1.9k
阅读时长:   11 分
阅读次数:  
  • Lyapunov指数，混沌系统的基本特点就是系统对初始值的极端敏感性，两个相差无几的初值所产生的轨迹，随着时间的推移按指数方式分离。

Lyapunov exponents

最大李雅普诺夫指数

Lorenz系统族模型

model

方法

matlab

Z=[];  
a=10; 
c=-5.75; 
d0=1e-7; 
bs = linspace(0,300,301);
transient = 50;
for b=bs
    params = [a,b,c];
    lsum=0; 
    x=1;y=1;z=1;   % #初始基准点
    x1=1;y1=1;z1=1+d0;  % #初始偏离点
    for i=1:100  
        [T1,Y1]=ode45(@(t,X) Lorenz(t,X,params),[0 1],[x;y;z]); 
        [T2,Y2]=ode45(@(t,X) Lorenz(t,X,params),[0 1],[x1;y1;z1]); 
        n1=length(Y1);n2=length(Y2); 
        x=Y1(n1,1);y=Y1(n1,2);z=Y1(n1,3);   
        x1=Y2(n2,1);y1=Y2(n2,2);z1=Y2(n2,3);   
        d1=sqrt((x-x1)^2+(y-y1)^2+(z-z1)^2);     
        % #新的偏离点在上一次计算的两轨迹末端的连线上，且距离仍等于d0
        x1=x+(d0/d1)*(x1-x);   
        y1=y+(d0/d1)*(y1-y);    
        z1=z+(d0/d1)*(z1-z); 
        % #舍弃暂态过程的数据，因为初始基准点不一定在吸引子上
        if i> transient
            lsum=lsum+log(d1/d0);  
        end  
    end  
    Z=[Z lsum/(i-transient)]; 
end 

plot(bs,Z,'-k');  
title('')  
xlabel('parameter b'),ylabel('Largest Lyapunov Exponents'); 
grid on;

function dX = Lorenz(t,X,params) 

a = params(1);
b = params(2);
c = params(3);

x=X(1); 
y=X(2); 
z=X(3);

dX = zeros(3,1);
dX(1)=a*(y-x);
dX(2)=-c*y-x*z;
dX(3)=x*y-b;

end

image-20221111105451392

fortran

module Lorenz
    implicit none
    real,parameter :: h=0.001
    integer,parameter :: MaxT=2500000,N=3,M=1,T_trans=2000000
    real(kind=8) :: x1(M,N),x2(M,N)
    real :: tau(M),b
contains

    subroutine fnf(xx,fx,i)
        implicit none
        real(kind=8) :: xx(N),fx(N)
        integer :: t,i,j
        real :: a,c
        tau=1.0
        a=10.0
        c=-5.75
        fx(1)=tau(i)*a*(xx(2)-xx(1))
        fx(2)=tau(i)*(-xx(1)*xx(3)-c*xx(2))
        fx(3)=tau(i)*(xx(1)*xx(2)-b)
        return
    end subroutine fnf


    subroutine rk4(x)
        implicit none
        integer :: i
        real(kind=8) :: x(M,N),xx(N),fx(N)
        real(kind=8) :: kx1(N),kx2(N),kx3(N),kx4(N)
        do i=1,M,1
            xx=x(i,:)
            call fnf(xx,fx,i)
            kx1=h*fx
            xx=x(i,:)+0.5*kx1
            call fnf(xx,fx,i)
            kx2=h*fx
            xx=x(i,:)+0.5*kx2
            call fnf(xx,fx,i)
            kx3=h*fx
            xx=x(i,:)+kx3
            call fnf(xx,fx,i)
            kx4=h*fx
            x(i,:)=x(i,:)+(kx1+2.0*kx2+2.0*kx3+kx4)/6.0
        end do
        return
    end subroutine rk4


end module Lorenz



program main
    use Lorenz
    implicit none
    integer :: t,i
    real(kind=8) :: d1=0.0,d0,lam,lamda
    open(10,file="lamda.txt")
    d0=1.0e-7
    !初值基准点
    do i=0,300,1
        lam=0.0
        lamda=0.0
        b=real(i)
        x1(1,1)=1.0
        x1(1,2)=1.0
        x1(1,3)=1.0
        x2=x1
        x2(1,3)=x1(1,3)+d0
        do t=1,MaxT,1
            call rk4(x1)
            call rk4(x2)
            d1=sqrt((x2(1,1)-x1(1,1))**2+(x2(1,2)-x1(1,2))**2+(x2(1,3)-x1(1,3))**2)
            !新的偏离点在上一次计算的两轨迹末端的连线上，且距离仍等于d0
            x2(1,1)=x1(1,1)+(d0/d1)*(x2(1,1)-x1(1,1))
            x2(1,2)=x1(1,2)+(d0/d1)*(x2(1,2)-x1(1,2))
            x2(1,3)=x1(1,3)+(d0/d1)*(x2(1,3)-x1(1,3))
            if(t>T_trans) then
                lam=lam+log(d1/d0)
            end if
        end do
        lamda=lam/((MaxT-T_trans-1))/h
        write(10,*) b,lamda
!        write(*,*) b,lamda
    end do
    close(10)
end program main

结果

李雅普诺夫指数谱

Rossler系统模型

model

matlab

% 求解LE代码
% 计算Rossler吸引子的Lyapunov指数
clear;
yinit=[1,1,1];
orthmatrix=[1 0 0;
    0 1 0; 
    0 0 1];
a=0.15;
b=0.20;
c=10.0;
y=zeros(12,1);
% 初始化输入
y(1:3) = yinit;
y(4:12) = orthmatrix;
tstart = 0; % 时间初始值 
tstep = 1e-3; % 时间步长
wholetimes = 1e5; % 总的循环次数
steps = 10; %每次演化的步数
iteratetimes = wholetimes/steps; % 演化的次数
mod = zeros(3,1);
lp = zeros(3,1);
% 初始化三个Lyapunov指数
Lyapunov1 = zeros(iteratetimes,1);
Lyapunov2 = zeros(iteratetimes,1);
Lyapunov3 = zeros(iteratetimes,1);
for i=1:iteratetimes
    tspan = tstart:tstep:(tstart+tstep*steps);
    [T,Y] = ode45(@(tspan,y) Rossler_Ly(tspan,y),tspan,y);
    %取积分得到的最后一个时刻的值
    y=Y(size(Y,1),:);
    %重新定义起始时刻
    tstart = tstart + tstep*steps;
    y0=[y(4) y(7) y(10);
        y(5) y(8) y(11);
        y(6) y(9) y(12)];
    % 正交化
    y0=ThreeGS(y0);
    % 取三个向量的模
    mod(1) = sqrt(y0(:,1)'*y0(:,1));
    mod(2) = sqrt(y0(:,2)'*y0(:,2));
    mod(3) = sqrt(y0(:,3)'*y0(:,3));
    y0(:,1) = y0(:,1)/mod(1);
    y0(:,2) = y0(:,2)/mod(2);
    y0(:,3) = y0(:,3)/mod(3);
    lp = lp+log(abs(mod));
    % 三个Lyapunov指数
    Lyapunov1(i) = lp(1)/(tstart);
    Lyapunov2(i) = lp(2)/(tstart);
    Lyapunov3(i) = lp(3)/(tstart);
    y(4:12) = y0';
end
% 作Lyapunov指数谱图
i=1:iteratetimes;
plot(i,Lyapunov1,Lyapunov2,Lyapunov3);
Le1=Lyapunov1(end),Le2=Lyapunov2(end),Le3=Lyapunov3(end)
title('Lyapunov 指数谱图')
function dX = Rossler_Ly(t,X)
% Rossler 吸引子,用来计算Lyapunov指数
% a=0.15,b=0.20,c=10.0
% dx/dt = -y-z
% dy/dt = x+ay
% dz/dt = b+z(x-c)
a=0.15;
b=0.20;
c=10.0;
x=X(1);y=X(2);z=X(3);
%Y 的三个列向量为相互正交的单位向量
Y=[X(4),X(7),X(10);
    X(5),X(8),X(11);
    X(6),X(9),X(12)];
% 输出向量初始化，必不可少
dX=zeros(12,1);
% Rossler 吸引子
dX(1) = -y-z;
dX(2) = x+a*y;
dX(3) = z*(x-c)+b;
%Rossler 吸引子的Jaxobi矩阵
Jaco=[0 -1 -1;
    1 a 0;
    z 0 x-c];
dX(4:12) = Jaco*Y;
end
% G-S 正交化
function A = ThreeGS(V) %V为3*3向量
v1 = V(:,1);
v2 = V(:,2);
v3 = V(:,3);
a1 = zeros(3,1);
a2 = zeros(3,1);
a3 = zeros(3,1);
a1 = v1;
a2 = v2-((a1'*v2)/(a1'*a1))*a1;
a3 = v3-((a1'*v3)/(a1'*a1))*a1-((a2'*v3)/(a2'*a2))*a2;
A=[a1,a2,a3];
end

result

fortran

python

Lorenz系统模型

model

Lorenz模型(sun2009文章)

image-20221111144934935

Lorenz系统族模型

model

!------------Calculate Lyaponov Exponent Spectrum of Lorenze System----------------------
!-----------------------Lorenze Equation:-------------------
!------------------------dx/dt=aa*(y-x)----------------------
!------------------------dy/dt=cc*x-x*z-y--------------------
!------------------------dz/dt=x*y-bb*z----------------------
module head
    implicit none
    real(kind=8), parameter :: aa=10.0D0, bb=100.0D0, cc=-5.75D0, dt=0.001D0
    integer(kind=4), parameter :: st=1e5,n=1e6
    integer(kind=4), parameter :: m=3
contains

    subroutine initial(in)
        implicit none
        real(kind=8), dimension(m) :: in
        call random_seed()
        call random_number(in)
    end subroutine

    subroutine fnf1(in,out)
        implicit none
        real(kind=8), dimension(m) :: in, out
        out(1)=dt*(aa*(in(2)-in(1)))
        out(2)=dt*(-cc*in(2)-in(1)*in(3))
        out(3)=dt*(in(1)*in(2)-bb)
    end subroutine

    subroutine fnf2(in,out,jj)
        implicit none
        real(kind=8), dimension(m,m) :: in, out
        real(kind=8), dimension(m,m) :: jj
        integer(kind=4) :: i
        out=dt*matmul(jj,in) 
    end subroutine

end module

program main
    use head
    implicit none
    integer(kind=4) :: step,i,j,ii
    real(kind=8), dimension(m) :: fnorm, flya, fnorm1,fnorm2
    real(kind=8), dimension(m) :: y,k1,k2,k3,k4
    real(kind=8), dimension(m,m) :: jj,w,v,k11,k22,k33,k44
    real(kind=8), external ::funcdj
    real(kind=8), external ::funcfs

    open(unit=10,file='flya.dat')
    open(unit=11,file='x-y.dat')

    fnorm=0.0D0
    flya=0.0D0

    do i=1,m,1
        do j=1,m,1
            if(i==j) then
                w(i,j)=1.0D0
            else
                w(i,j)=0.0D0
            end if
        end do
    end do

    v=0.0D0

    do ii=2200,2200,1
        call initial(y)
        !cc=-5.70
        do step=1,n,1

            !-----------Evolution of Single Lorenz system-------------------------
            call fnf1(y,k1)
            call fnf1(y+K1/2,K2)
            call fnf1(y+K2/2,K3)
            call fnf1(y+K3,K4)
            y=y+(k1+2.0D0*k2+2.0D0*k3+k4)/6.0d0

            !---------Jacobian Matrix of Single Lorenz System-------------------
            jj(1,1)=-aa
            jj(1,2)=aa
            jj(1,3)=0.0D0
            jj(2,1)=-y(3)
            jj(2,2)=-cc
            jj(2,3)=-y(1)
            jj(3,1)=y(2)
            jj(3,2)=y(1)
            jj(3,3)=0.0D0
            !-------------------------------------------------------------------

            !-----------Evolution of Jacobian Matrix----------------------------
            call fnf2(w,k11,jj)
            call fnf2(w+K11/2,K22,jj)
            call fnf2(w+K22/2,K33,jj)
            call fnf2(w+K33,K44,jj)
            w=w+(k11+2.0D0*k22+2.0D0*k33+k44)/6.0d0
            !---------------------------------------------------------------------

            !-----------Gram-Schmit Normalization-----------------------------------
            do j=2,m,1
                do i=1,j-1,1
                    w(:,j)=w(:,j)-funcdj(w(:,j),w(:,i))/funcfs(w(:,i))*w(:,i)
                end do
            end do
            !-------------------------------------------------------------------------

            !-----------Calculate the module of every vector direction----------------
            do j=1,m,1
                fnorm(j)=funcfs(w(:,j))
                fnorm(j)=sqrt(fnorm(j))
                w(:,j)=w(:,j)/fnorm(j) !---------Renormalization every vector
                fnorm(j)=log(fnorm(j))
            end do
            !--------------------------------------------------------------------------

            !----------Calculate Lyaponov Exponent of every vector direction------------
            if(step>st) then
                do i=1,m,1
                    flya(i)=flya(i)+fnorm(i)
                end do
            end if
            !----------------------------------------------------------------------------

            if(step>n-st) then
                write(11,'(f15.9,2x,f15.9,2x,f15.9)') y(1),y(2),y(3)
            end if
        end do


        flya=flya/(n-st)/dt

        write(10,'(f15.9,2x,f15.9,2x,f15.9,2x,f15.9)') cc,flya(1),flya(2),flya(3)
    end do
    stop
end program

!--------------------calculate dot metrix of two vector: a.b-------------------
real(kind=8) function funcdj(u,v)
    use head
    implicit none
    real(kind=8), dimension(m) :: u, v
    integer(kind=4) :: i

    funcdj=0.0d0
    do i=1,m,1
        funcdj=funcdj+u(i)*v(i)
    end do

    return
end function
!--------------------------------------------------------------------------------

!------------------calculate module of one vector direction---------------------
real(kind=8) function funcfs(u)
    use head
    implicit none
    real(kind=8), dimension(m) :: u
    integer(kind=4) :: i

    funcfs=0.0d0
    do i=1,m,1
        funcfs=funcfs+u(i)*u(i)
    end do

    return
end function
!-------------------------------------------------------------------------------

image-20221114185637365

image-20221114185742338

image-20221114185811739

文章作者: rep-rebirth
文章链接: https://zhangwenhao.icu/2022/11/11/lyapunov/
版权声明: 本博客所有文章除特別声明外，均采用 CC BY 4.0 许可协议。转载请注明来源 rep-rebirth !
系统动力学 Lyapunov 混沌
赏
评论
评论
  目录