超定方程组的最小二乘解.txt

在生产和实验中

#include<iostream.h>
#include<iomanip.h>
#include<math.h>
#define m 4//n<m
#define n 2
void line_fun(double a[n][n],double b[n],double p[n])//列主高斯消元求解线性方程组
{
	double x[n];//线性方程组的解  
	double mm[n];//辅助乘数
	double T[n],BT=0;//存储换行临时变量
	double temp=0;
	double t;//最大列主元
	int tap1=0,tap2=0;//最大列主元下标
	int i,j,s;
	for(s=0;s<n;s++)
	{
		t=a[s][s];//赋初值
		for(int p=s;p<n;p++)//选取最大列主元
		{
			if(fabs(a[p][s])>t)
			{
				t=a[p][s];
				tap1=p;
				tap2=s;
			}
		}
		if(t==0)
		{
			cout<<"失败"<<endl;
			break;
		}
		if(tap1!=tap2)//换行
		{
			for(int i=s;i<n;i++)
			{
				T[i]=a[tap1][i];
				a[tap1][i]=a[s][i];
				a[s][i]=T[i];
			}
			BT=b[tap1];//给B换行
			b[tap1]=b[s];
			b[s]=BT;
		}
		for(i=s+1;i<n;i++)//消元
		{
			mm[i]=a[i][s]/a[s][s];
			for(j=s;j<n;j++)
			{
				a[i][j]=a[i][j]-mm[i]*a[s][j];
			}
			b[i]=b[i]-mm[i]*b[s];
		}
	}
	x[n-1]=b[n-1]/a[n-1][n-1];//回代步
	for(i=(n-1)-1;i>=0;i--)
	{
		for(j=i+1;j<n;j++)
		{
			temp+=a[i][j]*x[j];
		}
		x[i]=(b[i]-temp)/a[i][i];
        temp=0;
	}
	for(i=0;i<n;i++)
	{
		p[i]=x[i];
	}
}
void main()
{
	double A[m][n]={{2,4},{3,-5},{1,2},{2,1}};//系数矩阵
	double b[m]={11,3,6,7};
	double ATA[n][n];
	double ATy[n];
	double x[n];//方程组的最小二乘解
	int i,j,k;
	for(i=0;i<n;i++)
	{
		ATy[i]=0;
		for(j=0;j<n;j++)
		{
			ATA[i][j]=0;
		}
	}
	cout<<"超定方程组："<<endl;
	for(i=0;i<m;i++)
	{
		cout<<A[i][0]<<"x";
		cout.setf(ios::showpos);
		cout<<A[i][1]<<"="<<resetiosflags(ios::showpos)<<b[i]<<endl;
	}
	//ATAx=ATy
	for(i=0;i<n;i++)//计算ATA[n][n]
	{
		for(j=0;j<n;j++)
		{
			for(k=0;k<m;k++)
			{
				ATA[i][j]+=A[k][i]*A[k][j];
			}
		}
	}
	cout<<endl<<"ATA["<<n<<"]["<<n<<"]="<<endl;
	for(i=0;i<n;i++)//输出ATA[n][n]
	{
		for(j=0;j<n;j++)
		{
			cout<<setw(15)<<ATA[i][j];
		}
		cout<<endl;
	}
	for(i=0;i<n;i++)//利用A[m][n]和b[m]来计算ATy[n]
	{
			for(j=0;j<m;j++)
			{
				ATy[i]+=A[j][i]*b[j];
				
			}
	}
	cout<<endl<<"ATy["<<n<<"]="<<endl;
	for(i=0;i<n;i++)//输出ATy[n]
	{
		cout<<setw(15)<<ATy[i];
	}
	line_fun(ATA,ATy,x);//求解线性方程组
	cout<<endl<<"方程组的最小二乘解："<<endl;
	cout<<"X=(x,y)=("<<x[0]<<","<<x[1]<<")"<<endl;
	/*for(i=0;i<n;i++)
	{
		cout<<setw(10)<<"x("<<i+1<<")="<<x[i]<<endl;
	}*/

}