关于最大Lyapunov指数随参数变化的曲线绘制问题！

无水1324

看了上面的程序，没有看出来问题。等回学校之后调试一下再说

咕噜噜

说实话我没有看得很明白，因为觉得怪怪的，具体怎么也说不好
所以问清楚了好仔细想想

octopussheng

咕噜噜你好，首先先致意诚挚的谢意！

这个程序是仿照let程序中的一个demo——lorenz系统而编写的！不知道这个程序你有没有。

问题1：输出是[dx1;dx2;dx3;dx4]
问题2：dx1=x2,因为定义状态向量的时候是这样定义的：[x,x',y,y']
所以dx1=x2,dx2=......,dx3=x4,dx4=.......
问题3：Q = [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)];
%由于输入数据是一个行向量，因此需要将输入数据按照指定的顺序进行排列

octopussheng

我把demo中的范例程序也发上来，可以对照着看一下的。

function OUT=lorenzeq(t,X)
%LORENZEQ  Lorenz equation
%          (a 3rd-order continuous autonomous system):
%
%               dx/dt = SIGMA*(y - x)
%               dy/dt = RHO*x - y -x*z
%               dz/dt= x*y - BETA*z
%
%        In this demo, SIGMA = 16, RHO = 45.92, BETA = 4
%        Initial conditions: x(0) = 1, y(0) = 1, z(0) = 1;
%        Reference values: LE1 = 1.497, LE2 = 0.00, LE3 = -22.46, LD = 2.07
%
%        The reference values are from the following references:
%
%        [1] 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.
%
%        [2] Keith Briggs, "An Improved Method for Estimating Liapunov
%            Exponents of Chaotic Time Series," Phys. Lett. A, Vol. 151,
%            pp. 27-32, Nov. 1990.

%        by Steve Wai Kam SIU, Jun. 29, 1998.

%PARAMETERS
SIGMA = 16;
RHO = 45.92;
BETA = 4;

%Rearrange input data in desired format
%Note: the input data is a column vector
x=X(1);y=X(2);z=X(3);
Q=[X(4), X(7), X(10);
    X(5), X(8), X(11);
    X(6), X(9), X(12)];

%Lorenz equation
dx=SIGMA*(y-x);
dy=-x*z+RHO*x-y;
dz=x*y-BETA*z;

DX1=[dx;dy;dz];        %Output data

%Linearized system
J=[-SIGMA, SIGMA,     0;
     RHO-z,    -1,    -x;
         y,     x, -BETA];
  
%Variational equation   
F=J*Q;

%Output data must be a column vector
OUT=[DX1; F(:)];

咕噜噜

呵呵，不好意思有事才上来
按照你的答复和我的理解是一致的，可是感觉你的程序有那么点问题
Q矩阵是需要输入的数据，如果不是程序过程中求得的话最好直接写出来更合适
我帮你看看这个程序阿明天，希望能帮到你

[ 本帖最后由 咕噜噜 于 2007-7-8 21:22 编辑 ]

octopussheng

实在是太感谢了！呵呵！

期望早日解决这个问题啊！

咕噜噜

^_^，我把Q输入具体参数就能得到结果了，呵呵
其他还真是看不出问题啊

octopussheng

Q输入具体参数？？？？

从他的这个例子中来看，Q的数据不用管的啊？

不知道你得出来的结果如何啊？

咕噜噜

:@o 不用管？不会吧，那你的Q打算怎么办，没有定义啊你这样
%Note: the input data is a column vector
x=X(1);y=X(2);z=X(3);
Q=[X(4), X(7), X(10);
    X(5), X(8), X(11);
    X(6), X(9), X(12)];
看看这个，你的程序里面写的，数据要输入的
我随便输入的

咕噜噜

相关程序如下
function OUT = L_2( t,  X)
X=[1 5 2 6 5 2 3 5 4 9 5 4 5 12 75 85 6 39 65 20];
%%%%%%%%%%%%%%%%%%     X中前四个元素对应于两自由度非线性悬架系统中的4个状态变量，X是一个具有16个元素的向量
x1 = X(1);  x2 = X(2);  x3 = X(3);  x4 = X(4);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       基本参数
ms=337;
mt=50;
k0=-2316.4;
k1=12394;
k2=-73696;
k3=3170400;
c1=1228;
c2=1210;
k4=98000;
k5=1850000;
B10=k0/ms;
B11=k1/ms;
B12=k2/ms;
B13=k3/ms;
C11=c1/ms;
C12=c2/ms;
B20=k0/mt;
B21=k1/mt;
B22=k2/mt;
B23=k3/mt;
B24=k4/mt;
B25=k5/mt;
C21=c1/mt;
C22=c2/mt;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A=0.02;
omiga=6;
t=linspace(0,1,1000);
for i=1:1000
zr(i)=A*sin(omiga*t(i)); %%%%    激励
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dx1 = x2;
dx2 = -B10-B11*(x1-x3)-B12*(x1-x3).^2-B13*(x1-x3).^3-C11*(x2-x4)-C12*(x2-x4).^2;
dx3 = x4;
dx4 = B20+B21*(x1-x3)+B22*(x1-x3).^2+B23*(x1-x3).^3+C21*(x2-x4)+C22*(x2-x4).^2+B24*(x3-zr(i))+B25*(x3-zr(i)).^2;

Q = [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)];
J = [                0                            1                                    0                                    0;
    -B11-2*B12*(x1-x3)-3*B13*(x1-x3).^2  -C11-2*C12*(x2-x4)          B11+2*B12*(x1-x3)+3*B13*(x1-x3).^2              C11+2*C12*(x2-x4);
                     0                            0                                    0                                    1;
     B21+2*B22*(x1-x3)+3*B23*(x1-x3).^2   C21+2*C22*(x2-x4)  -B21-2*B22*(x1-x3)-3*B23*(x1-x3).^2+B24+2*B25*(x3-zr(i))  -C21-2*C22*(x2-x4)];
F = J*Q;
end
OUT = [dx1; dx2; dx3; dx4; F(:)];
不知道是不是你想要得结果，输入太长了，可是我觉得我把时间t循环了那，不知道对不对

[ 本帖最后由 咕噜噜 于 2007-7-9 10:27 编辑 ]

octopussheng

不知道你看上面那个lorenz系统的例子没有，我是看着let程序的帮助写程序的，你这样写我觉得好像有问题啊！

octopussheng

用let算上面的lorenz系统的LE时，Q矩阵中的各个参数都没有定义啊！
出来的结果和理论值很接近的啊！

Q矩阵应该是没有问题的吧

咕噜噜

原帖由 octopussheng 于 2007-7-9 11:02 发表

                               
登录/注册后可看大图

不知道你看上面那个lorenz系统的例子没有，我是看着let程序的帮助写程序的，你这样写我觉得好像有问题啊！

即便是那个lorenz系统的例子也是需要输入Q的阿，那个lorenz系统的例子你不输入Q就可以正常计算了吗？

octopussheng

对啊，不用输入Q的，直接就可以计算的，你可以试试

咕噜噜

我试过了，不可以，难道我电脑有问题