洛谷 P3389:[模板] 高斯消元法 ← 高斯-约旦消元法(Gauss-Jordan Elimination)
【题目来源】
https://www.luogu.com.cn/problem/P3389
【题目描述】
给定一个线性方程组,对其求解。
【输入格式】
第一行,一个正整数 n。
第二至 n+1 行,每行 n+1 个整数,为 a1,a2,…,an 和 b,代表一组方程。
【输出格式】
共 n 行,每行一个数,第 i 行为 xi(四舍五入保留 2 位小数)。
如果不存在唯一解或无解,在第一行输出 No Solution。
【输入样例】
3
1 3 4 5
1 4 7 3
9 3 2 2
【输出样例】
-0.97
5.18
-2.39
【数据范围】
1≤n≤100,|ai|≤10^4,|b|≤10^4。
保证数据若有解则所有解均满足 |xi|≤10^3,且 xi±10^(-6) 和 xi 四舍五入后的结果相同(即不会因为较小的精度误差导致四舍五入后的结果不同)。
【算法分析】
● 高斯消元法(Gaussian Elimination) 是一种用于求解线性方程组的经典算法,也是线性代数中最基础的工具之一。它通过一系列初等变换,将增广矩阵化为行阶梯形矩阵或行最简形矩阵,从而求出方程组的解。
● 高斯消元法的本质是系统地消除未知数,将复杂的方程组转化为更容易求解的形式。整个过程分为两个阶段:
(1)前向消元(Forward Elimination):将矩阵化为上三角矩阵(行阶梯形)。
(2)回代求解(Back Substitution):从最后一个方程开始,逐个求解未知数。
● 核心代码解析
(一)选主元
int idx=r; for(int i=r; i<n; i++) { if(fabs(a[i][c])>fabs(a[idx][c])) idx=i; }从第 r 行开始,往下找第 c 列上绝对值最大的那个元素。
idx 记录绝对值最大的行号。
为什么要找绝对值最大的?—— 减小浮点误差。
(二)判断是否跳过该列?
if(fabs(a[idx][c])<eps) continue;如果找到的最大值都接近于 0,说明这一列全是 0。
这个未知数在当前剩余的方程中不存在,直接跳过该列(它可能是自由变量)。
continue 跳过后面的代码,进入下一列,但 r 不变(还没找到新的主元)。
(三)交换到当前行
swap(a[r],a[idx]);把找到的主元行(第 idx 行)换到第 r 行。
swap(a[r], a[idx]) 交换的是整行,包括常数项。
(四)归一化
double div=a[r][c]; for(int j=c; j<=n; j++) a[r][j]/=div;div 就是主元的值。
把第 r 行的所有元素(从第 c 列到第 n 列,含常数项)都除以 div。
结果:主元变成了 1,这一行其他元素相应缩放。
(五)消去其他行
for(int i=0; i<n; i++) { if(i!=r && fabs(a[i][c])>eps) { double val=a[i][c]; for(int j=c; j<=n; j++) { a[i][j]-=val*a[r][j]; } } }对所有行(不只是下面的行),如果第 c 列不为 0,就消掉它。
val = a[i][c] 是第 i 行第 c 列的当前值。
a[i][j] -= val * a[r][j] 就是用第 r 行(归一化后主元为 1)去消。
消完后,第 i 行的第 c 列变成 0。
注意:这里消的是所有行(包括上面的),这是高斯-约旦消元法的特点,直接得到对角矩阵,不用回代。
【算法代码】
#include <bits/stdc++.h> using namespace std; const int N=105; const double eps=1e-8; double a[N][N]; int n; int gauss() { int r=0; for(int c=0; c<n; c++) { //1.Select pivot int idx=r; for(int i=r; i<n; i++) { if(fabs(a[i][c])>fabs(a[idx][c])) idx=i; } if(fabs(a[idx][c])<eps) continue; swap(a[r],a[idx]); //2.Normalize the current row, and the pivot becomes 1 double div=a[r][c]; for(int j=c; j<=n; j++) a[r][j]/=div; //3.Eliminate all rows below for(int i=0; i<n; i++) { if(i!=r && fabs(a[i][c])>eps) { double val=a[i][c]; for(int j=c; j<=n; j++) { a[i][j]-=val*a[r][j]; } } } r++; } for(int i=r; i<n; i++) { //No solution if(fabs(a[i][n])>eps) return -1; } if(r<n) return 0; //Infinite solutions return 1; //Unique solution } int main() { cin>>n; for(int i=0; i<n; i++) { for(int j=0; j<=n; j++) { cin>>a[i][j]; } } int t=gauss(); if(t==-1 || t==0) cout<<"No Solution\n"; else { for(int i=0; i<n; i++) { printf("%.2lf\n",a[i][n]); } } return 0; } /* in: 3 1 3 4 5 1 4 7 3 9 3 2 2 out: -0.97 5.18 -2.39 */
【参考文献】
https://www.luogu.com.cn/problem/solution/P3389
