[2018.12.18]高斯消元
其实并没有想象中那么复杂。
就是一个简单的解方程过程。
众所周知,\(n\)个不等价的\(n\)元一次方程可以确定一组解(或者确定方程无解)。
高斯消元就是为了求得这组解。
以一个四元一次方程组
\(2x_1+x_2-5x_3+2x_4=12\)
\(-2x_1-x_2+x_3+x_4=14\)
\(x_1+6x_2-3x_3+x_4=20\)
\(3x_1+x_2+4x_3-7x_4=21\)
为例。
方便操作,我们将系数放入矩阵,常数项放在这个矩阵右侧。
\(\left[ \begin{array}{cccc} 2 & 1 & -5 & 2\\-2 & -1 & 1 & 1\\1 & 6 & -3 & 1\\3 & 1 & 1 & -1\\ \end{array}\right]\)\(\left[ \begin{array}{c}12\\14\\20\\21\end{array}\right]\)
首先,我们用第1行使后3行的第1个系数为0。
对于第2行,我们直接把它加上第1行;
对于第3行,我们把它加上第1行的\(-\frac{1}{2}\)倍;
对于第4行,我们把它加上第1行的\(-\frac{3}{2}\)倍;
于是我们得到了
\(\left[ \begin{array}{cccc} 2 & 1 & -5 & 2\\0 & 0 & -4 & 3\\0 & \frac{11}{2} & -\frac{1}{2} & 0\\0 & -\frac{1}{2} & \frac{17}{2} & -4\\ \end{array}\right]\)\(\left[ \begin{array}{c}12\\26\\14\\3\end{array}\right]\)
然后,我们用第2行使后2行的第2个系数为0。
但是我们发现第2行第2个系数本身为0。
这时我们要交换它和之后某一第2个系数不为0的行,如果不存在则方程无解/有无数组解。
因为第3行难算我们交换第2,4行,得到
\(\left[ \begin{array}{cccc} 2 & 1 & -5 & 2\\0 & -\frac{1}{2} & \frac{17}{2} & -4\\0 & \frac{11}{2} & -\frac{1}{2} & 0\\0 & 0 & -4 & 3\\ \end{array}\right]\)\(\left[ \begin{array}{c}12\\3\\14\\26\end{array}\right]\)
然后正常工作,矩阵变为
\(\left[ \begin{array}{cccc}2&1&-5&2\\0&-\frac{1}{2}&\frac{17}{2}&-4\\0&0&93&-44\\0&0&-4&3\end{array}\right]\)\(\left[ \begin{array}{c}12\\3\\47\\26\end{array}\right]\)
数没选好的后果即将出现
然后用第3行使最后一行的第3个系数变为0。
然后为了方便计算(没方便到哪里去),交换第3,4行。
\(\left[ \begin{array}{cccc}2&1&-5&2\\0&-\frac{1}{2}&\frac{17}{2}&-4\\0&0&-4&3\\0&0&93&-44\end{array}\right]\)\(\left[ \begin{array}{c}12\\3\\26\\47\end{array}\right]\)
然后...
\(\left[ \begin{array}{cccc}2&1&-5&2\\0&-\frac{1}{2}&\frac{17}{2}&-4\\0&0&-4&3\\0&0&0&\frac{147}{4}\end{array}\right]\)\(\left[ \begin{array}{c}12\\3\\26\\\frac{1303}{2}\end{array}\right]\)
好在我们结束了解方程。
根据第4行,我们知道了\(x_4=\frac{2606}{147}\);
知道了\(x_4\),根据第3行,我们知道了\(x_3\);
知道了\(x_3,x_4\),根据第3行,我们知道了\(x_2\);
知道了\(x_2,x_3,x_4\),根据第3行,我们知道了\(x_1\)。
由于太难算已放弃
解完了。
模板题代码:
#include<bits/stdc++.h>
using namespace std;
struct equation{
double n[110],x;
}e[110];
int n;
double ans[110];
void Work(int x,int y){
double v=-e[y].n[x]/e[x].n[x];
for(int i=1;i<=n;i++)e[y].n[i]+=e[x].n[i]*v;
e[y].x+=e[x].x*v;
}
void Solve(){
for(int i=1;i<=n;i++){
for(int j=i;j<=n;j++){
if(e[j].n[i]){
swap(e[j],e[i]);
break;
}
}
if(!e[i].n[i])puts("No Solution"),exit(0);
for(int j=i+1;j<=n;j++)Work(i,j);
}
ans[n]=e[n].x/e[n].n[n];
double k=0;
for(int i=n;i>=1;i--,k=0){
for(int j=i+1;j<=n;j++)k+=ans[j]*e[i].n[j];
ans[i]=(e[i].x-k)/e[i].n[i];
}
}
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++)scanf("%lf",&e[i].n[j]);
scanf("%lf",&e[i].x);
}
Solve();
for(int i=1;i<=n;i++)printf("%.2lf\n",ans[i]);
return 0;
}
练习题:
[JSOI2008]球形空间产生器
->Luogu
->BZOJ
->题解
