线性方程组的直接法

求解线性方程组的最基本而古老的直接法即为消去法，利用一系列递推公式计算有限步即可得到方程组的精确解，当然有时也会有舍入误差等，也可能引致严重精度问题。系数矩阵设置为float a[ ]，方程组的结果为float y[ ]，未知数矩阵为float x[ ],自然y[ ]、x[ ]都只有一列，下面介绍两种经典的解法。

1>Jordan消去法

将每列的系数消成只剩一个，最后一一对应，即求解完毕。

for(int i=0;i<row;i++)
{
	//选主元素的过程
	largest=0;
	for(int r=i;r<row;r++)
		if(a[r][i]>largest)
		{
			largest=a[r][i];
			locate=r;
		}
	if(locate!=i)
	{
		ExchangeRow(a,locate,i,i);
		ExchangeRow(b,locate,i,1);
	}
	//将当前运算行系数变化
	b[i]=b[i]/a[i][i];
	for(int j=line-1;j>=i;j++)a[i][j]=a[i][j]/a[i][i];		
	//将其它行进行消去
	for(int k=0;k<row;k++)
	{
		if(k==i)continue;
		for(int j=i+1;j<line;j++)//其它行只需从第i+1列开始即可
			a[k][j]=a[k][j]-a[i][j]*a[k][i];
		b[k]=b[k]-b[i]*a[k][i];
		a[k][i]=0;
	}
}

 

2>Gauss消去法

将系数矩阵消成三角形（这里以上三角为例），再将其进行回代处理。

for(i=0;i<row;i++)
{
	//选主元素的过程
	largest=0;
	for(int r=i;r<row;r++)
		if(a[r][i]>largest)
		{
			largest=a[r][i];
			locate=r;
		}
	if(locate!=i)
	{
		ExchangeRow(a,locate,i,i);
		ExchangeRow(b,locate,i,1);
	}
	//将当前计算行系数变化
	b[i]=b[i]/a[i][i];
	for(int j=line-1;j>=i;j++)a[i][j]=a[i][j]/a[i][i];		
	//行与列只需从第i+1列开始即可
	for(int k=i+1;k<row;k++)
	{
		for(int j=i+1;j<line;j++)
			a[k][j]=a[k][j]-a[i][j]*a[k][i];
		b[k]=b[k]-b[i]*a[k][i];
		cout<<b[k]<<endl;
		a[k][i]=0;
	}
}
//上三角的回代求解
for(i=row-1;i>=0;i--)
{
	sum=0;
	for(int j=i+1;j<line;j++)sum=sum+a[i][j]*x[j];
	x[i]=b[i]-sum;
}

上述的选主元素即为防止出现除数为0或舍入误差等情况的出现而设定的步骤，下面的图将显示系数矩阵在两种消元方法中的变化情况：

由已知可得，Jordan消去法的乘除法计算工作量为∑(n-k+1)*n=(n^3)/2，而Gauss消去法的乘除法计算工作量为(n^3)/3，由此可见后者的方便与快捷~

 

3>三对角矩阵的追赶法

对于某些特殊的矩阵，发现其只在主对角线与其相邻的两条次对角线上，如果用上述两种消去法不仅效率低下而且浪费资源，所以我们在上述的消去法前提下进行专门的改进。由于Gauss消去的快捷，所以我们也选择消成一个上三角，再进行回代处理。首先源于对角占优的改进，以前我们是选择单列最大值进行选主元素，这里我们改为主对角线元素为最大值。

u[0]=m[0][1]/m[0][0];
y[0]=y[0]/m[0][0];
for(i=1;i<row;i++)
{
	u[i]=m[i][i+1]/(m[i][i]-u[i-1]*m[i][i-1]);//u当中的i只有row-1个
	y[i]=(y[i]-y[i-1]*m[i][i-1])/(m[i][i]-u[i-1]*m[i][i-1]);
	cout<<u[i]<<'\t'<<y[i]<<'\n';
}
//回代过程
x[row-1]=y[row-1];
for(i=row-2;i>=0;i--)
	x[i]=y[i]-u[i]*x[i+1];

由上易知，追赶法的乘除计算工作量约为5n。

 

4>Cholesky、Doolittle、Crout三角分解
通过Gauss消去法，我们看到了三角矩阵的求解方便快捷，所以如果能把系数矩阵转化成三角形式的话，我们的求解速度会加快许多。Cholesky选择将一个系数矩阵转化成为L*D*~L（~L表示L的转置矩阵，D则为对角阵）的形式，而Doolittle与Crout则选择将系数矩阵转化成为L*U的形式（Crout转化中的U为上三角，Doolittle转化中的L为下三角），为了使分解具有唯一性，其对角主元素必须单位化。

建立矩阵类MATRIX，方程组为a*x=b，a、b、x都为类的对象。

//Cholesky矩阵分解
for(int i=0;i<a.row;i++)
{
	L.matrix[i][i]=1;//单位化
	for(int j=0;j<i;j++)
	{
		sum=0;
		for(int k=0;k<j-1;k++)
			sum+=D.matrix[k][k]*L.matrix[i][k]*L.matrix[j][k]/D.matrix[j][j];
		L.matrix[i][j]=(a.matrix[i][j]-sum)/D.matrix[j][j];
	}
	sum=0;
	for(j=0;j<i-1;j++)
		sum+=D.matrix[j][j]*L.matrix[i][j]*L.matrix[i][j];
	D.matrix[i][i]=a.matrix[i][i]-sum;
}

Cholesky分解公式：Lij=(aij-∑Dk*Lik*Ljk)/Dj，Di=aii-∑Dk*Lik*Lik。

//Doolittle矩阵分解
for(int i=0;i<a.row;i++)
{
	L.matrix[i][i]=1;//单位化
	for(int j=0;j<i-1;j++)
	{
		sum=0;
		for(int k=0;k<j-1;k++)
			sum+=L.matrix[i][k]*U.matrix[k][j];
		L.matrix[i][j]=(a.matrix[i][j]-sum)/U.matrix[j][j];
	}
	for(j=i;j<a.line;j++)//注意j从i开始
	{
		sum=0;
		for(int k=0;k<i-1;k++)
			sum+=L.matrix[i][k]*U.matrix[k][j];
		U.matrix[i][j]=a.matrix[i][j]-sum;
	}
}

Doolittle分解公式：Lij=(aij-∑Lik*Ukj)/Ujj，Uij=aij-∑Lik*Ukj。

//Crout矩阵分解
for(int i=0;i<a.row;i++)
{
	U.matrix[i][i]=1;//单位化
	for(int j=0;j<i-1;j++)
	{
		sum=0;
		for(int k=0;k<j-1;k++)
			sum+=L.matrix[i][k]*U.matrix[k][j];
		L.matrix[i][j]=a.matrix[i][j]-sum;
	}
	for(j=i;j<a.line;j++)//注意j从i开始
	{
		sum=0;
		for(int k=0;k<i-1;k++)
			sum+=L.matrix[i][k]*U.matrix[k][j];
		U.matrix[i][j]=(a.matrix[i][j]-sum)/L.matrix[i][i];
	}
}

Crout分解公式：Lij=aij-∑Lik*Ukj，Uij=(aij-∑Lik*Ukj)/Lii。