最小二乘法

  • 格式:doc
  • 大小:34.50 KB
  • 文档页数:4

#include

#include

void main()

{

int i,j,k;

float m=0,n=0;

float u[250],A[248][3],B[3][248],a[3][1],b[1][3];/*u为输入,A相当于系数矩阵,B是A的转置,a相当于递推式中的hui,b是a的转置*/

float y[250]={10,9}; /*y为输出*/

float O[3][1]={{0},{0},{0}};

float O1[3][1]={{0},{0},{0}}; /*中间变量*/

float Q[3][1]={{0},{0},{0}}; /*Q为P与a的乘积*/

float K[3][1]={{0},{0},{0}};

float M[1][3]={0,0,0}; /*M为b与P的乘积*/

float P[3][3]={{100000},{0,100000},{0,0,100000}};

float P1[3][3]={{0},{0},{0}}; /*中间变量*/

float P2[3][3]={{0},{0},{0}}; /*P2为Q与M的乘积*/

float e,f,g; /*差分方程系数*/

printf("input e,f,g:");

scanf("%f,%f,%f",&e,&f,&g);

/*产生输入u*/

for (i=0;i<=249;i++){

u[i]=1;

}

/*产生输出y*/

for (i=2;i<=249;i++){

y[i]=-1*e*y[i-1]-1*f*y[i-2]+1*g*u[i];

}

/*形成矩阵A*/

for (i=0;i<=247;i++){

A[i][0]=-y[i+1];

A[i][1]=-y[i];

A[i][2]=u[i+2];

}

/*B为A的转置矩阵*/

for (i=0;i<=2;i++)

for(j=0;j<=247;j++)

B[i][j]=A[j][i];

/*输出系数矩阵A*/

for (i=0;i<=247;i++){

for (j=0;j<=2;j++)

printf("%10.4f",A[i][j]);

printf("\n");

}

/*输出矩阵B,B为A的转置*/

for (i=0;i<=2;i++){

for (j=0;j<=247;j++)

printf("%10.4f",B[i][j]);

printf("\n");

}

/*开始递推求e,f,g*/

for (i=0;i<=247;i++){

/*每循环一次取一组实测数据*/

for (j=0;j<=2;j++){

a[j][0]=B[j][i]; /*相当于递推公式中的hui*/

b[0][j]=A[i][j];

}

for (j=0;j<=2;j++){

for (k=0;k<=2;k++){

Q[j][0]+=P[j][k]*a[k][0]; /*相当于求递推公式中的P乘以hui*/

M[0][j]+=b[0][k]*P[k][j]; /*相当于求递推公式中的hui的转置乘以P*/

}

}

/*for (j=0;j<=2;j++)

printf("%15.4f",Q[j][0]);

printf("\n");*/

for (j=0;j<=2;j++){

m+=b[0][j]*Q[j][0]; /*相当于求递推公式中的hui的转置乘以P乘以hui*/

n+=b[0][j]*O[j][0]; /*相当于求递推公式中的hui的装置乘以theta*/

}

m=m+1;

/*printf("%15.4f\n",m);

printf("%15.4f\n",n);*/

for (j=0;j<=2;j++){

K[j][0]=Q[j][0]/m; /*相当于求递推公式中的K(n+1)*/

/*printf("K[j][0]=%15.4f",K[j][0]);*/

}

for (j=0;j<=2;j++){

O1[j][0]=O[j][0]+K[j][0]*(y[i+2]-n);

O[j][0]=O1[j][0]; /*相当于求递推公式中的theta(n+1)*/

/*printf("%15.4f",O[j][0]);*/

}

for (j=0;j<=2;j++){

for (k=0;k<=2;k++){

P2[j][k]+=Q[j][0]*M[0][k];

P1[j][k]=P[j][k]-P2[j][k]/m;

P[j][k] =P1[j][k]; /*相当于求递推公式中的P(n+1)*/

/*printf("P[j][k]=%15.4f",P[j][k]);*/

}

}

for (j=0;j<=2;j++)

for (k=0;k<=2;k++)

P2[j][k]=0;

for (j=0;j<=2;j++){

Q[j][0]=0;

M[0][j]=0;

}

m=0;

if (y[i+2]-n<0.000001) /*循环结束判断条件*/

break;

n=0;

}

printf("%d\n",i);

for (j=0;j<=2;j++)

printf("%f\n",O[j][0]);

system("pause");

}