线性代数代码实现(七)求解线性方程组(C++)

线性代数代码实现(七)求解线性方程组(C++),第1张

线性代数代码实现(七)求解线性方程组(C++)

前言:

        上次博客,我写了一篇关于定义矩阵除法并且代码的文章。矩阵除法或许用处不大,不过在那一篇文章中,我认为比较好的一点是告诉了大家一种计算方法,即:若矩阵  已知且可逆,矩阵  已知,并且  ,求解矩阵 B 。我认为这种初等行变换的方法还是挺好的。

        在本篇文章中,我和大家探讨一下线性代数里面一个重要的知识——线性方程组及其解法。

一、线性代数知识回顾:

我们先探讨一下二元一次方程组的解法:

        相信这个解法大家已经很熟悉了,将第一个式子的 -2 倍加到第二个式子上,就可以消掉第二个式子中的 x 了,然后第二个式子就只剩下 y 一个未知数了,就可以轻松解出 y ,然后将 y 代入第一个式子,就可以轻松解出 x ,到此,二元一次方程组求解完成。

        从几何意义上看,这两个方程对应着二维平面上的两条直线,并且这两条直线交于一点(1,1),因此有唯一解。

        有人可能会想到,不是所有的二元一次方程组都有解,因此我们看一看下面的例子: 

         大家看到,这个这个将第二个式子中两个未知数都化没了,并且观察这个方程组可知,它有无数解,其实容易看到,第二个式子就是第一个式子的 2 倍,实际上,这两个式子是等价的,第二个式子所描述的信息完全可以用第一个式子描述,也就是说这个方程组中 “有效方程” 的个数为 1 ,而只靠一个式子无法固定两个未知数,因此方程组有无数解。

        从几何意义上看,这两个式子就是二维平面中的两个直线,而这两个直线是重合的,也就是说,两个直线相交的点有无数个,即:方程组有无数解。

        那么,二元一次方程组有没有可能无解呢?看看下面的例子:

         可以看到,消元后,第二个方程变为了一个不可能成立的式子。容易看出,这两个方程是矛盾的,因此是无解的。

        从几何意义上看,这两个二维平面上的直线平行且不重合,没有交点,因此无解。

        从上面二元一次方程组的回顾中,我们可以将方程组推广到 n 维情形:

         同样按照上面的消元方法,不过这次稍微复杂,首先消去第 2 个到第 n 个方程的未知数  然后消去第 3 个到第 n 个方程的未知数  ,以此类推,直到消去最后一个方程的未知数  ,然后自下而上,先求出  ,代入倒数第二个方程,求出  ,代入倒数第三个方程,依次类推,便可求出方程组的解。大家仔细看看消元的过程,是不是和初等行变换十分像,我们其实可以将上述方程组这样表示:

设                                  

上述方程组可以表示为:

为了表示方便,我们引入增广矩阵:

         当我们对未知数规定书写顺序,那么增广矩阵就和线性方程组就 一 一 对应了,对方程组进行消元就等价于对增广矩阵  进行初等行变换,方程组的有效方程的个数就等于增广矩阵的秩,也等于系数矩阵的秩,也就是说,如果系数矩阵满秩,那么方程组有唯一解。当系数矩阵不满秩时,若没有矛盾方程出现,则方程组有无数解,若有矛盾方程出现,则方程组无解。

        从几何意义来看,这 n 个方程组相当于 n 维空间中 n 个 n-1 维的平面。n 个平面交于一点等价于方程组中的 “有效方程” 的个数为 n ,等价于系数矩阵  满秩,等价于增广矩阵  的秩为 n ,等价于方程组有唯一解。

二、算法描述:

        输入系数矩阵  以及 向量  ,求解向量  ,使得  。

        1. 求出增广矩阵

        2. 将增广矩阵初等行变换,化为上三角形,在这个过程中,判断增广矩阵的秩是否是 n ,如果不是 n ,则说明方程组无解或有无数解。

        3. 自下而上依次计算出  。

         详情请看下面代码:

三、代码实现:

        类定义为:

class Mat
{
public:
	int m = 1, n = 1; //行数和列数
	double mat[N][N] = { 0 };  //矩阵开始的元素

	Mat() {}
	Mat(int mm, int nn)
	{
		m = mm; n = nn;
	}
	void create();//创建矩阵
	void Print();//打印矩阵

	void augmat(Mat a, Mat b);//求矩阵 a 和向量 b 的增广矩阵
	bool solve(Mat a, Mat b); //求线性方程组的解
};

        其中solve函数求解线性方程组,其定义如下:

