基于RHR算法的求Lyapunov程序
- Lya1=[];Lya2=[];Lya3=[];
- V=eye(3);
- S=V;b1=0;
- a=0.4;c=0.2;gama=3.5;
- b=4.0;
- h=0.01;
- x(1)=0.1;y(1)=0;z(1)=0;n=0;
- while z<=200
- n=n+1;
- k1=h*y(n);
- m1=h*(-sin(x(n))-a*y(n)+b*cos(gama*z(n)).*sin(x(n))+c);
- k2=h*(y(n)+m1/2);
- m2=h*(-sin(x(n)+k1/2)-a*(y(n)+m1/2)+b*cos(gama*(z(n)+h/2)).*sin(x(n)+k1/2)+c);
- k3=h*(y(n)+m2/2);
- m3=h*(-sin(x(n)+k2/2)-a*(y(n)+m2/2)+b*cos(gama*(z(n)+h/2)).*sin(x(n)+k2/2)+c);
- k4=h*(y(n)+m3);
- m4=h*(-sin(x(n)+k3)-a*(y(n)+m3)+b*cos(gama*(z(n)+h)).*sin(x(n)+k3)+c);
- x(n+1)=x(n)+(k1+2*k2+2*k3+k4)/6;
- y(n+1)=y(n)+(m1+2*m2+2*m3+m4)/6;
- z(n+1)=n*h;
- J = [0 1 0;
- b*cos(gama*z(n+1))*cos(x(n+1))-cos(x(n+1)) -a -b*gama*sin(gama*z(n+1))*sin(x(n+1));
- 0 0 0];
- J=eye(3)+h*J;
- B=J*V*S;
- [V,S,U]=svd(B);
- a_max=max(diag(S));
- S=(1/a_max)*S;
- b1=b1+log(a_max);
- Lyapunov1=(log(diag(S))+b1)/(n*h);
- Lya1=[Lya1,Lyapunov1(1,:)];
- Lya2=[Lya2,Lyapunov1(2,:)];
- Lya3=[Lya3,Lyapunov1(3,:)];
- end
- Lyapunov1
- n=1:20001;
- plot(n,Lya1,'k',n,Lya2,'k',n,Lya3,'k')
- %grid on
- axis([0,30001,-0.8,0.5])
- title('Lyapunov exponents of Warship')
- xlabel('n'),ylabel('LCE')
复制代码
|