leastsquare.h

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

#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;
}