最小二乘法
- 格式: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");
}