洛谷P3389 【模板】高斯消元法
P3389 【模板】高斯消元法
题目背景
Gauss消元
题目描述
给定一个线性方程组,对其求解
输入输出格式
输入格式:
第一行,一个正整数 n
第二至 n+1行,每行 n+1 个整数,为a1,a2⋯an 和 b,代表一组方程。
输出格式:
共n行,每行一个数,第 i行为 xi (保留2位小数)
如果不存在唯一解,在第一行输出"No Solution".
输入输出样例
输入样例#1:
3
1 3 4 5
1 4 7 3
9 3 2 2
输出样例#1:
-0.97
5.18
-2.39
说明
1≤n≤100,∣ai∣≤104,∣b∣≤104
sol:究竟是什么东西(感觉像小学奥数)
我们首先确定一个方程组作为例子
x-2y+3z=6
4x-5y+6z=12
7x-8y+10z=21
先将它转化为矩阵
1 -2 3 6
4 -5 6 12
7 -8 10 21
解决这个方程组
我们会希望它变成如下形式
1 0 0 a
0 1 0 b
0 0 1 c
这样就可以表示为x=a,y=b,z=cx=a,y=b,z=c
我们使用高斯消元,就要一步一步将每个未知数约去。
这种方法是以列为单位消去的
首先我们将第一列转化为1 0 0的形式
在这里要注意一下,我们往往是将这个系数绝对值最大的方程转移到被减的这一行,这样就可以减小误差
所以我们先将矩阵变成这样
7 -8 10 21
4 -5 6 12
1 -2 3 6
然后将正在处理的方程式化简,让正被处理的系数化1
1 -8/7 10/7 3
4 -5 6 12
1 -2 3 6
然后使用加减法将第二个与第三个方程组的第一个系数化0
1 -8/7 10/7 3
0 -3/7 2/7 0
0 -6/7 11/7 3
然后这时候第一列就被化简完成
同理我们化去第二行与第三行,步骤如下:
1.化简第二行
1 -8/7 10/7 3
0 1 -2/3 0
0 -6/7 11/7 3
2.用第一行减第二行×(-8/7),第三行减第二行×(-6/7)
1 0 2/3 3
0 1 -2/3 0
0 0 1 3
3.不需要化简第三行,所以直接用第一行减去第三行×2/3,第二行减去第三行×(-2/3)
1 0 0 1
0 1 0 2
0 0 1 3
最后我们就得到了一组解x=1,y=2,y=3x=1,y=2,y=3。所以高斯消元其实是运用了小学解方程组的加减法的呢。
注意是一列列消去的
#include <bits/stdc++.h> using namespace std; typedef int ll; inline ll read() { ll S=0; bool f=0; char ch=' '; while(!isdigit(ch)) { f|=(ch=='-'); ch=getchar(); } while(isdigit(ch)) { S=S*10+(ch-'0'); ch=getchar(); } return (f)?(-S):(S); } #define R(x) x=read() inline void write(ll x) { if(x<0) { putchar('-'); x=-x; } if(x<10) { putchar(x+'0'); return; } write(x/10); putchar(x%10+'0'); return; } #define W(x) write(x),putchar(' ') #define Wl(x) write(x),putchar('\n') const int N=105; const double eps=1e-8; int n; double a[N][N]; inline bool Gauss() { int i,j,k; for(i=1;i<=n;i++) { k=i; for(j=i+1;j<=n;j++) if(fabs(a[j][i])>fabs(a[k][i])) k=j; //找到最大的数 if(fabs(a[k][i])<eps) return false; if(i!=k) for(j=i;j<=n+1;j++) swap(a[k][j],a[i][j]); //对换一行或一列,属于找最大当前系数的其中一步.(这样就可以只处理当前行的系数啦!) double Div=a[i][i]; for(j=i;j<=n+1;j++) a[i][j]/=Div; for(j=1;j<=n;j++) if(j!=i) { Div=a[j][i]; for(k=i;k<=n+1;k++) { a[j][k]-=a[i][k]*Div; } } //把这列除这个数外都消成0 } return true; } int main() { int i,j; R(n); for(i=1;i<=n;i++) { for(j=1;j<=n+1;j++) scanf("%lf",&a[i][j]); } if(!(Gauss())) { puts("No Solution"); } else { for(i=1;i<=n;i++) printf("%.2lf\n",a[i][n+1]); } return 0; } /* input 3 1 3 4 5 1 4 7 3 9 3 2 2 output -0.97 5.18 -2.39 */
附上更加符合上面教程的代码(完全按照上面写的)
#include <bits/stdc++.h> using namespace std; typedef int ll; inline ll read() { ll s=0; bool f=0; char ch=' '; while(!isdigit(ch)) { f|=(ch=='-'); ch=getchar(); } while(isdigit(ch)) { s=(s<<3)+(s<<1)+(ch^48); ch=getchar(); } return (f)?(-s):(s); } #define R(x) x=read() inline void write(ll x) { if(x<0) { putchar('-'); x=-x; } if(x<10) { putchar(x+'0'); return; } write(x/10); putchar((x%10)+'0'); return; } #define W(x) write(x),putchar(' ') #define Wl(x) write(x),putchar('\n') const double eps=1e-9; const int N=105; int n; double a[N][N],b[N]; inline void Debug() { int i,j; for(i=1;i<=n;i++) { for(j=1;j<=n;j++) printf("%.2lf ",a[i][j]); printf("%.2lf",b[i]); puts(""); } puts(""); } inline void Gauss(int n) { int i,j,k; double Div; for(i=1;i<=n;i++) { int Pos=i; for(j=i+1;j<=n;j++) if(fabs(a[Pos][i])<fabs(a[j][i])) Pos=j; if(fabs(a[Pos][i])<eps) { puts("No Solution"); exit(0); } if(Pos!=i) { swap(a[i],a[Pos]); swap(b[i],b[Pos]); } Div=a[i][i]; for(j=i;j<=n;j++) a[i][j]/=Div; b[i]/=Div; for(j=1;j<=n;j++) if(j!=i) { Div=a[j][i]; for(k=i;k<=n;k++) { a[j][k]-=Div*a[i][k]; } b[j]-=Div*b[i]; } // Debug(); } } int main() { int i,j; R(n); for(i=1;i<=n;i++) { for(j=1;j<=n;j++) scanf("%lf",&a[i][j]); scanf("%lf",&b[i]); } Gauss(n); for(i=1;i<=n;i++) printf("%.2lf\n",b[i]); return 0; } /* input 3 1 3 4 5 1 4 7 3 9 3 2 2 output -0.97 5.18 -2.39 */
河田は河田、赤木は赤木……。
私は誰ですか。教えてください、私は誰ですか。
そうだ、俺はあきらめない男、三井寿だ!