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

Posted 2021-10-25 喜欢吃糖糖

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

声明：仅为个人笔记，相关内容仅供参考。

---

目录

• 线性方程直接求法——高斯消元法(一）
• • • 解决方法思想：
    • • 一、Gausss顺序消元法：
      • 二、Gauss-Jordan消元法：
    • 代码实现：
    • 注意事项：

---

对于低阶的线性方程组，可用克莱姆法则计算其解，对高阶线性方程组，计算量巨大，应采用数值计算方法中的直接法和迭代法。这里介绍直接法。
直接法——不考虑计算过程中的舍入误差，经有限次数的运算便可求得方程组的准确解的方法，克莱姆法则求解即为直接法。克莱姆法则渐近复杂度为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),相比克莱姆法则计算量已经减少很大，但是计算阶数比较高的方程时计算量仍很大。
以上问题都有解决方法，后面都会讲到。