lyapunov_rhr.txt

现我做了一个基于RHR算法的李氏指数计算方法

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(;
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')