请高手帮忙看看我的程序有什么问题，谢谢！

yina_111

发现原来那个程序在求最大的LE指数时的迭代有点问题，进行了修改，但是运行的时候还是有问题，请高人，大侠帮帮忙，我现在急用这个程序，非常感谢！
Le=zeros(401,2);
i=0;f=[1 0; 0 1];
for q=0.4:0.01:1.4
    i=i+1;
    x=0.83;y=0.55;
    z=0;w=0;
    for j=1:100
        x=x+q*(2*x*(1-x)-0.6*x*y);
        y=y+q*(-0.5+0.6*x*y);
    end
    for k=1:200
        x=x+q*(2*x*(1-x)-0.6*x*y);
        y=y+q*(-0.5+0.6*x*y);
        J=[1+q*(2-4*x)-0.6*y q*0.6*x;
          q*0.6*y  1+q*(2-4*x)-0.6*y];
      f=f*J
        z=log(abs(max(eig(f))));
    end
    Le(i,1)=q;
    Le(i,2)=z/200;
end
plot(Le(:,1),Le(:,2));
hold on;
n=0.4:0.0001:0.47;
plot(n,0,'k');
xlabel('q');
ylabel('Le');

yina_111

我现在急用这个程序，请高人指点一下
非常感谢！！！

11m

我是新手，有个问题。J第二行第二列1+q*(2-4*x)-0.6*y是怎么得到的？难道不是对y求导吗？那应该是1+q*0.6*x吧。

11m

z=log(abs(max(eig(f))))根据你付的图，应该是z=log(max(abs(eig(f))))。

liliangbiao

如果我没有看错的话 应该是计算二维离散系统的LEs,我给你一个，以Henon为例：当然了 你还可以进一步修改你的程序，我觉得你的这个程序除了细节上的错误之外，没有别的毛病，思想很正确！

yina_111

原帖由 liliangbiao 于 2008-11-15 08:29 发表

                               
登录/注册后可看大图

如果我没有看错的话 应该是计算二维离散系统的LEs,我给你一个，以Henon为例：当然了 你还可以进一步修改你的程序，我觉得你的这个程序除了细节上的错误之外，没有别的毛病，思想很正确！

太感谢了，发到我的邮箱里吧
yina_111@163.com

万分感谢

yina_111

谢谢11m，给我指出的问题，但是主要问题不是在这个地方，
当运行程序后
出现了一个问题
??? Error using ==> eig
NaN or Inf prevents convergence.

Error in ==> Fjcly1 at 17
        z=log(max(abs(eig(f))));

我不清楚这什么原因导致的，即使写成z=log(max(abs(eig(f))))，出现的问题还是一样的

yina_111

您如果有时间，能帮我看看问题出现在什么地方么，万分感谢！！！

11m

原来的映射是
x(n+1)=x(n)+q*(2*x(n)(1-x(n))-0.6*x(n)*y(n))
y(n+1)=y(n)+q*(-0.5+0.6*x(n)*y(n))
对吗？
我计算了，不动点是（0.5，10/6）
jacobian应该是
J=[1+q*(2-4*x-0.6*y) -q*0.6*x;
          q*0.6*y  1+q*0.6*x]
我不理解你为什么要用 x=0.83;y=0.55当初值，我用0，0当初值也是发散，所以就用0.5，1.67当初值，程序如下
Le=zeros(401,2);
i=0;f=[1 0; 0 1];
for q=0.4:0.01:1.4
    i=i+1;
    x=0.5;y=1.667;
    z=0;w=0;
    for j=1:100
        x=x+q*(2*x*(1-x)-0.6*x*y);
        y=y+q*(-0.5+0.6*x*y);
    end
for k=1:200
        x=x+q*(2*x*(1-x)-0.6*x*y);
        y=y+q*(-0.5+0.6*x*y);
        J=[1+q*(2-4*x-0.6*y) -q*0.6*x;
          q*0.6*y  1+q*0.6*x];
        f=f*J;
z=log(abs(max(eig(f))));
    end
    Le(i,1)=q;
    Le(i,2)=z/200;
end
plot(Le(:,1),Le(:,2));
hold on;
n=0.4:0.0001:0.47;
plot(n,0,'k');
xlabel('q');
ylabel('Le');
得到了两条直线两条曲线。我也不明白线的意思。不过都是大于0。

11m

n=0.4:0.0001:0.47;
plot(n,0,'k');
这两行是什么意思？只画一条直线？

11m

%             x(n+1)=x(n)+q*(2*x(n)(1-x(n))-0.6*x(n)*y(n))
%             y(n+1)=y(n)+q*(-0.5+0.6*x(n)*y(n))
            
Le=zeros(101,2);
i=0;x=zeros(301,1);y=zeros(301,1);
for q=0.4:0.01:1.4
    i=i+1;
    x(1)=0.5;y(1)=1.67;
           for n=1:100
               x(n+1)=x(n)+q*(2*x(n)*(1-x(n))-0.6*x(n)*y(n));
               y(n+1)=y(n)+q*(-0.5+0.6*x(n)*y(n));
           end
           F=[1 0;0 1];
           for n=101:300
               x(n+1)=x(n)+q*(2*x(n)*(1-x(n))-0.6*x(n)*y(n));
               y(n+1)=y(n)+q*(-0.5+0.6*x(n)*y(n));
               J=[1+q*(2-4*x(n)-0.6*y(n)) -q*0.6*x(n);
                          q*0.6*y(n)  1+q*0.6*x(n)];
               F=F*J;
            end
     z=log(max(abs(eig(F))));
     Le(i,1)=q;
     Le(i,2)=z/200;
end
plot(Le(:,1),Le(:,2));
这样能画出来。
但是如果把初始值设为   x(1)=0;y(1)=0;则提示NaN。所以我想可能（0.5，10/6）是不稳定不动点。具体我也不清楚，也可能你的迭代式本身有错误。

11m

x=x+q*(2*x*(1-x)-0.6*x*y);
y=y+q*(-0.5+0.6*x*y);
这样写是不对滴。应该加n和n+1才对。

yina_111

首先非常感谢11m对我的帮助和指点，
x=x+q*(2*x*(1-x)-0.6*x*y);
y=y+q*(-0.5+0.6*x*y);

这么写没有错误，

加n和n+1也对

[ 本帖最后由 yina_111 于 2008-11-18 16:09 编辑 ]

yina_111

不知道你的这个程序里面的
  for n=101:300
               x(n+1)=x(n)+q*(2*x(n)*(1-x(n))-0.6*x(n)*y(n));
               y(n+1)=y(n)+q*(-0.5+0.6*x(n)*y(n));
               J=[1+q*(2-4*x(n)-0.6*y(n)) -q*0.6*x(n);
                          q*0.6*y(n)  1+q*0.6*x(n)];
               F=F*J;
            end
n怎么从101开始呢？

yina_111

我运行了一下这个程序，虽然画出了图
但是和我要得差别太大了，我觉得还是有问题