bianshi.c

用c语言编程实现系统辨识

#include<stdio.h>
#include<malloc.h>

double comput_D(double *q, int c);                                               /*求矩阵的行列式，p为地址，c为行*/
double * Turn(double * q ,int row ,int line );                                   /*求矩阵转置*/
double * multiply(double * M1,int row1,int line1,double * M2,int row2,int line2);/*求矩阵的乘积*/
double Creat_M(double *p,int m,int n, int k);                                    /*求代数余子式*/
double * Qiuni(double * q,int row );                                             /*求矩阵的逆*/

void main()                                                            
{
	int i,j;
	int n,m;                                                                    /*n系统的阶次,设定起始阶次为2,m为p矩阵的总数据*/
    double sum;                                                                 /*矩阵的行列式值*/
	double * start;                                                             /*起始地址*/   
	double * p,* Y,* E,*A;                                                      /*p为输入u(0),u(1),u(2)...u(3n-2)和输出数据y(0),y(1),y(2)....y(3n-1)所构成的2n×2n矩阵 ,
						                                                          Y为[y(n),y(n+1)....y(3n+1)]的转置, E为残差,A为参数矩阵*/
	/********************************/
	/*初始化输入输出并赋值*/
	/********************************/
	
	double y[41]={1,0.5,-0.435,1.1649,4.4645,6.667,5.7681,1.1378,-7.4876,-16.8768,
		-16.6648,-11.9831,-2.253,5.4578,7.6442,7.4855,-2.1909,-11.5583,-14.1677,-11.7752,
        -5.9268,0.5166,7.1186,10.8402,7.5873,1.172,-4.1787,-6.6867,-2.8530,4.0195,
		6.3653,5.4713,6.4315,4.091,0.3473,0.9958,3.6766,4.9766,0.7394,-10.1271,-15.1661};
 
	double u[41]={-0.3,0.2,0.99,0.16,0.68,-0.87,0.22,-3.1,0.76,-0.75,
		-0.24,1.18,-1.92,1.81,-1.75,-1.49,0.28,-0.51,-0.08,-0.28,
		0.2,0.9,-0.28,-0.63,0.31,-0.47,0.55,1.05,0.25,-1.06,
		1.4,0.34,-0.76,0.84,0.98,0.19,0.19,-2.18,-1.39,1.36,-0.1};


    /**************开辟空间*****************/
	p=(double *)calloc(40*40,sizeof(double));
	Y=(double *)calloc(40*40,sizeof(double));
	E=(double *)calloc(40*40,sizeof(double));
	A=(double *)calloc(40*40,sizeof(double));

    start=p;
	/*开始循环，设本系统的初始阶次为2阶*/	
    for(n=1;n<=2;n++)
	{
		m=2*n;	
		/*************p矩阵赋值********************/		
		/********先给矩阵赋输出数据y（t）**********/                
		for(i=0;i<2*n;i++)                 		
		{
			for(j=0;j<n;j++)			
			{
				*(start+i*2*n+j)=y[j+i];
			}
		}    
		/********给矩阵赋值输入数据u（t）*********/	
		for(i=0;i<2*n;i++)
		{
			for(j=n;j<2*n;j++)
			{
				*(start+i*2*n+j)=u[j-n+i];
			}
		}



		for(i=0;i<2*n;i++)
		{
			for(j=0;j<2*n;j++)
				printf("%f",*(start+i*2*n+j));
			printf("\n");
		}
		sum=comput_D(p,m);  
		/*由递退方法求你n可知，判断sum是否为0，如果sum不为0，矩阵p为满秩，则n即为所要求的阶次，否则继续循环*/	    
		if(sum!=0)
		{
			printf("\nThe rank of this systerm is:%d",n);
			printf("\nThe value of Matrix is :%f:",sum);
		}
	}
    free(p);
	free(Y);
	free(E);
	free(A);
}


