高斯消元详解

应该算详细了

这是例题:

P3389 【模板】高斯消元法 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)

给你一个这样的方程:

\begin{cases} a_{1, 1} x_1 + a_{1, 2} x_2 + \cdots + a_{1, n} x_n = b_1 \\ a_{2, 1} x_1 + a_{2, 2} x_2 + \cdots + a_{2, n} x_n = b_2 \\ \cdots \\ a_{n,1} x_1 + a_{n, 2} x_2 + \cdots + a_{n, n} x_n = b_n \end{cases}

你要求出

x_1, x_2, \cdots, x_n 的值

以这样一个方程为例:

\begin{cases} 2x_1 +4x_2 + x_3 = 13 \\ 3x_1 + 4 x_2 + 5x_3 = 26 \\ 2x_1 + 3x_2 + 9x_3 = 35\end{cases}

一步一步转化:

将第一个式子中的 x_1 项的系数化为 1

\begin{cases} x_1 +2x_2 + \frac{x_3}2 = \frac{13}2 \\ 3x_1 + 4 x_2 + 5x_3 = 26 \\ 2x_1 + 3 x_2 + 9x_3 = 35\end{cases}

消掉它下面式子中的 x_1

\begin{cases} x_1 +2x_2 + \frac{x_3}2 = \frac{13}2 \\ -2 x_2 + \frac{7x_3} 2 = \frac{13}2 \\ - x_2 + 8 x_3 = 22\end{cases}

将第二个式子中的 x_2 项的系数化为 1

\begin{cases} x_1 +2x_2 + \frac{x_3}2 = \frac{13}2 \\ x_2 + -\frac{7x_3} 4= -\frac{13}4 \\ - x_2 + 8 x_3 = 22\end{cases}

消掉它下面式子中的 x_2

\begin{cases} x_1 +2x_2 + \frac{x_3}2 = \frac{13}2 \\ x_2 + -\frac{7x_3} 4 = -\frac{13}4 \\ \frac{25x_3} 4 = \frac{75}4\end{cases}

现在可以发现第三个式子只剩下一个元了,求得

x_3 = 3

x_3 = 3 回带到第二个式子

x_2 + -\frac{21} 4 = -\frac{13}4

求得

x_2 = 2

x_3 = 3x_2 = 2 回带到第一个式子

x_1 +2x_2 + \frac{x_3}2 = \frac{13}2

求得

x_1 = 1

故原方程解为

\begin{cases} x_1 = 1 \\ x_2 = 2 \\ x_3 = 3 \end{cases}

虽然很麻烦,但对于计算机来说的确是行之有效的算法

核心就是

把原方程转化成这样

\begin{cases} a_{1, 1} x_1 + a_{1, 2} x_2 + a_{3, 1} x_3 + \cdots + a_{1, n} x_n = b_1 \\ 0 + a_{2, 2} x_2 + a_{3, 2} x_2 + \cdots + a_{2, n} x_n = b_2 \\ 0 + 0 +a_{3, 3} x_3+ \cdots + a_{3, n} x_n = b_3 \\ \cdots \\ 0 +0 + 0 + \cdots + a_{n, n} x_n = b_n \end{cases}

然后从最后一个式子回溯求解

时间复杂度 O(n^3)

这里我们可以再优化一下

让主元是最大的来避免误差

可以以看看这一篇 专栏 尝试用期望证明之(反正我不会)

代码

#include <iostream>
#include <math.h> 
using namespace std;
int n;
double a[1001][1002], ans[1001];
const double eps = -1e5;
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]);
	for(int i = 1; i <= n; i++)//枚举每个x项(从x1到xn)
	{
		int x = i;//优化:让主元是最大的来避免误差
		for (int j = i + 1; j <= n; j++)
			if (fabs(a[x][i]) < fabs(a[j][i]))
				x = j;
		if (fabs(a[x][i]) <= eps)//无解输出-1
		{
			printf("-1");
			return 0; 
		}
		if (x != i)//保证主元最大
			swap(a[x], a[i]);
		double d = a[i][i];
		for(int j = i; j <= n + 1; j++)//系数化为1
            a[i][j] /= d;
        for (int j = i + 1; j <= n; j++)//消掉它下面式子中的xi项
		{
			d = a[j][i];
			for (int k = i; k <= n + 1; k++)
				a[j][k] -= a[i][k] * d;
		}
	}
	ans[n] = a[n][n + 1];//xn的解
	for(int i = n - 1; i >= 1; i--)//回溯求x1到xn-1的解
	{
        ans[i] = a[i][n + 1];
        for(int j = i + 1; j <= n; j++)
            ans[i] -= a[i][j] * ans[j];
    }
    for(int i = 1; i <= n; i++)//输出解
        printf("x%d=%.2lf\n", i, ans[i]);
	return 0;
}

显然,这代码判不了无解

7 个赞

@邓君羽1 tj有些爆了

2 个赞

今天在学

已经疯了~