利用列主元高斯消元法求线性方程组及其算法实现

学过线性代数的朋友都知道有多种方法能用来解线性方程组,今天我给大家介绍的方法是,列主元高斯消元法及其算法实现
如何解线性方程组?
相信大家在初中就学过解方程组,如下面做这个题目

我们求解的时候,就是用消元的方法,即通过两个式子相减,可以消去一个未知数,进行两次就可以得到两个只含相同未知数的方程,这个时候再将这两个式子相减,又可以消去未知数,接下来的步骤就是求出一个未知数,然后代入先前的方程,可以求另外一个,然后再次回代,就可以求出来三个未知数了,更多未知数的方程也是可以通过类似的方法去求解的,但是一个方程组可以求解的条件是:未知数的数量不大于方程组方程的个数,因此后面的算法若是方程数小于未知数是无法求解的
为什么要用列主元消元法
在线性代数里面,想必课本上面有一个可以直接求解线性方程组的方法:克拉默法则,大家想详细了解的可以到这个网页看看link,既然克拉默法则也可以直接求解,那为什么还要用你这个消元法呢?这就牵扯到一个效率问题了,它虽然是求解线性方程组的直接方法,但却是不实用的方法!为什么这么说呢,一个n阶线性方程组,需要计算n+1个n阶行列式,即使在不计加减的情况下,其时间复杂度至少为n!(n²-1),一个算法的时间复杂度如果达到了n!级别,就不能说是一个好的算法了,何况还是n!(n²-1),而用列主元高斯消元法,n阶的线性方程组,第一次需要将第一列的n-1个元素消去,这要进行n-1次除法以及(n-1)n次乘法(因为是增广矩阵,所以共有n+1列,除去第一列,每列有n个元素需要做乘法),第二则分别是n-2和(n-2)(n-1),然后将其消元为一个上三角矩阵,回代的过程,求解第k个未知量,每求解一个未知量就需要一次除法和n-k次乘法,经过求和得到下图结果(不知道图片怎么调整,大伙将就将就)

上面说的是第一个原因,还有第二个原因,就是可以防止在消元的过程中,数据溢出!这里我们要清楚一个概念,什么是列主元?列主元便是在这一列里面,其绝对值最大的那个元素,那我们为什么要选列主元呢?就是为了防止数据溢出,若按照普通的消元方法,若你用来消元的位置上的这个系数很小,而有一个被消元的系数很大,比如一个是0.00000001,一个是99999,那么这样消元,后一个数除以前一个数的结果就溢出了
关于高斯消元法我们就说到这里,接下来开始上代码!
初始化线性方程组

double** init_Matrix(int r, int c)
{
	double** p = new double* [r];
	int d = c + 1;
	for (int i = 0; i < r; i++)
	{
		p[i] = new double[d];
		memset(p[i], 0, sizeof(double) * d);
	}
	cout << "请输入线性方程组对应的增广矩阵:" << endl;
	for (int i = 0; i < r; i++)
	{
		for (int j = 0; j < d;j++) 
		{
			cin >> p[i][j];
		}
	}
	return p;
}

这里我们将用户输入的线性方程组对应的增广矩阵进行初始化,函数返回一个二级指针(动态分配的二维数组的地址)
寻找每一列的主元素,返回其行下标

int Find_coMax(int co_index, int r, int c, double **a)
{
	if (co_index>=c||co_index<0)
	{
		fputs("查找位置不合理!\n", stderr);
		return -1;
	}
	int max = co_index;
	for (int i = co_index; i < c-1; i++)
	{
		if (fabs(a[i][co_index])<fabs(a[i+1][co_index]))
		{
			max = i + 1;
		}
	}
	return max;
}

这里找列主元,就是要找到绝对值最大的那一个,因此需要用到fabs()函数求绝对值
将列主元所在的那一行与当前行进行交换

void swap_ro(int co_swap, int co_need, int c, double **a)
{
	if (co_swap==co_need)//无须交换
	{
		return;
	}
	else
	{
		double tmp = 0;
		for (int i = co_swap; i <c+1 ; i++)
		{
			tmp = a[co_swap][i];
			a[co_swap][i] =a[co_need][i] ;
			a[co_need][i] = tmp;
		}
	}
}

注:当前行的意思是,需要用这一行的第一个不为0的元素与其他行进行消去的行
开始消元

void chang_index_co(int row_index, int r,int col, double** a)
{
	double c = 0;
	for (int i = row_index; i < r-1; i++)
	{
		c = 0;
		c = -(a[i + 1][row_index] / a[row_index][row_index]);
		for (int j = row_index; j < col+1; j++)
		{
			a[i + 1][j] += a[row_index][j] * c;
		}
	}
}

因为要使当前行的第一个不为0的元素都消去,我们乘上去的c每消一次元都要进行初始化
开始回代计算各个未知数

void calculate(int r, int c, double** a, double* p)
{
	int k = r - 1, n = c - 1;
	double tmp, q;
	for (int i = k; i >= 0; i--)
	{
		tmp = 0, q = 0;
		for (int j = n; j > i; j--)
		{
			q += a[i][j] * p[j];
		}
		p[i] = (a[i][c] - q) / a[i][i];
	}
}

