線性方程直接求法——高斯消元法(一）

宣告：僅為個人筆記，相關內容僅供參考，

對于低階的線性方程組，可用克萊姆法則計算其解，對高階線性方程組，計算量巨大，應采用數值計算方法中的直接法和迭代法，這里介紹直接法，
直接法——不考慮計算程序中的舍入誤差，經有限次數的運算便可求得方程組的準確解的方法，克萊姆法則求解即為直接法，克萊姆法則漸近復雜度為O(n*n!)雖然很穩定但是計算量太大了！我們這里引入其他的方法，過多的文字贅述也大可不必，這里簡單講解一下實作思路就行，

解決方法思想：

一、Gausss順序消元法：

1、化為三角矩陣：將一個非奇異矩陣經過初等行變換得到上或下三角矩陣
2、回代方程求解：上三角從Xn開始回代方程求解，上三角則是X1開始

二、Gauss-Jordan消元法：

相對高斯消元法，它最后得到的線性方程組更易求解，它得到的是一個簡化行列式，
1、化為三角矩陣：將一個非奇異矩陣經過初等行變換得到上或下三角矩陣
2、化為對角型：除對角元素其他都消去為0

示例：求解下面這個方程：


代碼實作：

//注：下面的代碼中并沒有把A矩陣和b矩陣放在一起合成增廣矩陣而是單獨存盤的，其實都一樣
#include<iomanip>
#include<iostream>
using namespace std;
//一個增廣矩陣輸出的函式
void print(double A[], double b[], int n);
//Gauss順序消元法
void GaussSequence_x(double a[], double b[], int n)
{
	double temp;
	//下面用了一個單獨的X陣列來存盤X的值，目的是為了方便理解
	//但其實也沒有這個必要，下面的 GaussSequence_b函式可理解、參考
	double* x = new double[n];
	//第一步：通過初等行變換將矩陣轉換為上三角矩陣：
	for (int i = 0; i < n - 1; i++) //N-1次消元
	{
		//第i列將首元下的n-i個元素消去
		//即在二維陣列中表示為i+1到n-1行i列
		for (int j = i + 1; j < n; j++)
		{
			//先取與首元倍數
			temp = a[j * n + i] / a[i * n + i];
			//第j這一行元素做初等行變換
			for (int k = i; k < n; k++)
			{
				a[j * n + k] -= temp * a[i * n + k];
			}
			//增廣矩陣做初等行變換，b矩陣也要做
			b[j] -= temp * b[i];
		}
	}
	cout << "上三角矩陣為：" << endl; print(a, b, n);
	//回代方程計算
	//X從N-1開始計算，第一次回圈便可以可到下面這行代碼的效果
	//x[n - 1] = b[n - 1] / a[(n - 1) * n + (n - 1)];
	//將外層回圈i的起始條件改為n-2上面的這行代碼可以先寫在外面，更加直觀，但沒必要
	for (int i = n - 1; i >= 0; i--)
	{
		x[i] = b[i];
		//N-1行后面操作的每行算的X[i]都要用b[i]減去
		//前面操作算出的X和對應A元素乘積和
		for (int j = n - 1; j > i; j--)
		{
			x[i] -= x[j] * a[i * n + j];
		}
		//最后減去乘積之和之后再除以對角線元素就是對應的X值了
		x[i] /= a[i * n + i];
	}
	for (int i = 0; i < n; i++)
	{
		cout << "X[" << i + 1 << "]= " << x[i] << endl;
	}
	for (int i = 0; i < 20 * n; i++) cout << "_"; cout << endl;
	delete[]x;
}
//Gauss順序消元法
//和上面是一樣的，唯有不同就是沒有單獨創建一個X陣列，這里自行理解就行
void GaussSequence_b(double a[], double b[], int n)
{
	double temp;
	//第一步：通過初等行變換將矩陣轉換為上三角矩陣：
	for (int i = 0; i < n - 1; i++) //N-1次消元
	{
		for (int j = i + 1; j < n; j++)
		{
			temp = a[j * n + i] / a[i * n + i];
			for (int k = i; k < n; k++)
			{
				a[j * n + k] -= temp * a[i * n + k];
			}
			b[j] -= temp * b[i];
		}
	}
	cout << "上三角矩陣為：" << endl; print(a, b, n);
	for (int i = n - 1; i >= 0; i--)
	{
		for (int j = n - 1; j > i; j--)
		{
			b[i] -= b[j] * a[i * n + j];
		}
		b[i] /= a[i * n + i];
	}
	for (int i = 0; i < n; i++)
	{
		cout << "X[" << i + 1 << "]= " << b[i] << endl;
	}
	for (int i = 0; i < 20 * n; i++) cout << "_"; cout << endl;
}
//Gauss-Jordan消元法
//這里也不單獨創建X陣列來存盤未知數，
//而是和上面函式一樣用b陣列來放解得的X值
void GaussJordan_b(double a[], double b[], int n)
{
	double temp;
	//第一步：通過初等行變換將矩陣轉換為上三角矩陣：
	for (int i = 0; i < n - 1; i++) //N-1次消元
	{
		for (int j = i + 1; j < n; j++)
		{
			temp = a[j * n + i] / a[i * n + i];
			for (int k = i; k < n; k++)
			{
				a[j * n + k] -= temp * a[i * n + k];
			}
			b[j] -= temp * b[i];
		}
	}
	cout << "上三角矩陣為：" << endl; print(a, b, n);
	//化為對角型
	for (int i = n - 1; i >= 0; i--)
	{
		b[i] /= a[i * n + i];
		for (int j = 0; j < i; j++)
		{
			b[j] -= b[i] * a[j * n + i];
		}
	}
	print(a, b, n);
	for (int i = 0; i < n; i++)
	{
		cout << "X[" << i + 1 << "]= " << b[i] << endl;
	}
	for (int i = 0; i < 20*n; i++) cout << "_"; cout << endl;
}
int main()
{
	int n = 3;
	double* b = new double[n] { 3, 1, -7};
	double* A = new double[(double)n * n]{ 2,2,3,4,7,7,-2,4,5};
	GaussSequence_x(A, b, n);
	//GaussSequence_b(A, b, n);
	//GaussJordan_b(A, b, n);
	delete[]A;
	delete[]b;
	system("pause");
	return 0;
}
void print(double A[], double b[], int n)
{
	for (int i = 0; i < n; i++)
	{
		for (int j = 0; j < n; j++)
		{
			cout <<left<<setw(15)<< A[i * n + j];
		}
		cout << left << setw(15) << b[i] << endl;
	}
	for (int i = 0; i < 20 * n; i++) cout << "_"; cout << endl;
}

單次運行結果：

注意事項：

(1)每次運算時，必須保證對角線上的元素不為0(即運算中的分母不為0)，否則演算法無法繼續進行，
(2)首元如果絕對值很小，由于第k次運算中在分母位置，因此作除數帶來的舍入誤差有時超過預期，從而影響演算法的穩定性，d(X/Y)=(YdX-XdY)/Y^2
(3)漸進復雜度為O(n^3),相比克萊姆法則計算量已經減少很大，但是計算階數比較高的方程時計算量仍很大，
以上問題都有解決方法，后面都會講到，