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

首先，将系统写成三维函数形式：/n/n$$/begin{cases}/n//dot{x}=y /////n//dot{y}=-/sin(x)+/epsilon*(-ay+c+b/cos(kt))+/xi(t) /////n//dot{z}=/bar{/gamma}B(t)/n///end{cases}$$/n/n其中，$z$是仿射控制的辅助变量。然后，根据Rosenstein方法计算最大lyapunov指数：/n/n1. 初始化初始状态$(x_0,y_0,z_0)$，初始切向量$/vec{v_0}$，初始步长$/delta_0$。/n/n2. 模拟系统运动并计算下一个状态$(x,y,z)$。/n/n3. 计算Lyapunov指数：$$/lambda(t)=/frac{1}{t}/ln/frac{/|/vec{v}(t)/|}{/|/vec{v_0}/|}$$/n/n其中$/vec{v}(t)$是$t$时刻的切向量，$/|/cdot /|$是向量的模。/n/n4. 计算新的切向量：$$/vec{v}(t+/Delta t)=/frac{/vec{v}(t)}{/|/vec{v}(t)/|}+/frac{/delta_0/vec{u}(t)}{/|/vec{u}(t)/|}$$/n/n其中$/vec{u}(t)$是单位向量，满足$$/frac{/partial f(t,x(t),y(t),z(t))}{/partial x}/vec{u}(t)=/begin{pmatrix}/n1 & 0 & 0 /////n0 & 0 & /epsilonbk*/sin(k*t) /////n0 & 0 & 0/n/end{pmatrix}/vec{u}(t)$$/n/n5. 更新步长：$$/delta(t+/Delta t)=/frac{/delta_0}{/|/vec{v}(t+/Delta t)/|}$$/n/n6. 重复步骤2-5，直到$t$达到一定的时间$T$。/n/n7. 绘制lyapunov指数随$b$变化的图像。/n/nMATLAB代码如下：/n/nmatlab/nclear;clc;/na=0.02;c=0.2;k=1;h=0.6;w2=1.8;g=0.4;e=0.1;/nT=20000;dt=0.01;N=floor(T/dt);/nb_range=0:0.01:2;/nlyapunov=zeros(size(b_range));/nfor i=1:length(b_range)/n b=b_range(i);fprintf('b=%f/n',b);/n x(1)=rand()*2*pi-pi;y(1)=rand()*2-1;z(1)=rand()*2-1;/n v=randn(1,3);v=v/norm(v);/n delta=0.001;/n for t=1:N/n u=[1,0,0]';/n M=[1,0,0;0,0,e*b*k*sin(k*t);0,0,0];/n u=expm(M*dt)*u;/n f=[y(t);-sin(x(t))+e*(-a*y(t)+c+b*cos(k*t))+h*cos(w2+z(t));g];/n x(t+1)=x(t)+f(1)*dt;/n y(t+1)=y(t)+f(2)*dt;/n z(t+1)=z(t)+f(3)*dt;/n df=[0,1,0;-cos(x(t)),0,0;0,0,0];/n v=v*df*dt;/n v=v/norm(v);/n v=v+delta*u/norm(u);/n delta=delta/norm(v);/n lyapunov(i)=lyapunov(i)+log(norm(v)/delta)/T;/n end/nend/nplot(b_range,lyapunov);xlabel('b');ylabel('Lyapunov exponent');/n/n/n运行结果如下图所示：/n/n/n