多点最小二乘法平面方程拟合计算

1、-精选文档 -平面方程拟合计算平面方程的一般表达式为：AxByCzD0,( C0)zAxByDCCC记： a0A , a1B , a2DCCC则： za0 xa1 ya2平面方程拟合 :对于一系列的n 个点 ( n3)：( xi , yi , zi ), i 0,1, , n1要用点 ( xi , yi, zi ), i0,1, , n1拟合计算上述平面方程，则使：n1z 2Sa0 x a1 y a2i0最小。要使得 S 最小，应满足：S0, k 0,1,2ak2(a0 xia1 yia 2zi ) xi0即：2(a0 xia1 yia2zi ) yi02(a0 xia1 yia2zi )0a

2、0xi 2a1xi yia2xixi zi有， a0xi yia1yi 2a2yiyi zia0xia1yia2 nzixi 2xi yixia0xi zixi yi2yia1yi zi或，yixiyina2zi解上述线形方程组，得：a 0 , a1, a2可编辑-精选文档 -即： za0 xa1 ya2其程序代码如下：#include stdafx.h#include #include #include #define MAX 10void Inverse(double *matrix1,double *matrix2,int n,double d);double Determinant(d

3、ouble* matrix,int n);double AlCo(double* matrix,int jie,int row,int column);double Cofactor(double* matrix,int jie,int row,int column);int _tmain(int argc, _TCHAR* argv)double array123,Y3;double A,B,C;A = B = C = 0.0;ZeroMemory(array,sizeof(array);ZeroMemory(Y,sizeof(Y);for (int i = 0;i 12;i+)for (i

4、nt j = 0;j 3;j+)可编辑-精选文档 -arrayij = (double)rand();for (int i = 0; i 12;i+)arrayi0 = 1.0;/ 设计了 12 个最简单的数据点，x = 1平面上的点，double *Matrix3,*IMatrix3;for (int i = 0;i 3;i+)Matrixi= new double3;IMatrixi = new double3;for (int i = 0;i 3;i+)for (int j = 0;j 3;j+)*(Matrixi + j) = 0.0;for (int j = 0;j 3;j+)可编辑

5、-精选文档 -for (int i = 0;i 12;i+)*(Matrix0 + j) += arrayi0*arrayij;*(Matrix1 + j) += arrayi1*arrayij;*(Matrix2 + j) += arrayi2*arrayij;Yj -= arrayij;double d = Determinant(Matrix,3);if (abs(d) 0.0001)printf(n矩阵奇异 );getchar();return -1;Inverse(Matrix,IMatrix,3,d);for (int i = 0;i 3;i+)A += *(IMatrix0 +

6、i)*Yi;B += *(IMatrix1 + i)*Yi;C += *(IMatrix2 + i)*Yi;可编辑-精选文档 -printf(n A = %5.3f,B = %5.3f,C= %5.3f,A,B,C);for (int i = 0;i 3;i+)delete Matrixi;delete IMatrixi;getchar();return 0;void Inverse(double *matrix1,double *matrix2,int n,double d)int i,j;for(i=0;in;i+)matrix2i=(double *)malloc(n*sizeof(do

7、uble);for(i=0;in;i+)for(j=0;jn;j+)*(matrix2j+i)=(AlCo(matrix1,n,i,j)/d);double Determinant(double* matrix,int n)可编辑-精选文档 -double result=0,temp;int i;if(n=1)result=(*matrix0);elsefor(i=0;in;i+)temp=AlCo(matrix,n,n-1,i);result+=(*(matrixn-1+i)*temp;return result;double AlCo(double* matrix,int jie,int

8、row,int column)double result;if(row+column)%2 = 0)result = Cofactor(matrix,jie,row,column);else result=(-1)*Cofactor(matrix,jie,row,column);return result;可编辑-精选文档 -double Cofactor(double* matrix,int jie,int row,int column)double result;int i,j;double* smallmatrMAX-1;for(i=0;ijie-1;i+)smallmatri= new doublejie - 1;for(i=0;irow;i+)for(j=0;jcolumn;j+)*(smallmatri+j)=*(matrixi+j);for(i=row;ijie-1;i+)for(j=0;jcolumn;j+)*(smallmatri+j)=*(matrixi+1+j);for(i=0;irow;i+)for( j=column;jjie-1;j+)*(smallmatri+j)=*(matrixi+j+1);for(i=row;ijie-1;i+)for(j=column;jjie-