用高斯消元法求解多元线性方程组(C++)

用高斯消元法求解多元线性方程组

多元线性方程组

\qquad 形式为:
{ a 11 x 1 + a 12 x 2 + ⋯ + a 1 n x n = b 1 a 21 x 1 + a 22 x 2 + ⋯ + a 2 n x n = b 2 ⋮ a n 1 x 1 + a n 2 x 2 + ⋯ + a n n x n = b n \begin{cases} a_{11}x_1+a_{12}x_2+\dots+a_{1n}x_n=b_1 \\ a_{21}x_1+a_{22}x_2+\dots+a_{2n}x_n=b_2 \\ \vdots \\ a_{n1}x_1+a_{n2}x_2+\dots+a_{nn}x_n=b_n \\ \end{cases} a11x1+a12x2++a1nxn=b1a21x1+a22x2++a2nxn=b2an1x1+an2x2++annxn=bn
\qquad 求解 x 1 , x 2 , … , x n x_1,x_2,\dots,x_n x1,x2,,xn.

\qquad 解的可能一共有三种:
{ 1. 有唯一解 2. 无解 3. 有无穷多组解 \begin{cases} 1.有唯一解 \\ 2.无解 \\ 3.有无穷多组解 \\ \end{cases} 1.有唯一解2.无解3.有无穷多组解

解法:高斯消元法

初等行变换操作:

(以下操作是对原方程组的等价变形

  1. 把某一行乘一个非零的数;(等式两边同时乘一个非0的数)
  2. 交换某两行;(交换某两个方程的位置)
  3. 把某行的若干倍加到另一行上去。(某个方程乘若干倍加到另一个方程上去)

步骤:

枚举每一列:

  1. 找到剩下行中当前这一列绝对值最大的一行;
  2. 操作2: 将这一行换到最上面去;
  3. 操作1: 将该行的第一个数变成1,即每个系数都等比缩小(或放大);
  4. 操作3: 用这一行将下面所有行的当前列消成0。

举例:

1.有唯一解

{ − x 1 − x 2 + 2 x 3 = 7 2 x 1 + x 2 − 3 x 3 = − 9 x 1 + 2 x 2 − x 3 = − 6 \begin{cases} -x_1-x_2+2x_3=7 \\ 2x_1+x_2-3x_3=-9 \\ x_1+2x_2-x_3=-6 \\ \end{cases} x1x2+2x3=72x1+x23x3=9x1+2x2x3=6

写成矩阵的形式:

-1-127
21-3-9
12-1-6

对于第 1 列

  1. 先找到 2 2 2 最大;
  2. 将第 2 2 2 行换到第 1 1 1 行;
21-3-9
-1-127
12-1-6
  1. 使第一个数变成 1 1 1,同时除以 2 2 2
10.5-1.5-4.5
-1-127
12-1-6
  1. 用这一行将下面所有行的当前列消成 0 0 0:第 2 2 2 行减第 1 1 1 行的 − 1 -1 1 倍;第3行减第1行的 1 1 1 倍。
10.5-1.5-4.5
0-0.50.52.5
01.50.5-1.5

对于第 2 列

  1. 第 2 行开始找,第 2 列绝对值最大的是1.5;
  2. 将第 3 3 3 行换到第 2 行
10.5-1.5-4.5
01.50.5-1.5
0-0.50.52.5
  1. 使第一个数变成 1 1 1,同时除以 1.5 1.5 1.5
10.5-1.5-4.5
011/3-1
0-0.50.52.5
  1. 用这一行将下面所有行的当前列消成 0 0 0:第 3 3 3 行减第 1 1 1 行的 − 0.5 -0.5 0.5 倍。
10.5-1.5-4.5
011/3-1
002/32

对于第 3 列

  1. 第 3 行开始找,第 3 列绝对值最大 的是 2 3 \frac{2}{3} 32
  2. 第 3 行换到第 3 行,相当于不变;
  3. 使第 3 行第一个数变成 1 1 1,同时除以 2 3 \frac{2}{3} 32
10.5-1.5-4.5
011/3-1
0013
  1. 没有下一行,跳过此步。

\qquad 所有有唯一解的方程组,经过此 4 4 4 步的循环操作之后,行数 r r r 会枚举完每一行,最终行列式会变成上面的三角形的形式。也就是将方程组化简为这样:
{ x 1 + 0.5 x 2 − 1.5 x 3 = − 4.5 x 2 + 1 3 x 3 = − 1 x 3 = 3 \begin{cases} x_1+0.5x_2-1.5x_3=-4.5 \\ \qquad \qquad x_2+\frac{1}{3}x_3=-1 \\ \qquad \qquad \qquad \quad x_3=3 \\ \end{cases} x1+0.5x21.5x3=4.5x2+31x3=1x3=3
\qquad 我们就可以倒着一步步推出每个解。最终每个解就是矩阵每行的最后一个数字,即 x i = a [ i ] [ n + 1 ] x_i=a[i][n+1] xi=a[i][n+1]

2.无解

{ x 1 = 0 x 3 = 1 x 3 = 0 \begin{cases} x_1=0 \\ x_3=1 \\ x_3=0 \\ \end{cases} x1=0x3=1x3=0

其对应的矩阵表示为:

1000
0011
0010

\qquad 枚举到第 2 列时,从第 2 行看起,会发现绝对值最大的是0;

\qquad 于是直接从第 3 列第 2 行看起,找到绝对值最大的第 2 行的 1换到第 2 行不必交换了。然后把下面的全部变成 0 0 0,也就是用第 3 行减去第 2 行。结果为:

1000
0011
000-1

也就等价于:
{ x 1 = 0 x 3 = 1 0 = − 1 \begin{cases} x_1=0 \\ x_3=1 \\ 0=-1 \\ \end{cases} x1=0x3=10=1
\qquad 由此我们看出,经过等价变形,我们得出了一个显然不成立的式子 0 = − 1 0=-1 0=1,所以这个方程组无解。

\qquad 也就是说,当出现 0 = 非零数 0=非零数 0=非零数 这样形式的式子的时候,就是方程组无解的时候。

3.有无穷组解

{ x 1 + x 3 = 2 x 1 + x 2 = 3 2 x 1 + 2 x 2 = 6 \begin{cases} x_1+x_3=2 \\ x_1+x_2=3 \\ 2x_1+2x_2=6 \\ \end{cases} x1+x3=2x1+x2=32x1+2x2=6

\qquad 首先我们会发现,(3)式是由(2)式乘2得到,相当于这个方程组只有两个式子,却有3个未知数,显然有无穷组解。

\qquad 我们看看程序运行过程是怎样的:

1012
1103
2206

对于第 1 列

\qquad 首先交换1、3行,将 x 1 x_1 x1 的系数变为 1 1 1,变成: { x 1 + x 2 = 3 x 1 + x 2 = 3 x 1 + x 3 = 2 \begin{cases} x_1+x_2=3 \\ x_1+x_2=3 \\ x_1+x_3=2 \\ \end{cases} x1+x2=3x1+x2=3x1+x3=2

\qquad 然后将2、3行的 x 1 x_1 x1 项消掉,即(2)式减(1)式,(3)式减((1)式,结果为: { x 1 + x 2 = 3 0 = 0 − x 2 + x 3 = − 1 \begin{cases} x_1+x_2=3 \\ 0=0 \\ -x_2+x_3=-1 \\ \end{cases} x1+x2=30=0x2+x3=1,即:

1103
0000
0-11-1

对于第 2 列

\qquad 首先交换2、3行,将 x 2 x_2 x2 的系数变为 1 1 1,变成: { x 1 + x 2 = 3 x 2 − x 3 = 1 0 = 0 \begin{cases} x_1+x_2=3 \\ x_2-x_3=1 \\ 0=0 \\ \end{cases} x1+x2=3x2x3=10=0;下面的 x 2 x_2 x2 项系数全部为 0 0 0,因此无需改动。即:

1103
01-11
0000

\qquad 枚举的行数 r r r 没有枚举到最后一行就已经完毕,因此属于非正常情况,又没有出现 0=非零数 这样的无解情况,因此是无穷组解。

代码

#include <iostream>
#include <cmath>

using namespace std;

const int N = 105;
const double eps = 1e-6;

int n;
double a[N][N];

int gauss() {			//高斯消元法
	int r, c, i, j;							//r为行,c为列
	for (r = 1, c = 1; c <= n; c++) {
		int t = r;
		for (i = r; i <= n; i++)			//循环找到当前列绝对值最大的行t
			if (fabs(a[i][c]) > fabs(a[t][c]))
				t = i;
		if (fabs(a[t][c]) < eps)			//为0
			continue;

		if (t != r)							//有更大的才交换
			for (i = c; i <= n + 1; i++)	//交换到当前行,从c列开始交换是因为前面都是0
				swap(a[r][i], a[t][i]);
		if (a[r][c] != 1)					//不是1才需要变成1
			for (i = n + 1; i >= c; i--)	//将第c列的数字变为1,第一个数字要最后变
				a[r][i] /= a[r][c];

		for (i = r + 1; i <= n; i++)		//将下面每行的第c列变成0
			if (fabs(a[i][c]) > eps)		//不是0的才变,是0就不变了
				for (j = n + 1; j >= c; j--)
					a[i][j] -= a[r][j] * a[i][c];
		r++;
	}
	if (r <= n) {							//说明不是正常情况,正常r=n+1
		for (i = r; i <= n; i++)
			if (fabs(a[i][n + 1]) > eps)
				return -1;					//无解
		return 2;							//无穷组解
	}
	for (i = n; i > 0; i--)					//倒着一步步推出每个解
		for (j = i + 1; j <= n; j++)
			a[i][n + 1] -= a[i][j] * a[j][n + 1];
	return 1;								//有唯一解
}

int main() {
	int i, j;
	scanf("%d", &n);
	for (i = 1; i <= n; i++)
		for (j = 1; j <= n + 1; j++)			//输入的是一个n*(n+1)的矩阵
			scanf("%lf", &a[i][j]);
	int t = gauss();
	if (t == -1) {								//无解
		puts("No solution");
	} else if (t == 2) {						//有无穷组解
		puts("Infinite group solutions");
	} else {									//有唯一解
		for (i = 1; i <= n; i++)
			printf("%.2lf\n", a[i][n + 1]);
	}
	return 0;
}

变式

1.解异或线性方程组

形式为:
{ a 11 x 1 ∧ a 12 x 2 ∧ ⋯ ∧ a 1 n x n = b 1 a 21 x 1 ∧ a 22 x 2 ∧ ⋯ ∧ a 2 n x n = b 2 ⋮ a n 1 x 1 ∧ a n 2 x 2 ∧ ⋯ ∧ a n n x n = b n \begin{cases} a_{11}x_1 \land a_{12}x_2\land \dots \land a_{1n}x_n=b_1 \\ a_{21}x_1 \land a_{22}x_2 \land \dots \land a_{2n}x_n=b_2 \\ \vdots \\ a_{n1}x_1 \land a_{n2}x_2 \land \dots \land a_{nn}x_n=b_n \\ \end{cases} a11x1a12x2a1nxn=b1a21x1a22x2a2nxn=b2an1x1an2x2annxn=bn
其中 ^ 表示异或运算(XOR),方程组的系数 a a a 和常数 b b b 均为 0 0 0 1 1 1,未知数可能的值也为 0 0 0 1 1 1

求解: x 1 , x 2 , … , x n x_1,x_2,\dots,x_n x1,x2,,xn.

思路

与一般的高斯消元法相同,只是稍作改动。

  1. 不再需要 double类型,使用 bool 类型即可;
  2. 不再需要等比例缩小至系数为 1 1 1,因为系数不是 0 0 0 就是 1 1 1
  3. 加减消元可以用 ^ 运算来代替,因为异或运算俗称不进位的加法运算
  4. 乘法运算可以用 & 运算来代替,众所周知位运算是最快的运算。
举例

{ x 1 ∧ x 2 = 1 x 2 ∧ x 3 = 0 x 1 = 1 \begin{cases} x_1 \land x_2=1 \\ x_2 \land x_3=0 \\ x_1=1 \\ \end{cases} x1x2=1x2x3=0x1=1

表示为:

1101
0110
1001

将第 3 3 3 行的第 1 1 1 列变成 0 0 0,就是将第 1 1 1 行与第 3 3 3 行的系数做异或运算再放入第 3 3 3 行,变成: { x 1 ∧ x 2 = 1 x 2 ∧ x 3 = 0 x 2 = 0 \begin{cases} x_1 \land x_2=1 \\ x_2 \land x_3=0 \\ x_2=0 \\ \end{cases} x1x2=1x2x3=0x2=0

1101
0110
0100

将第 3 3 3 行的第 2 2 2 列变成 0 0 0,就是将第 2 2 2 行与第 3 3 3 行的系数做异或运算再放入第 3 3 3 行,变成: { x 1 ∧ x 2 = 1 x 2 ∧ x 3 = 0 x 3 = 0 \begin{cases} x_1 \land x_2=1 \\ x_2 \land x_3=0 \\ x_3=0 \\ \end{cases} x1x2=1x2x3=0x3=0

1101
0110
0010

最后倒推回去即可。这里注意:**异或运算的逆运算还是异或运算。**因此要求 x 2 x_2 x2,只需求 x 2 = 0 ∧ x 3 = 0 ∧ 0 = 0 x_2=0 \land x_3=0 \land 0=0 x2=0x3=00=0 即可。

最终可得到,此方程组的解为: { x 1 = 1 x 2 = 0 x 3 = 0 \begin{cases} x_1=1 \\ x_2=0 \\ x_3=0 \\ \end{cases} x1=1x2=0x3=0

代码
#include <iostream>

using namespace std;

const int N = 105;

int n;
bool a[N][N];


int gauss() {			//高斯消元法
	int r, c, i, j;							//r为行,c为列
	for (r = 1, c = 1; c <= n; c++) {
		int t = r;
		for (i = r; i <= n; i++)			//循环找到当前列绝对值最大的行t
			if (a[i][c]){					//是1就可以
				t = i;
				break;
			}	
		if (!a[t][c])	continue;			//为0

		if (t != r)							//有更大的才交换
			for (i = c; i <= n + 1; i++)	//交换到当前行,从c列开始交换是因为前面都是0
				swap(a[r][i], a[t][i]);

		for (i = r + 1; i <= n; i++)		//将下面每行的第c列变成0
			if (a[i][c])					//是1才需要变
				for (j = n + 1; j >= c; j--)
					a[i][j] ^= a[r][j];
		r++;
	}
	if (r <= n) {							//说明不是正常情况,正常r=n+1
		for (i = r; i <= n; i++)
			if (a[i][n + 1])
				return -1;					//无解
		return 2;							//多组解
	}
	for (i = n; i > 0; i--)					//倒着一步步推出每个解
		for (j = i + 1; j <= n; j++)
			a[i][n + 1] ^= a[i][j] & a[j][n + 1];
	return 1;								//有唯一解
}

int main() {
	int i, j;
	scanf("%d", &n);
	for (i = 1; i <= n; i++)
		for (j = 1; j <= n + 1; j++)			//输入的是一个n*(n+1)的矩阵
			cin >> a[i][j];						//此处注意,如果想用scanf读入,就不能使用bool数组了,应换成int
	int t = gauss();
	if (t == -1) {								//无解
		puts("No solution");
	} else if (t == 2) {						//有多组解
		puts("Multiple sets of solutions");
	} else {									//有唯一解
		for (i = 1; i <= n; i++)
			printf("%d\n", a[i][n + 1]);
	}
	return 0;
}
  • 2
    点赞
  • 5
    收藏
    觉得还不错? 一键收藏
  • 2
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值