small_data_sets_lyapunov1.c

小数据量法求Lyapunov指数_龙格库塔

#include "mex.h"
#include "stdlib.h"
#include "math.h"

long int ABS_longint(long int x)
{
  if(x<0)
    return (-x);
  else
    return (x);
}
double ABS_double(double x)
{
  if(x<0)
    return (-x);
  else
    return (x);
}
void small_data_sets_Lyapunov1(double *data,double *Length_data,double *m,double *tau,double *P,double *delt_t,double *LY_length,double *LY)
{
  long int i,j,Count,j1,j2;
  long int Length_data1,m1,tau1,P1;
  double d_j_min,d_s,S_j_i,d_j_i,delt_t1;
  long int max_i,sum_max_i,largest_max_i,index_j;
  long int sum_former,former;
  long int row,col;
  double *Q,*Q_j_jj;
  long int *Q_max_i,*Q_index_j;
  
  Length_data1=(long int)(*Length_data);
  m1=(long int)(*m);
  tau1=(long int)(*tau);
  P1=(long int)(*P);
  delt_t1=*delt_t;
  row=Length_data1-(m1-1)*tau1;
  col=m1;
  
  Count=0;
  Q=(double *)calloc(row*col,sizeof(double));
  Q_max_i=(long int *)calloc(row,sizeof(long int));
  Q_index_j=(long int *)calloc(row,sizeof(long int));
  for(i=0;i<row;i++)
    for(j=0;j<col;j++)
      Q[Count++]=data[i+j*tau1]; //相空间重构
  
  sum_max_i=0;
  largest_max_i=0;
  for(j1=0;j1<row;j1++)
  {
    d_j_min=10000;//寻找j1与j2的最短距离
    for(j2=0;j2<row;j2++)
    {
      d_s=0;//记录两点间的距离
      if(ABS_longint(j1-j2)>P1)
      {
         for(i=0;i<col;i++)
           d_s=d_s+ABS_double(Q[j1*col+i]-Q[j2*col+i]);
         if(d_s<d_j_min)
         {
            d_j_min=d_s;
            index_j=j2;
          }
      }   
    }
    Q_index_j[j1]=index_j;
    //printf("%f", Q_index_j[j1]);
    if((row-index_j)<(row-j1))
      max_i=row-index_j;
    else
      max_i=row-j1;
    Q_max_i[j1]=max_i;
    //printf("%f", Q_max_i[j1]);
    sum_max_i=sum_max_i+max_i;
    //printf("%d", sum_max_i);
    if(largest_max_i<max_i)
       largest_max_i=max_i;
  }
  
  Count=0;
  Q_j_jj=(double *)calloc(sum_max_i,sizeof(double));
  for(j1=0;j1<row;j1++)
  {
    for(i=0;i<Q_max_i[j1];i++)
    {
      d_j_i=0;
      for(j=0;j<col;j++)
        d_j_i=d_j_i+(Q[(j1+i)*col+j]-Q[(Q_index_j[j1]+i)*col+j])*(Q[(j1+i)*col+j]-Q[(Q_index_j[j1]+i)*col+j]);
      Q_j_jj[Count]=sqrt(d_j_i);
      //printf("%f",Q_j_jj[Count]);
      Count=Count+1;
    }
  }
  free(Q);
  
  for(i=0;i<largest_max_i;i++)
  {
    S_j_i=0;
	sum_former=0;
	Count=0;
	for(j=0;j<row;j++)
	{
	  if(j==0)
	    former=0;
	  else
	    former=Q_max_i[j-1];
	  sum_former=sum_former+former;
	  if(i<Q_max_i[j])
	  {
	    if(Q_j_jj[sum_former+i]!=0)
	    {
	      S_j_i=S_j_i+log(Q_j_jj[sum_former+i]);
	      Count=Count+1;
	    }
	  }
	}
	LY[i]=S_j_i/(Count*delt_t1);
  }
  *LY_length=largest_max_i;
  free(Q_max_i);
  free(Q_index_j);
  free(Q_j_jj);
}
void mexFunction(
     int nlhs,
     mxArray *plhs[],
     int nrhs,
     const mxArray *prhs[]
)
{
	double *x1,*x2,*x3,*x4,*x5,*x6,*y1,*y2;
    
	int y1_mrows,y1_ncols;
	int y2_mrows,y2_ncols;

	//检查参量的个数
	if (nrhs!=6)
			mexErrMsgTxt("6 inputs required");
  
    x1=mxGetPr(prhs[0]);
	x2=mxGetPr(prhs[1]);
	x3=mxGetPr(prhs[2]);
	x4=mxGetPr(prhs[3]);
	x5=mxGetPr(prhs[4]);
	x6=mxGetPr(prhs[5]);
	
	  
  y1_mrows=1;
  y1_ncols=1;
  
  y2_mrows=1;
  y2_ncols=(*x2)-((*x3)-1)*(*x4);
  
      
  plhs[0]=mxCreateDoubleMatrix(y1_mrows,y1_ncols,mxREAL);
  plhs[1]=mxCreateDoubleMatrix(y2_mrows,y2_ncols,mxREAL);
      
	y1=mxGetPr(plhs[0]);
	y2=mxGetPr(plhs[1]);
  small_data_sets_Lyapunov1(x1,x2,x3,x4,x5,x6,y1,y2);	   			
}