tic
clear
clc
=lyapunov(4,@lorenz_ext,@ode15s,0,0.5,0.5,',10);
%plot(T,Res);
figure;
plot(T,Res(:,1));
title('Dynamics of Lyapunov exponents x');
xlabel('Time'); ylabel('Lyapunov exponents x');
figure;
plot(T,Res(:,3));
title('Dynamics of Lyapunov exponents y');
xlabel('Time'); ylabel('Lyapunov exponents y');
toc
function =lyapunov(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp);
%
% Lyapunov exponent calcullation for ODE-system.
%
% The alogrithm employed in this m-file for determining Lyapunov
% exponents was proposed in
%
% A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano,
% "Determining Lyapunov Exponents from a Time Series," Physica D,
% Vol. 16, pp. 285-317, 1985.参考的论文
%
% For integrating ODE system can be used any MATLAB ODE-suite methods.
% This function is a part of MATDS program - toolbox for dynamical system investigation
% See: http://www.math.rsu.ru/mexmat/kvm/matds/相关的网址
%
% Input parameters:输入参数
% n - number of equation微分方程个数
% rhs_ext_fcn - handle of function with right hand side of extended ODE-system.
% This function must include RHS of ODE-system coupled with
% variational equation (n items of linearized systems, see Example). 定义微分方程组的函数句柄名
% fcn_integrator - handle of ODE integrator function, for example: @ode45 所调用求微分方程的MATLAB函数名
% tstart - start values of independent value (time t)时间的初始值
% stept - step on t-variable for Gram-Schmidt renormalization
% procedure.进行Gram_Schmidt正交化的迭代次数
% tend - finish value of time时间的终止数值
% ystart - start point of trajectory of ODE system.求解微分方程的初值
% ioutp - step of print to MATLAB main window. ioutp==0 - no print,
% if ioutp>0 then each ioutp-th point will be
% print.每隔ioutp个点输出相应的时间和李雅普诺夫指数数值
%
% Output parameters:输出参数
% Texp - time values时间的数值
% Lexp - Lyapunov exponents to each time value.李雅普诺夫指数的数值
%
% Users have to write their own ODE functions for their specified
% systems and use handle of this function as rhs_ext_fcn - parameter.
%用户如果要计算其他方程，需要写出自己的常微分方程组来替换lorenz_ext文件内容得到相应的函数表达式。
% Example. Lorenz system:
% dx/dt = sigma*(y - x) = f1
% dy/dt = r*x - y - x*z = f2
% dz/dt = x*y - b*z = f3
%
% The Jacobian of system:
% | -sigmasigma0 |
% J = | r-z -1 -x |
% | y x -b |
%
% Then, the variational equation has a form:
%
% F = J*Y
% where Y is a square matrix with the same dimension as J.
% Corresponding m-file:
% function f=lorenz_ext(t,X)
% SIGMA = 10; R = 28; BETA = 8/3;
% x=X(1); y=X(2); z=X(3);
%
% Y= [X(4), X(7), X(10);
% X(5), X(8), X(11);
% X(6), X(9), X(12)];
% f=zeros(9,1);
% f(1)=SIGMA*(y-x); f(2)=-x*z+R*x-y; f(3)=x*y-BETA*z;
%
% Jac=[-SIGMA,SIGMA,0; R-z,-1,-x; y, x,-BETA];
%
% f(4:12)=Jac*Y;
%
% Run Lyapunov exponent calculation:
%
% =lyapunov(3,@lorenz_ext,@ode45,0,0.5,200,,10);
%
% See files: lorenz_ext, run_lyap.
%
% --------------------------------------------------------------------
% Copyright (C) 2004, Govorukhin V.N.
% This file is intended for use with MATLAB and was produced for MATDS-program
% http://www.math.rsu.ru/mexmat/kvm/matds/
% lyapunov.m is free software. lyapunov.m is distributed in the hope that it
% will be useful, but WITHOUT ANY WARRANTY.
%
%
% n=number of nonlinear odes非线性方程个数
% n2=n*(n+1)=total number of odes
%
n1=n; n2=n1*(n1+1);
%Number of steps步数
nit = round((tend-tstart)/stept);
% Memory allocation 内存分配
y=zeros(n2,1); cum=zeros(n1,1); y0=y;
gsc=cum; znorm=cum;
% Initial values初值
y(1:n)=ystart(:);
for i=1:n1 y((n1+1)*i)=1.0; end;
t=tstart;
% Main loop
for ITERLYAP=1:nit
% Solutuion of extended ODE system
= feval(fcn_integrator,rhs_ext_fcn,,y);
t=t+stept;
y=Y(size(Y,1),:);
for i=1:n1
for j=1:n1 y0(n1*i+j)=y(n1*j+i); end;
end;
%
% construct new orthonormal basis by gram-schmidt
%
znorm(1)=0.0;
for j=1:n1 znorm(1)=znorm(1)+y0(n1*j+1)^2; end;
znorm(1)=sqrt(znorm(1));
for j=1:n1 y0(n1*j+1)=y0(n1*j+1)/znorm(1); end;
for j=2:n1
for k=1:(j-1)
gsc(k)=0.0;
for l=1:n1 gsc(k)=gsc(k)+y0(n1*l+j)*y0(n1*l+k); end;
end;
for k=1:n1
for l=1:(j-1)
y0(n1*k+j)=y0(n1*k+j)-gsc(l)*y0(n1*k+l);
end;
end;
znorm(j)=0.0;
for k=1:n1 znorm(j)=znorm(j)+y0(n1*k+j)^2; end;
znorm(j)=sqrt(znorm(j));
for k=1:n1 y0(n1*k+j)=y0(n1*k+j)/znorm(j); end;
end;
%
% update running vector magnitudes
%
for k=1:n1 cum(k)=cum(k)+log(znorm(k)); end;
%
% normalize exponent
%
for k=1:n1
lp(k)=cum(k)/(t-tstart);
end;
% Output modification
if ITERLYAP==1
Lexp=lp;
Texp=t;
else
Lexp=;
Texp=;
end;
if (mod(ITERLYAP,ioutp)==0)
fprintf('t=%6.4f',t);
for k=1:n1 fprintf(' %10.6f',lp(k)); end;
fprintf('\n');
end;
for i=1:n1
for j=1:n1
y(n1*j+i)=y0(n1*i+j);
end;
end;
end;
function f=lorenz_ext(t,X)
%
%Lorenz equation
%
% dx/dt = SIGMA*(y - x)
% dy/dt = R*x - y -x*z
% dz/dt= x*y - BETA*z
%
% In demo run SIGMA = 10, R = 28, BETA = 8/3
% Initial conditions: x(0) = 0, y(0) = 1, z(0) = 0;
% Reference values for t=10 000 :
% L_1 = 0.9022, L_2 = 0.0003, LE3 = -14.5691
%
% See:
% K. Ramasubramanian, M.S. Sriram, "A comparative study of computation
% of Lyapunov spectra with different algorithms", Physica D 139 (2000) 72-86.
%
% --------------------------------------------------------------------
% Copyright (C) 2004, Govorukhin V.N.
% Values of parameters
%SIGMA = 10;
%R = 28;
%BETA = 8/3;
%x=X(1)
%y=X(2)
%z=X(3)
%赋初值
r=0.0003;%m
l1=0.130;%m
l2=0.308;%m
l4=0.520;%m
ls2=0.154;%m
ls3=0.0797;%m
m2=8.516;%kg
m3=1.161;%kg
%k=7.15e4;%N/mm
Js2=0.0616;%转动惯量kg.mm^2
Js3=0.005772;%转动惯量kg.mm^2
J4=0.007105;
if sqrt(X(1)^2+X(3)^2)>r
k=7.15e10;%N/m
else
k=0;
end
n=10;
theta1=n*pi*t;
theta11=n*pi;
theta12=0;
aa1=-l1*sin(theta1);
aa2=l4-l1*cos(theta1);
if aa1<0
theta20=-asin(l2/sqrt(aa1^2+aa2^2))-atan(aa2/aa1);
else
theta20=pi-asin(l2/sqrt(aa1^2+aa2^2))-atan(aa2/aa1);
end
l0=(l1*sin(theta1)+l2*sin(theta20))/cos(theta20);
a1=-(l2*cos(theta20)+l0*sin(theta20))/l0;
a2=-(l2*sin(theta20)-l0*cos(theta20))/l0;
b1=-cos(theta20)./l0;
b2=-sin(theta20)./l0;
theta201=(l1*theta11*sin(theta1-theta20))/(l1*sin(theta1)*cos(theta20)+(l4-l1*cos(theta1))*sin(theta20));
theta202=(l1*theta12*sin(theta1-theta20)-theta201^2*(l4-l1*cos(theta1))*cos(theta20)+theta201^2*l1*sin(theta1)*sin(theta20)+l1*theta11^2*cos(theta1-theta20))/(l1*sin(theta1)*cos(theta20)+(l4-l1*cos(theta1))*sin(theta20));
l01=(l1*theta11*cos(theta1)+l2*theta201*cos(theta20)+l0*theta201*sin(theta20))/cos(theta20);
l02=(l1*theta12*cos(theta1)-l1*theta11^2*sin(theta1)+l2*theta202*cos(theta20)-l2*theta201^2*sin(theta20)+2*l01*theta201*sin(theta20)+l0*theta202*sin(theta20)+l0*theta201^2*cos(theta20))/cos(theta20);
b11=(l0*theta201*sin(theta20)+l01*cos(theta20))/l0^2;
b12=(-2*b11*l0*l01+l02*cos(theta20)+l0*theta201^2*cos(theta20))/l0^2;
b21=(l01*sin(theta20)-l0*theta201*cos(theta20))/l0^2;
b22=(l02*sin(theta20)-l0*theta202*cos(theta20)+l0*theta201^2*sin(theta20)-2*l0*l01*b21)/l0^2;
a11=l2*b11-theta201*cos(theta20);
a12=l2*b12-theta202*cos(theta20)+theta201^2*sin(theta20);
a21=l2*b21-theta201*sin(theta20);
a22=l2*b22-theta202*sin(theta20)--theta201^2*cos(theta20);
%dtheta2=b1*X(1)+b2*X(3);
dl=a1*X(1)+a2*X(3);
dl1=a11*X(1)+a1*X(2)+a21*X(3)+a2*X(4);
dtheta21=b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4);
%theta2=theta20+b1*X(1)+b2*X(3);
A=Js2+Js3+J4+m2*ls2^2+m3*ls3^2;
B=(m2+m3)*(l0+dl)^2-2*m3*ls3*(l0+dl);
A1=(A+B)*b1^2-2*m2*ls2*a1*b1+(m2+m3)*a1^2;
B1=(A+B)*b1*b2-m2*ls2*(a1*b2+b1*a2)+(m2+m3)*a1*a2;
A2=(A+B)*b1*b2-m2*ls2*(a1*b2+b1*a2)+(m2+m3)*a1*a2;
B2=(A+B)*b2^2-2*m2*ls2*a2*b2+(m2+m3)*a2^2;
H=((A+B)*(m2+m3)-m2^2*ls2^2)/(l0^2);
C=2*((m2+m3)*(l0+dl)-m3*ls3)*(l01+dl1)*(theta201+dtheta21);
D=theta202+2*b11*X(2)+2*b21*X(4)+b12*X(1)+b22*X(3);
E=l02+2*a11*X(2)+2*a21*X(4)+a12*X(1)+a22*X(3);
F=((m2+m3)*(l0+dl)-m3*ls3)*(theta201+dtheta21)^2;
Y= [X(5), X(9),X(13), X(17);
X(6), X(10), X(14), X(18);
X(7), X(11), X(15), X(19);
X(8), X(12), X(16), X(20)];
f=zeros(20,1);
%Lorenz equation
%f(1)=SIGMA*(y-x);
%f(2)=-x*z+R*x-y/(x^2+5);
%f(3)=x*y-BETA*z/(y^2+5);
f(1)=X(2);
f(3)=X(4);
f(2)=-D*(l2*sin(theta20)-l0*cos(theta20))+E*sin(theta20)-C/H*(m2*ls2*(-sin(theta20)/(l0^2))+(m2+m3)*(-a2/l0))-F/H*((A+B)*(-sin(theta20)/(l0^2))+m2*ls2*(-a2/l0))-k*(1-r/sqrt(X(1)^2+X(3)^2))/H*(X(1)*B2-X(3)*B1);
f(4)=D*(l2*cos(theta20)+l0*sin(theta20))-E*cos(theta20)+C/H*(m2*ls2*(-cos(theta20)/(l0^2))+(m2+m3)*(-a1/l0))+F/H*((A+B)*(-cos(theta20)/(l0^2))+m2*ls2*(-a1/l0))+k*(1-r/sqrt(X(1)^2+X(3)^2))/H*(X(1)*A2-X(3)*A1);
%Linearized system
%Jac=[-SIGMA, SIGMA, 0;
% R-z, -1/(x^2+5), -x;
% y, x, -BETA/(y^2+5)];
%f=f(2),g=f(4),X1=X(1),X2=X(2),X3=X(3),X4=X(4)
%J21=diff(f,'X1')
%J22=diff(f,'X2')
%J23=diff(f,'X3')
%J24=diff(f,'X4')
%J41=diff(g,'X1')
%J42=diff(g,'X2')
%J43=diff(g,'X3')
%J44=diff(g,'X4')
I=l2*sin(theta20)-l0*cos(theta20);
I1=l2*cos(theta20)+l0*sin(theta20);
J=m2*ls2*(-sin(theta20)/(l0^2))+(m2+m3)*(-a2/l0);
J1=m2*ls2*(-cos(theta20)/(l0^2))+(m2+m3)*(-a1/l0);
K=(A+B)*(-sin(theta20)/(l0^2))+m2*ls2*(-a2/l0);
K1=(A+B)*(-cos(theta20)/(l0^2))+m2*ls2*(-a1/l0);
P=k*(1-r/sqrt(X(1)^2+X(3)^2))*(X(1)*B2-X(3)*B1);
P1=k*(1-r/sqrt(X(1)^2+X(3)^2))*(X(1)*A2-X(3)*A1);
Bx=2*((m2+m3)*(l0+dl)-m3*ls3)*a1;
By=2*((m2+m3)*(l0+dl)-m3*ls3)*a2;
Hx=Bx*(m2+m3)/l0^2;
Hy=By*(m2+m3)/l0^2;
Cx=2*(m2+m3)*a1*(l01+dl1)*(theta201+dtheta21)+2*((m2+m3)*(l0+dl)-m3*ls3)*(a11*(theta201+dtheta21)+(l01+dl1)*b11);
Cy=2*(m2+m3)*a2*(l01+dl1)*(theta201+dtheta21)+2*((m2+m3)*(l0+dl)-m3*ls3)*(a21*(theta201+dtheta21)+(l01+dl1)*b21);
Cx1=2*((m2+m3)*(l0+dl)-m3*ls3)*((l01+dl1)*b1+(theta201+dtheta21)*a1);
Cy1=2*((m2+m3)*(l0+dl)-m3*ls3)*((l01+dl1)*b2+(theta201+dtheta21)*a2);
Fx=(m2+m3)*a1*(theta201+dtheta21)^2+2*((m2+m3)*(l0+dl)-m3*ls3)*(theta201+dtheta21)*b11;
Fy=(m2+m3)*a2*(theta201+dtheta21)^2+2*((m2+m3)*(l0+dl)-m3*ls3)*(theta201+dtheta21)*b21;
Fx1=2*((m2+m3)*(l0+dl)-m3*ls3)*(theta201+dtheta21)*b1;
Fy1=2*((m2+m3)*(l0+dl)-m3*ls3)*(theta201+dtheta21)*b2;
Kx=Bx*(-sin(theta20)/(l0^2));
Ky=By*(-sin(theta20)/(l0^2));
K1x=Bx*(-cos(theta20)/(l0^2));
K1y=By*(-cos(theta20)/(l0^2));
B2x=Bx*b2^2;
B1x=Bx*b2*b1;
B2y=By*b2^2;
B1y=By*b2*b1;
A2x=Bx*b2*b1;
A1x=Bx*b1^2;
A2y=By*b2*b1;
A1y=By*b1^2;
Px=k*(1-r/sqrt(X(1)^2+X(3)^2))*(B2+B2x*X(1)-X(3)*B1x)+k*r*X(1)*(X(1)*B2-X(3)*B1)/((X(1)^2+X(3)^2)^(3/2));
Py=k*(1-r/sqrt(X(1)^2+X(3)^2))*(B2y*X(1)-X(3)*B1y-B1)+k*r*X(3)*(X(1)*B2-X(3)*B1)/((X(1)^2+X(3)^2)^(3/2));
P1x=k*(1-r/sqrt(X(1)^2+X(3)^2))*(A2+A2x*X(1)-X(3)*A1x)+k*r*X(1)*(X(1)*A2-X(3)*A1)/((X(1)^2+X(3)^2)^(3/2));
P1y=k*(1-r/sqrt(X(1)^2+X(3)^2))*(A2y*X(1)-X(3)*A1y-A1)+k*r*X(3)*(X(1)*A2-X(3)*A1)/((X(1)^2+X(3)^2)^(3/2));
J21=-b12*I+a12*sin(theta20)-(Cx*J+Fx*K+Kx*F+Px)/H+(C*J+F*K+P)*Hx/H^2;
J22=-2*b11*I+2*a11*sin(theta20)-(Cx1*J+Fx1*K)/H;
J23=-b22*I+a22*sin(theta20)-(Cy*J+Fy*K+Ky*F+Py)/H+(C*J+F*K+P)*Hy/H^2;
J24=-2*b21*I+2*a21*sin(theta20)-(Cy1*J+Fy1*K)/H;
J41=b12*I1-a12*cos(theta20)+(Cx*J1+Fx*K1+K1x*F+P1x)/H-(C*J1+F*K1+P1)*Hx/H^2;
J42=2*b11*I1-2*a11*cos(theta20)+(Cx1*J1+Fx1*K1)/H;
J43=b22*I1-a22*cos(theta20)+(Cy*J1+Fy*K1+K1y*F+P1y)/H-(C*J1+F*K1+P1)*Hy/H^2;
J44=2*b21*I1-2*a21*cos(theta20)+(Cy1*J1+Fy1*K1)/H;
%J21 =-b12*(l2*sin(theta20)-l0*cos(theta20))+a12*sin(theta20)-(2*(m2+m3)*a1*(l01+a11*X(1)+a1*X(2)+a21*X(3)+a2*X(4))+2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a11-2*m3*ls3*a11*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))-2*m3*ls3*(l01+a11*X(1)+a1*X(2)+a21*X(3)+a2*X(4))*b11)*(-m2*ls2*sin(theta20)/l0^2-(m2+m3)*a2/l0)-(m2+m3)*a1*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))^2/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*l0^2*(-(A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*sin(theta20)/l0^2-m2*ls2*a2/l0)-2*((m2+m3)*(l0+a1*X(1)+a2*X(3))-m3*ls3)*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*l0^2*(-(A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*sin(theta20)/l0^2-m2*ls2*a2/l0)*b11+((m2+m3)*(l0+a1*X(1)+a2*X(3))-m3*ls3)*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))^2/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)^2*l0^2*(-(A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*sin(theta20)/l0^2-m2*ls2*a2/l0)*(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a1-2*m3*ls3*a1)*(m2+m3)+((m2+m3)*(l0+a1*X(1)+a2*X(3))-m3*ls3)*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))^2/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a1-2*m3*ls3*a1)*sin(theta20)-k*r/(X(1)^2+X(3)^2)^(3/2)*X(1)/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*l0^2*(X(1)*((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*b2^2-2*m2*ls2*a2*b2+(m2+m3)*a2^2)-X(3)*((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*b1*b2-m2*ls2*(a1*b2+b1*a2)+(m2+m3)*a1*a2))+k*(1-r/(X(1)^2+X(3)^2)^(1/2))/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)^2*l0^2*(X(1)*((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*b2^2-2*m2*ls2*a2*b2+(m2+m3)*a2^2)-X(3)*((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*b1*b2-m2*ls2*(a1*b2+b1*a2)+(m2+m3)*a1*a2))*(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a1-2*m3*ls3*a1)*(m2+m3)-k*(1-r/(X(1)^2+X(3)^2)^(1/2))/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*l0^2*((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*b2^2-2*m2*ls2*a2*b2+(m2+m3)*a2^2+X(1)*(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a1-2*m3*ls3*a1)*b2^2-X(3)*(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a1-2*m3*ls3*a1)*b1*b2);
%J22 =-2*b11*(l2*sin(theta20)-l0*cos(theta20))+2*a11*sin(theta20)-(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a1-2*m3*ls3*a1*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))-2*m3*ls3*(l01+a11*X(1)+a1*X(2)+a21*X(3)+a2*X(4))*b1)*(-m2*ls2*sin(theta20)/l0^2-(m2+m3)*a2/l0)-2*((m2+m3)*(l0+a1*X(1)+a2*X(3))-m3*ls3)*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*l0^2*(-(A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*sin(theta20)/l0^2-m2*ls2*a2/l0)*b1;
%J23 =-b22*(l2*sin(theta20)-l0*cos(theta20))+a22*sin(theta20)-(2*(m2+m3)*a2*(l01+a11*X(1)+a1*X(2)+a21*X(3)+a2*X(4))+2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a21-2*m3*ls3*a21*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))-2*m3*ls3*(l01+a11*X(1)+a1*X(2)+a21*X(3)+a2*X(4))*b21)*(-m2*ls2*sin(theta20)/l0^2-(m2+m3)*a2/l0)-(m2+m3)*a2*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))^2/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*l0^2*(-(A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*sin(theta20)/l0^2-m2*ls2*a2/l0)-2*((m2+m3)*(l0+a1*X(1)+a2*X(3))-m3*ls3)*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*l0^2*(-(A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*sin(theta20)/l0^2-m2*ls2*a2/l0)*b21+((m2+m3)*(l0+a1*X(1)+a2*X(3))-m3*ls3)*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))^2/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)^2*l0^2*(-(A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*sin(theta20)/l0^2-m2*ls2*a2/l0)*(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a2-2*m3*ls3*a2)*(m2+m3)+((m2+m3)*(l0+a1*X(1)+a2*X(3))-m3*ls3)*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))^2/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a2-2*m3*ls3*a2)*sin(theta20)-k*r/(X(1)^2+X(3)^2)^(3/2)*X(3)/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*l0^2*(X(1)*((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*b2^2-2*m2*ls2*a2*b2+(m2+m3)*a2^2)-X(3)*((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*b1*b2-m2*ls2*(a1*b2+b1*a2)+(m2+m3)*a1*a2))+k*(1-r/(X(1)^2+X(3)^2)^(1/2))/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)^2*l0^2*(X(1)*((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*b2^2-2*m2*ls2*a2*b2+(m2+m3)*a2^2)-X(3)*((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*b1*b2-m2*ls2*(a1*b2+b1*a2)+(m2+m3)*a1*a2))*(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a2-2*m3*ls3*a2)*(m2+m3)-k*(1-r/(X(1)^2+X(3)^2)^(1/2))/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*l0^2*(X(1)*(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a2-2*m3*ls3*a2)*b2^2-(A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*b1*b2+m2*ls2*(a1*b2+b1*a2)-(m2+m3)*a1*a2-X(3)*(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a2-2*m3*ls3*a2)*b1*b2);
%J24 =-2*b21*(l2*sin(theta20)-l0*cos(theta20))+2*a21*sin(theta20)-(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a2-2*m3*ls3*a2*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))-2*m3*ls3*(l01+a11*X(1)+a1*X(2)+a21*X(3)+a2*X(4))*b2)*(-m2*ls2*sin(theta20)/l0^2-(m2+m3)*a2/l0)-2*((m2+m3)*(l0+a1*X(1)+a2*X(3))-m3*ls3)*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*l0^2*(-(A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*sin(theta20)/l0^2-m2*ls2*a2/l0)*b2;
%J41 =b12*(l2*cos(theta20)+l0*sin(theta20))-a12*cos(theta20)+(2*(m2+m3)*a1*(l01+a11*X(1)+a1*X(2)+a21*X(3)+a2*X(4))+2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a11-2*m3*ls3*a11*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))-2*m3*ls3*(l01+a11*X(1)+a1*X(2)+a21*X(3)+a2*X(4))*b11)*(-m2*ls2*cos(theta20)/l0^2-(m2+m3)*a1/l0)+(m2+m3)*a1*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))^2/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*l0^2*(-(A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*cos(theta20)/l0^2-m2*ls2*a1/l0)+2*((m2+m3)*(l0+a1*X(1)+a2*X(3))-m3*ls3)*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*l0^2*(-(A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*cos(theta20)/l0^2-m2*ls2*a1/l0)*b11-((m2+m3)*(l0+a1*X(1)+a2*X(3))-m3*ls3)*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))^2/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)^2*l0^2*(-(A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*cos(theta20)/l0^2-m2*ls2*a1/l0)*(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a1-2*m3*ls3*a1)*(m2+m3)-((m2+m3)*(l0+a1*X(1)+a2*X(3))-m3*ls3)*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))^2/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a1-2*m3*ls3*a1)*cos(theta20)+k*r/(X(1)^2+X(3)^2)^(3/2)*X(1)/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*l0^2*(X(1)*((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*b1*b2-m2*ls2*(a1*b2+b1*a2)+(m2+m3)*a1*a2)-X(3)*((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*b1^2-2*m2*ls2*a1*b1+(m2+m3)*a1^2))-k*(1-r/(X(1)^2+X(3)^2)^(1/2))/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)^2*l0^2*(X(1)*((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*b1*b2-m2*ls2*(a1*b2+b1*a2)+(m2+m3)*a1*a2)-X(3)*((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*b1^2-2*m2*ls2*a1*b1+(m2+m3)*a1^2))*(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a1-2*m3*ls3*a1)*(m2+m3)+k*(1-r/(X(1)^2+X(3)^2)^(1/2))/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*l0^2*((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*b1*b2-m2*ls2*(a1*b2+b1*a2)+(m2+m3)*a1*a2+X(1)*(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a1-2*m3*ls3*a1)*b1*b2-X(3)*(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a1-2*m3*ls3*a1)*b1^2);
%J42 =2*b11*(l2*cos(theta20)+l0*sin(theta20))-2*a11*cos(theta20)+(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a1-2*m3*ls3*a1*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))-2*m3*ls3*(l01+a11*X(1)+a1*X(2)+a21*X(3)+a2*X(4))*b1)*(-m2*ls2*cos(theta20)/l0^2-(m2+m3)*a1/l0)+2*((m2+m3)*(l0+a1*X(1)+a2*X(3))-m3*ls3)*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*l0^2*(-(A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*cos(theta20)/l0^2-m2*ls2*a1/l0)*b1;
%J43 =b22*(l2*cos(theta20)+l0*sin(theta20))-a22*cos(theta20)+(2*(m2+m3)*a2*(l01+a11*X(1)+a1*X(2)+a21*X(3)+a2*X(4))+2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a21-2*m3*ls3*a21*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))-2*m3*ls3*(l01+a11*X(1)+a1*X(2)+a21*X(3)+a2*X(4))*b21)*(-m2*ls2*cos(theta20)/l0^2-(m2+m3)*a1/l0)+(m2+m3)*a2*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))^2/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*l0^2*(-(A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*cos(theta20)/l0^2-m2*ls2*a1/l0)+2*((m2+m3)*(l0+a1*X(1)+a2*X(3))-m3*ls3)*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*l0^2*(-(A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*cos(theta20)/l0^2-m2*ls2*a1/l0)*b21-((m2+m3)*(l0+a1*X(1)+a2*X(3))-m3*ls3)*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))^2/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)^2*l0^2*(-(A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*cos(theta20)/l0^2-m2*ls2*a1/l0)*(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a2-2*m3*ls3*a2)*(m2+m3)-((m2+m3)*(l0+a1*X(1)+a2*X(3))-m3*ls3)*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))^2/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a2-2*m3*ls3*a2)*cos(theta20)+k*r/(X(1)^2+X(3)^2)^(3/2)*X(3)/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*l0^2*(X(1)*((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*b1*b2-m2*ls2*(a1*b2+b1*a2)+(m2+m3)*a1*a2)-X(3)*((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*b1^2-2*m2*ls2*a1*b1+(m2+m3)*a1^2))-k*(1-r/(X(1)^2+X(3)^2)^(1/2))/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)^2*l0^2*(X(1)*((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*b1*b2-m2*ls2*(a1*b2+b1*a2)+(m2+m3)*a1*a2)-X(3)*((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*b1^2-2*m2*ls2*a1*b1+(m2+m3)*a1^2))*(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a2-2*m3*ls3*a2)*(m2+m3)+k*(1-r/(X(1)^2+X(3)^2)^(1/2))/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*l0^2*(X(1)*(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a2-2*m3*ls3*a2)*b1*b2-(A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*b1^2+2*m2*ls2*a1*b1-(m2+m3)*a1^2-X(3)*(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a2-2*m3*ls3*a2)*b1^2);
%J44 =2*b21*(l2*cos(theta20)+l0*sin(theta20))-2*a21*cos(theta20)+(2*(m2+m3)*(l0+a1*X(1)+a2*X(3))*a2-2*m3*ls3*a2*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))-2*m3*ls3*(l01+a11*X(1)+a1*X(2)+a21*X(3)+a2*X(4))*b2)*(-m2*ls2*cos(theta20)/l0^2-(m2+m3)*a1/l0)+2*((m2+m3)*(l0+a1*X(1)+a2*X(3))-m3*ls3)*(theta201+b11*X(1)+b1*X(2)+b21*X(3)+b2*X(4))/((A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*(m2+m3)-m2^2*ls2^2)*l0^2*(-(A+(m2+m3)*(l0+a1*X(1)+a2*X(3))^2-2*m3*ls3*(l0+a1*X(1)+a2*X(3)))*cos(theta20)/l0^2-m2*ls2*a1/l0)*b2;
Jac=[0, 1, 0, 0;
J21,J22,J23,J24;
0 , 0, 0, 1;
J41,J42,J43,J44];
%Variational equation
f(5:20)=Jac*Y;
%Output data must be a column vector
wystar 发表于 2014-7-1 11:30
你是在计算Lorenz方程的Lyapunov指数的基础上改了一下，后面再计算Jac的时候，你没有必要直接 对各个项进行求导计算Jac；只需要利用matlab，其中有求Jac矩阵简单的调用语句；在你自己计算Lyapunov之前，你先编辑调用Jac语句，先把Jac算出来，让后再代入到Lorenz方程中的计算Lyapunov指数的JAc，进行替换，就可以了 我就是这样算出来的 而且不会出现错误你试试看
wystar 发表于 2014-7-1 11:30
而且你的是自治系统好像 计算更简单， 不用将非自治变为自治系统
而且你的是自治系统好像 计算更简单， 不用将非自治变为自治系统
我的也是非自治系统啊。分母中含有x,y，对x求偏导之后还含有x。非自治系统一定要化成自治系统吗？怎么化啊？
你是在计算Lorenz方程的Lyapunov指数的基础上改了一下，后面再计算Jac的时候，你没有必要直接 对各个项进行求导计算Jac；只需要利用matlab，其中有求Jac矩阵简单的调用语句；在你自己计算Lyapunov之前，你先编辑调用Jac语句，先把Jac算出来，让后再代入到Lorenz方程中的计算Lyapunov指数的JAc，进行替换，就可以了 我就是这样算出来的 而且不会出现错误你试试看
的确是在Lorenz方程的Lyapunov指数的基础上改了一下。只需要利用matlab，其中有求Jac矩阵简单的调用语句？怎么调用啊。
而且你的是自治系统好像 计算更简单， 不用将非自治变为自治系统
补充上面方程的。这方面的知识了解的不是很多，多谢指导。
wystar 发表于 2014-7-1 15:16
的确是在Lorenz方程的Lyapunov指数的基础上改了一下。只需要利用matlab，其中有求Jac矩阵简单的调用语句 ...
MATLAB中jacobian是用来计算Jacobi矩阵的函数。
syms r l f
x=r*cos(l)*cos(f);
y=r*cos(l)*sin(f);
z=r*sin(l);
J=jacobian(,) 调用计算 非常简单
mumianhua 发表于 2014-7-1 15:25
MATLAB中jacobian是用来计算Jacobi矩阵的函数。
syms r l f
x=r*cos(l)*cos(f);
一定必须是符号计算哦别忘记了 不然计算不出来哦
wystar 发表于 2014-7-1 15:21
补充上面方程的。这方面的知识了解的不是很多，多谢指导。
MATLAB中jacobian是用来计算Jacobi矩阵的函数。
syms r l f
x=r*cos(l)*cos(f);
y=r*cos(l)*sin(f);
z=r*sin(l);
J=jacobian(,) 调用计算 非常简单 你先试试 计算一下Jac矩阵，然后再代入到你 的Lya计算中去
MATLAB中jacobian是用来计算Jacobi矩阵的函数。
syms r l f
x=r*cos(l)*cos(f);
y=r*cos(l)*sin(f);
z=r*sin(l);
J=jacobian(,) 调用计算 非常简单 你先试试 计算一下Jac矩阵，然后再代入到你 的Lya计算中去
好的，谢谢！我试一下。
你先试试 计算一下Jac矩阵，然后再代入到你 的Lya计算中去。
=lyapunov(4,@lorenz_ext,@ode45,0,0.5,50,,10);
雅克比是在matlab中直接调用的，求解还有问题。然后我将求解器改为ode45，求解时间设为50s，只是取几个值验证结果是否正确；结果求出来的Lyapunov指数为无穷大，而且求解的特别缓慢。
t=5.0000 NaN NaN NaN NaN
t=10.0000 NaN NaN NaN NaN
t=15.0000 NaN NaN NaN NaN
t=20.0000 NaN NaN NaN NaN
t=25.0000 NaN NaN NaN NaN
t=30.0000 NaN NaN NaN NaN
t=35.0000 NaN NaN NaN NaN
t=40.0000 NaN NaN NaN NaN
t=45.0000 NaN NaN NaN NaN
t=50.0000 NaN NaN NaN NaN
Elapsed time is 5803.389644 seconds.
而且你的是自治系统好像 计算更简单， 不用将非自治变为自治系统
是不是由于我没有将它变为自治系统的原因而导致求解结果无穷大。
我现在调整了步长和求解时间，结果是下面的两幅图，结果对吗？
=lyapunov(4,@lorenz_ext,@ode45,0,0.001,1,,10);
为什么当求解时间增大的时候，Lyapunov指数就变成了无穷大？