【数学】高斯-约旦消元法

给定 nn 元一次方程组

{a1,1x1+a1,2x2++a1,nxn=b1a2,1x1+a2,2x2++a2,nxn=b2an,1x1+an,2x2++an,nxn=bn⎪ ⎪ ⎪⎪ ⎪ ⎪a1,1x1+a1,2x2++a1,nxn=b1a2,1x1+a2,2x2++a2,nxn=b2an,1x1+an,2x2++an,nxn=bn

请求出方程组的解的情况:

  • 无解;

  • 无穷多解;

  • 唯一解。

对于这样的问题,我们可以使用 高斯消元法 进行求解,当然高斯消元法有一个回代的过程,代码略长,而且精度较低。

所以我们隆重推出 高斯-约旦消元法 !!!

回顾一下我们是怎么手算的,一般用的都是 加减消元法,普通高斯和高斯-约旦用的都是加减消元。

在此之前,我们需要了解一下矩阵初等变换。

在线性代数中,矩阵初等行变换 是指以下三种变换类型 :

  1. 交换矩阵的两行;

  2. 用一个非零数 kk 乘矩阵的某一行所有元素;

  3. 把矩阵的某一行所有元素乘以一个数 kk 后加到另一行对应的同一列的元素上;

类似地,把以上的 改为 便得到 矩阵初等列变换 的定义。

矩阵初等行变换与初等列变换合称为 矩阵初等变换

若矩阵 AA 经过有限次的初等行变换变为矩阵 BB,则矩阵 AA 与矩阵 BB 行等价;若矩阵 AA 经过有限次的初等列变换变为矩阵 BB,则矩阵 AA 与矩阵 BB 列等价;若矩阵 AA 经过有限次的初等变换变为矩阵 BB,则矩阵 AA 与矩阵 BB 等价

当然列的用不着

首先有一个由系数构成的 n×nn×n 的矩阵

[a1,1a1,2a1,na2,1a2,2a2,nan,1an,2an,n]⎢ ⎢ ⎢ ⎢ ⎢a1,1a1,2a1,na2,1a2,2a2,nan,1an,2an,n⎥ ⎥ ⎥ ⎥ ⎥

然后是一个由常数构成的 n×1n×1 的列向量

[b1b2bn]⎢ ⎢ ⎢ ⎢b1b2bn⎥ ⎥ ⎥ ⎥

把它们放在一起构成一个 n×(n+1)n×(n+1) 的增广矩阵:

[a1,1a1,2a1,nb1a2,1a2,2a2,nb2an,1an,2an,nbn]⎢ ⎢ ⎢ ⎢ ⎢a1,1a1,2a1,nb1a2,1a2,2a2,nb2an,1an,2an,nbn⎥ ⎥ ⎥ ⎥ ⎥

我们遍历每一列,对于第 ii 列选出最大的 未处理过的行 上的数作为主元,意味着加减消元后除了主元这一行外其他行的第 ii 列都为 00

选最大的作为主元的原因是避免精度误差。

然后就是加减消元了。

举个例子

{2xy+z=14x+yz=5x+y+z=02xy+z=14x+yz=5x+y+z=0

增广矩阵就是

[211141151110]211141151110

先是第 11 列,选 44 为主元。

44 在第 22 行,将第 22 行与第 11 行交换。

[411521111110]411521111110

对于第 22 行,第一列上的 22441212

[411524×1211×121(1)×1215×121110]⎢ ⎢ ⎢411524×1211×121(1)×1215×121110⎥ ⎥ ⎥

化简得

[411503232321110]

对第 3 行的处理同理

[411503232320345454]

现在到了第 2 列,未处理过的是 2,3 行,选最大的 34 为主元。

将第 3 行与第 2 行交换

[411503454540323232]

消元得

[408320303454540044]

3 列,未处理过的是第 3 行,选 4 作主元。

[4004034000044]

这样就把原来的矩阵通过初等变换,使得系数矩阵只有主对角线位置的元素非零,即一个对角矩阵。

上面那个矩阵的意思是

{4x=434y=04z=4

所以再把系数除过去就得到

{x=1y=0z=1

请自行检验。

时间复杂度为 O(n3)

当然方程还可能是无解或无穷多解。

考虑一个一元一次方程 ax=b 的解的情况:

  1. 无解:a=0,b0
  2. 无穷多解:a=b=0
  3. 唯一解:a0
  4. 所以当发现某一列的主元为 0 时,因为主元是最大的,所以意味着这一列全都是 0,那么要么无解,要么有无穷多解。

P3389 【模板】高斯消元法

Code

#include <iostream>
#include <cstdio>
#include <cmath>
typedef double db;
using namespace std;

const int MAXN = 105;

int n;
db a[MAXN][MAXN];

bool Gauss_Jordan()
{
	for (int i = 1; i <= n; i++)
	{
		int mx = i;
		for (int j = i + 1; j <= n; j++)
		{
			if (fabs(a[j][i]) < fabs(a[mx][i]))
			{
				mx = j;
			}
		}
		if (mx != i)
		{
			swap(a[i], a[mx]);
		}
		if (!a[i][i])
		{
			return false;
		}
		for (int j = 1; j <= n; j++)
		{
			if (i != j)
			{
				db val = a[j][i] / a[i][i];
				for (int k = i + 1; k <= n + 1; k++)
				{
					a[j][k] -= a[i][k] * val;
				}
			}
		}
	}
	return true;
}

int main()
{
	scanf("%d", &n);
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n + 1; j++)
		{
			scanf("%lf", a[i] + j);
		}
	}
	if (Gauss_Jordan())
	{
		for (int i = 1; i <= n; i++)
		{
			printf("%.2lf\n", a[i][n + 1] / a[i][i]);
		}
	}
	else
	{
		puts("No Solution");
	}
	return 0;
}

关于判断无解和无穷多解

P2455 [SDOI2006]线性方程组

这下要具体到到底是无解还是无穷多解了。

这就是高斯-约旦的一个缺点:相比于普通高斯,它更难判断无解和无穷多解。

但是也是可以处理的。

具体地,我们用 r 来记录当前行,如果主元为 0,那么 continue 到下一列,但 r 不变;否则消元后令 rr+1

遍历完所有列后:

  • r=n,说明有唯一解;
  • r<n,说明第 r+1n 行系数矩阵全都是 0,若列向量全是 0,说明有无穷多解,否则无解。

Code

#include <iostream>
#include <cstdio>
#include <cmath>
typedef double db;
using namespace std;

const int MAXN = 55;

int n;
db a[MAXN][MAXN], ans[MAXN];

int Gauss_Jordan()
{
	int r = 1;
	for (int i = 1; i <= n; i++)
	{
		int mx = r;
		for (int j = r + 1; j <= n; j++)
		{
			if (fabs(a[j][i]) > fabs(a[mx][i]))
			{
				mx = j;
			}
		}
		if (mx != r)
		{
			swap(a[r], a[mx]);
		}
		if (!a[r][i])
		{
			continue;
		}
		for (int j = 1; j <= n; j++)
		{
			if (j != r)
			{
				db val = a[j][i] / a[r][i];
				for (int k = i + 1; k <= n + 1; k++)
				{
					a[j][k] -= a[r][k] * val;
				}
			}
		}
		r++;
	}
	if (--r < n)
	{
		for (int i = r + 1; i <= n; i++)
		{
			if (a[i][n + 1])
			{
				return -1;
			}
		}
		return 0;
	}
	for (int i = 1; i <= n; i++)
	{
		ans[i] = a[i][n + 1] / a[i][i];
		if (!ans[i])
		{
			ans[i] = 0;
		}
	}
	return 1;
}

int main()
{
	scanf("%d", &n);
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n + 1; j++)
		{
			scanf("%lf", a[i] + j);
		}
	}
	int res = Gauss_Jordan();
	if (res != 1)
	{
		printf("%d\n", res);
	}
	else
	{
		for (int i = 1; i <= n; i++)
		{
			printf("x%d=%.2lf\n", i, ans[i]);
		}
	}
	return 0;
}
posted @   mango09  阅读(849)  评论(0编辑  收藏  举报
编辑推荐:
· DeepSeek 解答了困扰我五年的技术问题
· 为什么说在企业级应用开发中,后端往往是效率杀手?
· 用 C# 插值字符串处理器写一个 sscanf
· Java 中堆内存和栈内存上的数据分布和特点
· 开发中对象命名的一点思考
阅读排行:
· DeepSeek 解答了困扰我五年的技术问题。时代确实变了!
· PPT革命!DeepSeek+Kimi=N小时工作5分钟完成?
· What?废柴, 还在本地部署DeepSeek吗?Are you kidding?
· 赶AI大潮:在VSCode中使用DeepSeek及近百种模型的极简方法
· DeepSeek企业级部署实战指南:从服务器选型到Dify私有化落地
-->
点击右上角即可分享
微信分享提示