高斯列主元消去法解線性方程組

最近在看慣導的東西，然後想要用C++解慣導控制方程，然後就重頭把C++解方程組這方面的知識回顧了一下，首先就是高斯列主元消去法，這個方法還算實用，這裡以3*3的矩陣為例，裡面注釋很詳細，各位小白可以參考參考hhh

//順序高斯消去法求線性方程組的解  AX=B,A為系數矩陣，B為常向量矩陣，X為矩陣的解
//為了減小誤差，采用列主元消去法，從第k列的a[k][k]及其以下的各元素中選取絕對值最大的元素，然後行變換将這兩行交換位置
//對于非奇異矩陣并且消元過程中第K步的主對角線上的元素不為0,就可以用高斯消元法來求解

#include<iostream>
const int n=3;
void Gauss_elimination(double a[n][n],double b[n])
{
	double x[n];//解的存儲數組
	int i,j,k;
	double aa[n];//左上的主對角元素為0時系數矩陣的行互相交換
	double bb;//左上的主對角元素為0時系數矩陣的行互相交換
	double mm;//找出最大的a[m][k]行的a[m][k]
	int m;//找出最大的a[m][i]行的行m值
	double c[n];//存儲初等行變換的系數，用于行的相減
	for(k=0;k<n-1;k++)	 //消元的整個過程如下，總共n-1次消元過程。
	{
		mm=a[k][k];m=k;bb=0;
		for(int i=1;i<n-k;i++)//找出最大的絕對值a[m][k]行的行m值
		{
			if(mm<abs(a[k+i][k]))
			{m=k+i;mm=abs(a[k+i][k]);}
		}
		if(mm==0)              //此種情況表明a[k][k]以及下面全部是0，則方程有多解或者無解，然後結束函數
		{
			std::cout<<"此方程無解或多解"<<std::endl;
			return;
		}
		bb=b[m];
		b[m]=b[k];
		b[k]=bb;
		for(int i=k;i<n;i++)//将最大的a[m][k]行移到最上面
		{
			aa[i]=a[m][i];
			a[m][i]=a[k][i];
			a[k][i]=aa[i];
		}
		for(i=k+1;i<n;i++)//求出第K次初等行變換的系數
			c[i]=a[i][k]/a[k][k];
		for(i=k+1;i<n;i++)      //第K次的消元計算
			{
				for(j=0;j<n;j++)
				{
					a[i][j]=a[i][j]-c[i]*a[k][j];
				}
				b[i]=b[i]-c[i]*b[k];
			}
	}
	std::cout<<"輸出變換後的系數矩陣："<<std::endl;//輸出變換後的系數矩陣及常向量矩陣
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			 std::cout<<a[i][j]<<"   ";
		}                                      
		std::cout<<"b["<<i+1<<"]="<<b[i]<<std::endl;
	}
	//先計算出最後一個未知數；
    x[n-1]=b[n-1]/a[n-1][n-1];
    //求出每個未知數的值
    for(i=n-2;i>=0;i--)
    {
        double sum=0;
        for(j=i+1;j<n;j++)
        {
            sum+=a[i][j]*x[j];
        }
        x[i]=(b[i]-sum)/a[i][i];
    }
    std::cout<<"此方程組的解為："<<std::endl;
    for(int i=0;i<n;i++)
	{
		std::cout<<"x["<<i+1<<"]="<<x[i]<<std::endl;
	}
}
           

主函數如下：

int main()
{
	double a[3][3]={1,1,1,0,4,-1,2,-2,1};
	double b[3]={6,5,1};
	Gauss_elimination(a,b);
	system("pause");
	return 0;
}
           

運作結果如下：