/*****************************************/
/*功能：求矩阵 c X c 的行列式*/
/*入口参数：矩阵首地址 q；矩阵行数 c*/
/*返回值：矩阵的行列式值*/
/*****************************************/
double comput_D(double *q,int c)  
{
     int i,j,m;         /*i--row; j--column*/
     int lop=0;
     double result=0;
     double mid=1;
   
     if (c!=1)
     {
         lop=(c==2)?1:c;     /*控制求和循环次数，若为2阶，则循环1次，否则为c次*/

         for(m=0;m<lop;m++)
         {
             mid=1;          /*顺序求和*/
             for(i=0,j=m;i<c;i++,j++)
                 mid = mid * ( *(q+i*c+j%c) );
             result+=mid;
         }

         for(m=0;m<lop;m++)
         {                       
             mid=1;          /*逆序相减*/
             for(i=0,j=c-1-m+c; i<c; i++,j--)
                 mid=mid * ( *(q+i*c+j%c));
             result-=mid;
         }
	 }
     else result=*q;
     return(result);
}

/*****************************************/
/*功能：求矩阵 row X line的转置*/
/*入口参数：矩阵首地址 q；矩阵行数 row,矩阵的列数line*/
/*返回值：矩阵的行列式值*/
/*****************************************/
double * Turn(double * q ,int row ,int line )
{ 
	int i,j;
	double * t;

	t=(double *)calloc(2*row*line,sizeof(double));	
	for(i=0;i<row;i++)
        for(j=0;j<line;j++)        
		  {
		     *(t+j*row+i)=*(q+i*line+j);

		  }
	return(t);
    free(t);
}
/*****************************************/
/*功能：求矩阵 A:row1行，line1列和矩阵B：row2行，line2列的乘积*/
/*入口参数：矩阵A首地址M1；矩阵B的首地址M2；A的行数row1,A的列数line1，B的行数row2,B的列数line2*/
/*返回值：矩阵的乘积值*/
/*****************************************/
double * multiply(double * M1,int row1,int line1,double * M2,int row2,int line2)
{
	double * C;
	int i,j,k,num3;
  	
	num3=2 * row1 * line2;
	C = (double *)calloc(num3, sizeof(double));
	for(i=0;i<row1;i++)          
	        for(j=0;j<line2;j++)
		       *(C+i*line2+j)=0;	
	for(i=0;i<row1;i++)
		{
			for(j=0;j<line2;j++)
			{
				for(k=0;k<line1;k++)
					*(C+i*line2+j)+=(*(	M1+i*line1+k))*(*(M2+k*line2+j));
			}
		}   
    return(C);
	free(C);
}
/*****************************************/
/*功能：求任何一个矩阵的逆*/
/*入口参数：矩阵首地址 q；矩阵行数row*/
/*返回值：方阵的逆*/
/*****************************************/
double * Qiuni(double * q,int row )
{
     double *p;                                      /*定义数组首地址指针变量*/
     int num;                                        /*定义矩阵行数row及矩阵元素个数*/
     int i,j;
     double determ;                                  /*定义矩阵的行列式*/
	 
     num=2 * row * row;
     p = (double *)calloc(num, sizeof(double));     /*分配内存单元*/
     determ=comput_D(q,row);                        /*求整个矩阵的行列式*/
     p=q + row * row;
     if (determ != 0)
     {
         for (i=0;i<row; i++)                      /*求逆矩阵*/
             for (j=0; j<row; j++)
                    *(p+j*row+i)=   Creat_M(q,i,j,row)/determ;   
            
         printf("The determinant is %G\n",determ);

         p=q + row * row;
         printf("\nThe inverse matrix is:\n");
         return(p);                               /*返回该矩阵*/
     }
     else
         return(0);
		 free( p );
}

/******************************************************/
/*功能：求k×k矩阵中元素A(mn)的代数余子式*/
/*入口参数：k×k矩阵首地址；元素A的下标m,n; 矩阵行数 k*/
/*返回值： k×k矩阵中元素A(mn)的代数余子式*/
/*****************************************************/
double Creat_M(double *p,int m,int n,int k)
{
     int len;
     int i,j;
     double mid_result=0;
     int quo=1;
     double *p_creat,*p_mid;

     len=(k-1)*(k-1);
     p_creat = (double *)calloc(len, sizeof(double));     
     p_mid=p_creat;
     for(i=0;i<k;i++)
         for(j=0;j<k;j++)
         {
             if (i!=m && j!=n)
                 *p_mid++ =* (p+i*k+j);            
         }
     quo = (m + n) %2==0 ? 1:-1;
     mid_result = (double ) quo * comput_D(p_creat,k-1);
     free(p_creat);
     return(mid_result);
}