解第i个未知数时,它后面的n-i个未知数都已经求出来了(因为转化成上三角矩阵,是从最后的元素往前面求解的),所以要见这n-1个解乘上其相应的系数进行求和,然后再用方程的值做差,再除以第i个未知数的系数即可
进行一次完整的列主元高斯消元解方程组(将前面的函数综合运用)

void Gauss_Maxtrix()
{
	cout << "请输入方程组的变量个数和方程个数:" << endl;
	int row = 0, col = 0;
	cin >> col >> row;
	if (col>row)
	{
		fputs("变量个数多于方程个数,无法求解!\n", stderr);
		return;
	}
	double** a = init_Matrix(row, col);
	int lin = 0;
	int Trow = row > col ? col : row;
	for (int i = 0; i <Trow ; i++)
	{
		lin = Find_coMax(i, row, col, a);
		if (lin==-1)
		{
			fputs("参数不合理!\n", stderr);
			return;
		}
		swap_ro(i, lin, col, a);
		chang_index_co(i, row, col,a);
	}
	double* p = new double[col];
	memset(p, 0, sizeof(double) * col);
	calculate(row,col,a,p);
	cout << "方程组的近似解为:" << endl;
	for (int i = 0; i < col; i++)
	{
		cout << p[i] << endl;
	}
	for (int i = 0; i < row; i++)
	{
		delete []a[i];
	}
	delete []a;
	delete []p;
}

大家在求解完以后,还是要养成释放动态内存的好习惯
完整代码

#include<iostream>
#include<cmath>
using namespace std;
/*
测试数据
2 2
2 3  5
1 -1 0

3 3
3 2 -3 -2
1 1  1  6
1 2 -1  2
*/
double** init_Matrix(int r, int c)
{
	double** p = new double* [r];
	int d = c + 1;
	for (int i = 0; i < r; i++)
	{
		p[i] = new double[d];
		memset(p[i], 0, sizeof(double) * d);
	}
	cout << "请输入线性方程组对应的增广矩阵:" << endl;
	for (int i = 0; i < r; i++)
	{
		for (int j = 0; j < d;j++) 
		{
			cin >> p[i][j];
		}
	}
	return p;
}
//当前列
int Find_coMax(int co_index, int r, int c, double **a)
{
	if (co_index>=c||co_index<0)
	{
		fputs("查找位置不合理!\n", stderr);
		return -1;
	}
	int max = co_index;
	for (int i = co_index; i < c-1; i++)
	{
		if (fabs(a[i][co_index])<fabs(a[i+1][co_index]))
		{
			max = i + 1;
		}
	}
	return max;
}
//co_need 值最大主元所在的行下标 co_swap最上面的行下标
void swap_ro(int co_swap, int co_need, int c, double **a)
{
	if (co_swap==co_need)//无须交换
	{
		return;
	}
	else
	{
		double tmp = 0;
		for (int i = co_swap; i <c+1 ; i++)
		{
			tmp = a[co_swap][i];
			a[co_swap][i] =a[co_need][i] ;
			a[co_need][i] = tmp;
		}
	}
}
void chang_index_co(int row_index, int r,int col, double** a)
{
	double c = 0;
	for (int i = row_index; i < r-1; i++)
	{
		c = 0;
		c = -(a[i + 1][row_index] / a[row_index][row_index]);
		for (int j = row_index; j < col+1; j++)
		{
			a[i + 1][j] += a[row_index][j] * c;
		}
	}
}
void calculate(int r, int c, double** a, double* p)
{
	int k = r - 1, n = c - 1;
	double tmp, q;
	for (int i = k; i >= 0; i--)
	{
		tmp = 0, q = 0;
		for (int j = n; j > i; j--)
		{
			q += a[i][j] * p[j];
		}
		p[i] = (a[i][c] - q) / a[i][i];
	}
}
void Gauss_Maxtrix()
{
	cout << "请输入方程组的变量个数和方程个数:" << endl;
	int row = 0, col = 0;
	cin >> col >> row;
	if (col>row)
	{
		fputs("变量个数多于方程个数,无法求解!\n", stderr);
		return;
	}
	double** a = init_Matrix(row, col);
	int lin = 0;
	int Trow = row > col ? col : row;
	for (int i = 0; i <Trow ; i++)
	{
		lin = Find_coMax(i, row, col, a);
		if (lin==-1)
		{
			fputs("参数不合理!\n", stderr);
			return;
		}
		swap_ro(i, lin, col, a);
		chang_index_co(i, row, col,a);
	}
	double* p = new double[col];
	memset(p, 0, sizeof(double) * col);
	calculate(row,col,a,p);
	cout << "方程组的近似解为:" << endl;
	for (int i = 0; i < col; i++)
	{
		printf("%f\t",p[i]);
	}
	cout<<endl;
	for (int i = 0; i < row; i++)
	{
		delete []a[i];
	}
	delete []a;
	delete []p;
}
int main(void) {
	Gauss_Maxtrix();
	system("pause");
	return 0;
}

大家还有什么更好的算法欢迎跟我一起探讨~