应该算详细了
这是例题:
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 = 3 和 x_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;
}
显然,这代码判不了无解