用C写的可以实现求解满秩线性方程组以及最小二乘曲线拟合的函数.rar

#define max1 10 ///拟合多项式的最高幂数或线性方程组的最大阶数 #define max2 1000///最小二乘拟合时的最大离散点数 ////最小二乘曲线拟合 BOOL CDataProcess::LeastSquare(double x[], double y[], int m, int n, double a[]) { //x,y -- X,Y两轴的坐标 //m 多项式的最高幂 //n 离散点数 //a 多项式的系数 int i,j,k,r; double f[MaxExp*2]; double A[MaxExp][MaxExp]; double Y[MaxExp]; double s=0; double max3,temp,xt; if(m>MaxExp) { printf("多项式幂数不能大于%d",MaxExp); exit(0); } if(n>MaxFitPt) { printf("输入点数不能大于%d",MaxFitPt); exit(0); } if(m>n) { m=n;//多项式的最高幂不能超过离散点数 } for(k=0;k<=2*m;k++)//计算x的次幂 { f[k]=0; for(i=0;i<n;i++) { s=pow(x[i],k); f[k]=f[k]+s; } } for(k=0;k<=m;k++)//计算Y { Y[k]=0; for(i=0;i<n;i++) { s=pow(x[i],k); Y[k]=Y[k]+y[i]*s; } } for(k=0;k<=m;k++)//把f数组赋给A，做正规方程组的系数矩阵 { for(j=0;j<=m;j++) { A[k][j]=f[k+j]; } } ///用高斯消去法解正规方程组 for(k=0;k<=m;k++) { r=k; max3=fabs(A[k][k]); for( i=k;i<=m;i++) { if(fabs(A[i][k])>MaxExp) { max3=fabs(A[i][k]); r=i; } } if(r!=k)//如果r和k不相等就说明需要进行两行交换 { for(int j=k;j<=m;j++) { temp=A[k][j];//下面三句是交换A数组 A[k][j]=A[r][j]; A[r][j]=temp; } temp=Y[k];//下面三句是交换Y数组 Y[k]=Y[r]; Y[r]=temp; } for(i=k+1;i<=m;i++) { xt=A[i][k]/A[k][k];//交换后置零 Y[i]=Y[i]-xt*Y[k]; for(j=k+1;j<=m;j++) A[i][j]=A[i][j]-xt*A[k][j]; A[i][k]=0; } } a[m]=Y[m]/A[m][m];//求的是最后一个x值 for (i=m-1;i>=0;i--)//求前面的几个x值 { temp=0; for (j=i+1;j<=m;j++) { temp=temp+A[i][j]*a[j]; } a[i]=(Y[i]-temp)/A[i][i]; } return TRUE; } BOOL CDataProcess::Gauss(double A[][MaxExp], double Y[], int m, double x[]) { double max3,temp,xt; int r,i,j; if(m>MaxExp) { printf("矩阵阶数不能大于%d",MaxExp); exit(0); } for(int k=0;k<=m;k++) { r=k; max3=fabs(A[k][k]); for( i=k;i<=m;i++) { if(fabs(A[i][k])>MaxExp) { max3=fabs(A[i][k]); r=i; } } if(r!=k)//如果r和k不相等就说明需要进行两行交换 { for(int j=k;j<=m;j++) { temp=A[k][j];//下面三句是交换A数组 A[k][j]=A[r][j]; A[r][j]=temp; } temp=Y[k];//下面三句是交换Y数组 Y[k]=Y[r]; Y[r]=temp; } for(i=k+1;i<=m;i++) { xt=A[i][k]/A[k][k];//交换后置零 Y[i]=Y[i]-xt*Y[k]; for(j=k+1;j<=m;j++) A[i][j]=A[i][j]-xt*A[k][j]; A[i][k]=0; } } x[m]=Y[m]/A[m][m];//求的是最后一个x值 for (i=m-1;i>=0;i--)//求前面的几个x值 { temp=0; for (j=i+1;j<=m;j++) { temp=temp+A[i][j]*x[j]; } x[i]=(Y[i]-temp)/A[i][i]; } return TRUE; }