随机系统最大Lyapunov指数随参数变化的MATLAB仿真

本文使用MATLAB仿真一个随机系统，该系统可以描述为：

$\xi(t)=h*cos(\omega_2+\bar{\gamma}*B(t))$，其中 $B(t)$ 是单位维纳过程。

$\ddot{x}+sin(x)-\epsilon*(-a*\dot{x}+c+bcos(kt))=\xi(t)$。

我们取 $a=0.02, c=0.2, k=1, h=0.6, w_2=1.8, \bar{\gamma}=0.4, \epsilon=0.1$，并使用Rosenstein方法计算最大Lyapunov指数随参数b的变化曲线。

为了便于仿真，我们将系统写成三维形式：

$\dot{x} = y$

$\dot{y} = -x + \gammahcos(w_2 + \gamma*y)$

$\dot{z} = x$

以下是MATLAB代码：

function [t,lyapunov]=lyapunov_exponent()
% Parameters
a=0.02; c=0.2; k=1; h=0.6; w2=1.8; gamma=0.4; epsilon=0.1;
b_range=0:0.01:2; % range of b values
dt=0.01; % time step
t0=0; tf=1000; % initial and final time
tol=1e-6; % tolerance for ODE solver
% Initial conditions
x0=0; xp0=0; y0=0;
lyapunov=zeros(length(b_range),1);
for i=1:length(b_range)
b=b_range(i);
% Define the ODE system
f=@(t,xp,x,y) [y+epsilon*a*xp-epsilon*c*x-epsilon*b*cos(k*t);-x+gamma*h*cos(w2+gamma*y);x];
% Solve the ODE system using Runge-Kutta 4
x=zeros(3,1);
xp=zeros(3,1);
y=zeros(3,1);
x(1)=x0;
xp(1)=xp0;
y(1)=y0;
for j=1:2
t=t0+(j-1)*dt;
k1=f(t,xp(j),x(j),y(j));
k2=f(t+dt/2,xp(j)+dt*k1(1)/2,x(j)+dt*k1(2)/2,y(j)+dt*k1(3)/2);
k3=f(t+dt/2,xp(j)+dt*k2(1)/2,x(j)+dt*k2(2)/2,y(j)+dt*k2(3)/2);
k4=f(t+dt,xp(j)+dt*k3(1),x(j)+dt*k3(2),y(j)+dt*k3(3));
xp(j+1)=xp(j)+(dt/6)*(k1(1)+2*k2(1)+2*k3(1)+k4(1));
x(j+1)=x(j)+(dt/6)*(k1(2)+2*k2(2)+2*k3(2)+k4(2));
y(j+1)=y(j)+(dt/6)*(k1(3)+2*k2(3)+2*k3(3)+k4(3));
end
% Calculate the tangent vector
tangent_vector=[x(2)-x(1); xp(2)-xp(1); y(2)-y(1)];
tangent_vector=tangent_vector/norm(tangent_vector);
% Perturb the initial conditions for the next iteration
x0=x(2)+tol*tangent_vector(2);
xp0=xp(2)+tol*tangent_vector(1);
y0=y(2)+tol*tangent_vector(3);
% Calculate the Lyapunov exponent
lyapunov(i)=1/(2*dt*i)*log(norm(tangent_vector)/tol);
end
% Plot the results
plot(b_range,lyapunov);
xlabel('b');
ylabel('Lyapunov exponent');
end

运行该函数，即可得到随参数b变化的最大Lyapunov指数图。

本文使用Rosenstein方法计算Lyapunov指数，该方法通过扰动初值并观察扰动随时间的增长来计算Lyapunov指数。代码中，我们首先使用四阶龙格库塔方法求解系统微分方程，然后通过扰动初始条件并计算两个解之间的距离来计算Lyapunov指数。

该仿真结果可以帮助我们理解参数b对系统动力学的影响，并为进一步研究该系统的混沌特性提供参考。