lyapunov_small.c

求取lyapunov指数的小数据量方法

#include <math.h>
#include "mex.h"
#include "stdio.h"
#include "stdlib.h"
#include "matrix.h"


int ABS(int x)
{
    int y;
    if (x>=0)
        y = x;
    else
        y = -x;
    return y;
}

int MIN(int x,int y)
{
    
    if (x<=y)
        return x;
    else
        return y;
}


// 计算数组最大值
int MAX_VECTOR(int *p_vector,
                  int len_vector)
{
    int i;    
    int max_value = *p_vector;
    for (i=0; i<len_vector; i++)
    {   if (*(p_vector+i)>max_value)
	{
		max_value = *(p_vector+i);
	}
    }
    return max_value;
}



double ** TwoArrayAlloc(int r,int c)
{  //两维内存动态分配函数
	double  *x, **y;
	int n;
	
	x=(double *)calloc(r*c,sizeof(double));
	if(x==NULL)
    {
		printf("\nFAilue in memory applying.");
		//AfxMessageBox("Failure in memory applying.");  
		exit(0);
    }
	y=(double **)calloc(r,sizeof(double *));
	for(n=0;n<=r-1;++n)
		y[n]=&x[c*n];
	return(y);
}

void TwoArrayFree(double **x)//2 维内存释放函数
{
	free(x[0]);
	free(x);
}



//求最短距离及向量
void min_dist(double **pdist,int n,int *place,double *pmindist)
{  int k,j;    
   double min_value ;
   
   for (k=0;k<n;k++)
   {
	   if (pdist[k][1]!=0) 
	   {min_value=pdist[k][1];j=k+1;
	   break;}
   }

  for (k=0; k<n; k++)
  {   
	  if ((pdist[k][1]!=0)&&(pdist[k][1]<min_value))
	  {	 min_value = pdist[k][1]; j=k+1; }
  }
  
  *place=j;
  *pmindist=min_value;

}

//---------------------------------------------------------------------------
void mexFunction (int nlhs, mxArray *plhs[],			// 输出参数个数，及输出参数数组
			 int nrhs, const mxArray *prhs[])	// 输入参数个数，及输入参数数组
{
    int m,tau,N,P;
	int M;
	int i,j,jj,k,h,idx_j=0;
    int q,max_i,*max_num,max_q;  
    double **dd,**d,*y,**Y,min_d,y_s,sum_d;
	
 	double *pdata;
	
	if (nrhs!=4) mexErrMsgTxt("需要4个参数!");  //检查输入参数的个数
	
    // 取得输入参数
    pdata = mxGetPr(prhs[0]);      // 时间序列（列向量）      
    tau = (int) *mxGetPr(prhs[1]);            // 最小嵌入维数 
	m = (int) *mxGetPr(prhs[2]);            // 最大嵌入维数 
    P = (int) *mxGetPr(prhs[3]);      // 时间延迟   
	N = mxGetM(prhs[0]);              // 序列长度
	
      // 为输出变量分配内存空间
	plhs[0]=mxCreateDoubleMatrix(N,1,mxREAL);
    
	// 取得输出参数指针
	y = mxGetPr(plhs[0]);
	 	
    M=N-(m-1)*tau;
    d=TwoArrayAlloc(M,M);
	max_num=(int*)malloc(M*sizeof(int));
    Y=TwoArrayAlloc(m,M);
	for (j=0;j<M;j++)  //相空间重构
	{	for (i=0;i<m;i++) 
			Y[i][j]=*(pdata+i*tau+j);
		
	}
			
	
    for (j=1;j<=M;j++) 
	{	k=0;
	    min_d=0.0;
        dd=TwoArrayAlloc(M,2);
	
        for (jj=1;jj<=M;jj++)         //寻找相空间中每个点的最近距离点，并记下该点下标
		{  
			if (ABS(j-jj)>P)
			{  dd[k][0]=jj;
		         sum_d = 0.0;
				 for (h=1;h<=m;h++)     
				 {   sum_d = sum_d + (Y[h-1][j-1]-Y[h-1][jj-1])*(Y[h-1][j-1]-Y[h-1][jj-1]);   }
				 dd[k][1]=sum_d;
				 k=k+1;
			}
		}
	
		min_dist(dd,k,&idx_j,&min_d);
			
		TwoArrayFree(dd);
		max_i=MIN((M-j),(M-idx_j));           //计算点j的最大演化时间步长i
		max_num[j-1]=max_i;
	
		for (i=1;i<=max_i;i++)              //计算点j与其最近邻点在i个离散步后的距离
		
		{
			
			sum_d = 0.0;
			for (h=1;h<=m;h++)     
			{   sum_d = sum_d + (Y[h-1][j-1+i]-Y[h-1][idx_j-1+i])*(Y[h-1][j-1+i]-Y[h-1][idx_j-1+i]);   }
			d[i-1][j-1]=sum_d;
		}
	
   }
	
	//对每个演化时间步长i，求所有的j的lnd(i,j)平均
	max_q=MAX_VECTOR(max_num,M);
        
	//y=(double*)malloc(max_q*sizeof(double));
	for (i=1;i<=max_q;i++) 
	{    q=0;
         y_s=0.0;
      for (j=1;j<=M;j++) 
	  { 
		  if (d[i-1][j-1]!=0)
		  {   q=q+1;
		     y_s=y_s+log(d[i-1][j-1]);
		  }
	  }
			*(y+i-1)=y_s/q;
	}
		
   TwoArrayFree(d);
   
    

}