bool Mat::solve(Mat a, Mat b) //a 为方阵 ,b 为列向量
//求线性方程组的解(ax=b ,求 x),矩阵 a 为方阵并且方程组有唯一解时返回 true
{
	if (a.n != a.m)//只求解是方阵时的情形
	{
		cout << "系数矩阵不是方阵" << endl;
		return false; //返回false
	}
	m = a.n; n = 1; //解向量中必定有 a.n( a.m )个分量,是 a.n * 1 的列向量
	Mat aa;
	aa.augmat(a, b); //求增广矩阵

	//下面代码将增广矩阵化为上三角矩阵,并判断增广矩阵秩是否为 n
	for (int i = 1; i <= aa.m; i++)
	{
		//寻找第 i 列不为零的元素
		int k;
		for (k = i; k <= aa.m; k++)
		{
			if (fabs(aa.mat[k][i]) > 1e-10) //满足这个条件时,认为这个元素不为0
				break;
		}
		if (k <= aa.m)//说明第 i 列有不为0的元素
		{
			//交换第 i 行和第 k 行所有元素
			for (int j = i; j <= aa.n; j++)//从第 i 个元素交换即可,因为前面的元素都为0
			{//使用aa.mat[0][j]作为中间变量交换元素
				aa.mat[0][j] = aa.mat[i][j]; aa.mat[i][j] = aa.mat[k][j]; aa.mat[k][j] = aa.mat[0][j];
			}
			double c;//倍数
			for (int j = i + 1; j <= aa.m; j++)
			{
				c = -aa.mat[j][i] / aa.mat[i][i];
				for (k = i; k <= aa.n; k++)
				{
					aa.mat[j][k] += c * aa.mat[i][k];//第 i 行 a 倍加到第 j 行
				}
			}
		}
		else //没有找到则说明系数矩阵秩不为 n ,说明方程组中有效方程的个数小于 n
		{
			cout << "系数矩阵奇异,线性方程组无解或有无数解" << endl;
			return false;
		}
	}
	//自下而上求解
	for (int i = a.m; i >= 1; i--)
	{
		mat[i][1] = aa.mat[i][aa.n];
		for (int j = a.m; j > i; j--)
		{
			mat[i][1] -= mat[j][1] * aa.mat[i][j];
		}
		mat[i][1] /= aa.mat[i][i];
	}
	return true;
}

线性方程组的求解就完成了,这种求解方法叫 “高斯消元法” 。下面附上完整代码:

#include
#include 
#define N 10
using namespace std;
class Mat
{
public:
	int m = 1, n = 1; //行数和列数
	double mat[N][N] = { 0 };  //矩阵开始的元素

	Mat() {}
	Mat(int mm, int nn)
	{
		m = mm; n = nn;
	}
	void create();//创建矩阵
	void Print();//打印矩阵

	void augmat(Mat a, Mat b);//求矩阵 a 和向量 b 的增广矩阵
	bool solve(Mat a, Mat b); //求线性方程组的解
};

void Mat::create()
{
	for (int i = 1; i <= m; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			cin >> mat[i][j];
		}
	}
}
void Mat::Print()
{
	for (int i = 1; i <= m; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			cout << mat[i][j] << "t";
		}
		cout << endl;
	}
}
void Mat::augmat(Mat a, Mat b)//求矩阵 a 和向量 b 的增广矩阵
{
	m = a.m; n = a.n + 1;
	for (int i = 1; i <= a.m; i++)
	{
		for (int j = 1; j <= a.n; j++)
		{
			mat[i][j] = a.mat[i][j];
		}
		mat[i][n] = b.mat[i][1];
	}
}

bool Mat::solve(Mat a, Mat b) //a 为方阵 ,b 为列向量
//求线性方程组的解(ax=b ,求 x),矩阵 a 为方阵并且方程组有唯一解时返回 true
{
	if (a.n != a.m)//只求解是方阵时的情形
	{
		cout << "系数矩阵不是方阵" << endl;
		return false; //返回false
	}
	m = a.n; n = 1; //解向量中必定有 a.n( a.m )个分量,是 a.n * 1 的列向量
	Mat aa;
	aa.augmat(a, b); //求增广矩阵

	//下面代码将增广矩阵化为上三角矩阵,并判断增广矩阵秩是否为 n
	for (int i = 1; i <= aa.m; i++)
	{
		//寻找第 i 列不为零的元素
		int k;
		for (k = i; k <= aa.m; k++)
		{
			if (fabs(aa.mat[k][i]) > 1e-10) //满足这个条件时,认为这个元素不为0
				break;
		}
		if (k <= aa.m)//说明第 i 列有不为0的元素
		{
			//交换第 i 行和第 k 行所有元素
			for (int j = i; j <= aa.n; j++)//从第 i 个元素交换即可,因为前面的元素都为0
			{//使用aa.mat[0][j]作为中间变量交换元素
				aa.mat[0][j] = aa.mat[i][j]; aa.mat[i][j] = aa.mat[k][j]; aa.mat[k][j] = aa.mat[0][j];
			}
			double c;//倍数
			for (int j = i + 1; j <= aa.m; j++)
			{
				c = -aa.mat[j][i] / aa.mat[i][i];
				for (k = i; k <= aa.n; k++)
				{
					aa.mat[j][k] += c * aa.mat[i][k];//第 i 行 a 倍加到第 j 行
				}
			}
		}
		else //没有找到则说明系数矩阵秩不为 n ,说明方程组中有效方程的个数小于 n
		{
			cout << "系数矩阵奇异,线性方程组无解或有无数解" << endl;
			return false;
		}
	}
	//自下而上求解
	for (int i = a.m; i >= 1; i--)
	{
		mat[i][1] = aa.mat[i][aa.n];
		for (int j = a.m; j > i; j--)
		{
			mat[i][1] -= mat[j][1] * aa.mat[i][j];
		}
		mat[i][1] /= aa.mat[i][i];
	}
	return true;
}


int main()
{
	Mat a(3, 3), b(3, 1);
	cout << "请输入 " << a.m << "*" << a.n << " 的矩阵:" << endl;
	a.create();
	cout << "请输入 " << b.m << "*" << b.n << " 的列向量:" << endl;
	b.create();

	Mat x;
	if (x.solve(a, b))//求线性方程组的解向量
	{
		cout << "解向量如下:" << endl << endl;
		x.Print();
	}
	return 0;
}

若有不足之处,欢迎大家指正!

欢迎分享,转载请注明来源:内存溢出

原文地址: http://outofmemory.cn/zaji/5432463.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2022-12-11
下一篇 2022-12-11

发表评论

登录后才能评论

评论列表(0条)

保存