求教:关于最大lyapunov指数的问题

ldyw

我做出的最大lyapunov指数和分岔图老是对应不上，就是说单周期的时候，最大lyapunov指数也大于0，程序我用算例试过，也没问题。是不是和某个参数的选取有关啊，求高手指教！！！！！！

[ 本帖最后由 ldyw 于 2007-4-25 20:43 编辑 ]

hohoo

程序贴出来 看看

shenyongjun

原帖由 ldyw 于 2007-4-25 20:40 发表
单周期的时候，最大lyapunov指数也大于0，

那肯定是程序或者算法的问题！
据个人所知，目前流行的计算最大Lyapunov指数的方法，有些只能用于混沌系统，这样会得到比较准确的结果。相反，如果用这些方法计算非混沌的非线性系统，肯定会得到错误的结果。就像楼主所说的！

ldyw

模型的方程如下：
double fun(double t, double y[N],double d[N],double xy)
{
double  m1,m2,d1,d2;
double k11,k12,k13,k14,k22,k21,k23,k24,k31,k32,k33,k34,k41,k42,k43,k44;
double wan,E1,E2,E,kr,delta,a,w3,w4,q1,q2,m11,m22;
double f;
double w,w1,w2;
f=0.12;
w=xy;
d1=d2=0.08;
m1=m2=102;
a=1;
q1=q2=0.2;
w3=q1/2;
w4=q2/2;
wan=61360/a/a/a/2;
k11=21*wan;k12=k21=3*a*wan;k13=k31=9*wan;k14=k41=-3*a*wan;
k22=13*a*a*wan;k23=k32=3*a*wan;k24=k42=-a*a*wan;k33=15*wan;
k34=k43=-9*a*wan;k44=7*a*a*wan;         
w1=sqrt(k11/m1);
w2=sqrt(k33/m2);
delta=0.003;
E1=E2=0.3;


d[0]=y[8];
d[1]=y[9];
d[2]=y[10];
d[3]=y[11];
d[4]=y[12];
d[5]=y[13];
d[6]=y[14];
d[7]=y[15];

E=sqrt(y[2]*y[2]+y[6]*y[6]);
if(E>1||E==1)
kr=4*k11;
else
kr=0;



d[8]=-2*d1*w1*y[8]/w-k11*y[0]/m1/w/w-k12*y[1]/m1/w/w/delta-k13*y[2]/m1/w/w-k14*y[3]/m1/w/w/delta+E1*cos(t);
d[9]=2*y[13]-k21*y[0]*delta/m1/w/w/w3/w3-k22*y[1]/m1/w/w/w3/w3-k23*y[2]*delta/m1/w/w/w3/w3-k24*y[3]/m1/w/w/w3/w3;
d[10]=-2*d2*w2*y[10]/w-k31*y[0]/m2/w/w-k32*y[1]/m2/w/w/delta-k33*y[2]/m2/w/w-k34*y[3]/m2/w/w/delta+E1*cos(t)-(E-1)*kr*(y[2]-f*y[6])/E/m2/w/w;
d[11]=2*y[15]-k41*y[0]*delta/m2/w/w/w4/w4-k42*y[1]/m2/w/w/w4/w4-k43*y[2]*delta/m2/w/w/w4/w4-k44*y[3]/m2/w/w/w4/w4;
d[12]=-2*d1*w1*y[12]/w-k11*y[4]/m1/w/w+k12*y[5]/m1/w/w/delta-k13*y[6]/m1/w/w+k14*y[7]/m1/w/w/delta+E1*sin(t);
d[13]=-2*y[9]+k21*y[4]*delta/m1/w/w/w3/w3-k22*y[5]/m1/w/w/w3/w3+k23*y[6]*delta/m1/w/w/w3/w3-k24*y[7]/m1/w/w/w3/w3;
d[14]=-2*d2*w2*y[14]/w-k31*y[4]/m2/w/w+k32*y[5]/m2/w/w/delta-k33*y[6]/m2/w/w+k34*y[7]/m2/w/w/delta+E1*sin(t)-(E-1)*kr*(f*y[2]+y[6])/E/m2/w/w;
d[15]=-2*y[11]+k41*y[4]*delta/m2/w/w/w4/w4-k42*y[5]/m2/w/w/w4/w4+k43*y[6]*delta/m2/w/w/w4/w4-k44*y[7]/m2/w/w/w4/w4;
return(1);
}
lyapunov指数程序如下：
double  Lyapunov_index(double y[N], double h,double xy )
{
  int k;
  long i,j,mm,zhou_shu,per_shu,zhou_shu0,zhou_shu1;
  double d0=0.0000001, lnd0, d, sum_ln_di=0.0, t;
  double z[N], y0[N], z0[N],hx,ht,dd=0;

  zhou_shu1=50000;
  for(k=0;k<N;K++)
      y0[k]=y[k];

   for(k=0; k<N; k++)
   {                                       
    z[k]=z0[k]=y[k]+d0;        
  }

   dd=0;                                       
   for(k=0;k<N;K++)
    dd=dd+fabs((y[k]-z[k])*(y[k]-z[k]));
   d0=sqrt(dd);
   lnd0=log(d0);
   
  t=0;
  for(j=0;j<ZHOU_SHU1;J++)
   {
    t=0;
    t=t+j*h;                        
    rkutta(t,h,y,xy);        
    rkutta(t,h,z,xy);
   
    dd=0;
    for(k=0;k<N;K++)
     {                                                   
    dd=dd+fabs((y[k]-z[k])*(y[k]-z[k]));                     
    }
    d=sqrt(dd);        
    sum_ln_di+= log(d);         
    for (k=0; k<N; d;
      
  }                          


  d = sum_ln_di/(zhou_shu1);

  d = (d - lnd0)/h;

   return(d);
}

请帮忙看看那里有问题，是方程那里有问题还是程序，谢谢。

我作出的结果： 分岔图从单周期进入多周期过程，lyapunov指数没有在分岔点对应正负的变化，
在单周期时就存在大于0的情况了。

yejet

c语言？看不懂