高斯消元

Number 数论

高斯消元

百度百科：数学上，高斯消元法（或译：高斯消去法），是线性代数规划中的一个算法，可用来为线性方程组求解。但其算法十分复杂，不常用于加减消元法，求出矩阵的秩，以及求出可逆方阵的逆矩阵。不过，如果有过百万条等式时，这个算法会十分省时。一些极大的方程组通常会用迭代法以及花式消元来解决。当用于一个矩阵时，高斯消元法会产生出一个“行梯阵式”。高斯消元法可以用在电脑中来解决数千条等式及未知数。亦有一些方法特地用来解决一些有特别排列的系数的方程组。
消元法是将方程组中的一方程的未知数用含有另一未知数的代数式表示，并将其代入到另一方程中，这就消去了一未知数，得到一解；或将方程组中的一方程倍乘某个常数加到另外一方程中去，也可达到消去一未知数的目的。消元法主要用于二元一次方程组的求解。
消元原理：
1)两方程互换，解不变；
2)一方程乘以非零数k，解不变；
3)一方程乘以数k加上另一方程，解不变
我们在线性代数中求解的时候常常把矩阵转化为三角形方阵来求解（也就是下三角形都是 0），然后我们从下往上逐步求解方程组：
①：把方程组转化成增广矩阵
②根据消元的方法从上到下依次进行消元：
ax + by = c;
dx + ey = f;
=>ax + by = c;
(d - (d/a)*a) x + (e - (d/a)*b) y = f - (d/a)*c；
=>ax + by = c;
0 x + (e - (d/a)*b) y = f - (d/a)*c；

上面介绍的其实是顺序消去法，即:
将Ax=b按照从上至下、从左至右的顺序化为上三角方程组，中间过程不对矩阵进行交换。
顺序消去法虽然编程操作简单，但是存在以下两方面限制:
(1).每次运算时，必须保证对角线上的元素不为0(即运算中的分母不为0)，否则算法无法继续进行。
(2).即使的值不为零，但如果绝对值很小，由于第k次运算中在分母位置，因此作除数会引起很大的误差，从而影响算法的稳定性。
正是由于顺序消去法会因为的值过小而引入计算误差，为了减少计算过程中舍入误差对方程组求解的影响，因此是否可以选择绝对值尽可能大的主元作为除数，基于这种思想就有了高斯消去法的改进型:列主元消去法和全主元消去法。

(1).列主元消去法基本思想:
在第k步消元前，先找出k行下所有第k列元素最大的非零元素a(r,k)，将第r行与第k行进行整行交换，这样既不影响原方程的解，也可以将绝对值最大的a(r,k）作为主元，放在除数的位置上，尽可能减小引入误差。
全主元消去法与列主元消去法类似，只不过列主元消去法是从第k列中选取一个最大的元素，与第k行进行交换。而全主元消去法是从第k行第k列开始的右下角矩阵中所有元素中选取一个最大的元素作为主元，同时交换r行与c列，从而保证稳定性。
上述引用的博客地址：

//列主元消去法
#include<iostream>
using namespace std;
double b[101][101];
int n;
double ans[101];
void swap(int x,int y){
	double t;
	for(int i=0;i<=n;i++){
		t = b[x][i];
		b[x][i] = b[y][i];
		b[y][i] = t;
	}
}
void Guass(){
	int k = 0;
	while(k<n-1){
		double Max = b[k][k];
		int maxi = k;
		for(int i=k+1;i<n;i++){
			if(b[i][k]>Max){
				maxi = i;
				Max = b[i][k];
			}
		}
		swap(k,maxi);
		for(int i=k+1;i<n;i++){
			double xi = -1 * b[i][k] / b[k][k];
			for(int j=k;j<=n;j++){
				b[i][j] += xi*b[k][j];
			}
		}
		k++;
		printf("第 %d 轮\n",k);
		for(int i=0;i<n;i++){
			for(int j=0;j<=n;j++){
				cout << b[i][j] << " ";
			}
			cout << endl;
		}	
	}
	for(int i=n-1;i>=0;i--){
		double sum = 0;
		for(int j=i+1;j<n;j++){
			sum += b[i][j]*ans[j];
		}
		ans[i] = (b[i][n]-sum)/b[i][i];
	}
	cout << "ans: ";
	for(int i=0;i<n;i++){
		cout << ans[i] << "  ";
	}
}
int main(){
	cin >> n;
	for(int i=0;i<n;i++){
		for(int j=0;j<=n;j++){
			cin >> b[i][j];
		}
	}
	Guass();
	return 0;
}