「學習筆記」高斯消元
點擊查看目錄
目錄
- 「學習筆記」高斯消元
- 演算法
- 思路
- 代碼
- 例題:[SDOI2006]線性方程組
- 思路
- 代碼
- 練習題
- [JSOI2008]球形空間產生器
- 思路
- 代碼
- [USACO10HOL]Driving Out the Piggies G
- 思路
- 代碼
注意:本文內的高斯消元均是「高斯-約旦消元法(\(\text{Gauss-Jordan Elimination}\))」,
演算法
思路
我們舉個例子:
\[\begin{cases} x_2+x_3&=3\\ x_1+2x_3&=9\\ 2x_1+x_2+x_3&=13 \end{cases} \]我們把它的所有系數化為矩陣 \(A\):
\[\left[\begin{matrix} 0&1&1\\ 1&0&2\\ 2&1&1\\ \end{matrix} \middle| \begin{matrix} 3 \\ 9 \\ 13 \end{matrix} \right] \]然后開始推:
- 首先遍歷第 \(i\) 列的第 \(i\) 行至 \(n\) 行(因為前 \(i-1\) 行已經形成陣型了,不能破壞),找出這一列最大的系數,并將那一行交換至第 \(i\) 行,(當前 \(i=1\))
-
然后判斷解,即判斷 \(A_{i,i}\) 是否為零,若為零則無唯一解,否則有唯一解,(這個例子顯然有唯一解)
-
系數化為一,即將 \(A_{i,j}(i<j\le n+1)\) 除以 \(A_{i,i}\),再將 \(A_{i,i}\) 除以自己得 \(1\),
- 然后用加減消元法消掉一部分元:
- 重復上述步驟:
所有操作步驟完成,最后得到解:
\[\begin{cases} x_1=5\\ x_2=1\\ x_3=2\\ \end{cases} \]時間復雜度 \(\Theta(n^3)\),
注意:輸入不要用快讀,雖然有的題輸入只有整數,但很多題都不是!
代碼
點擊查看代碼
inline void Gauss(){
_for(i,1,n){
ll mx=i;
_for(j,i+1,n)if(fabs(a[j][i])>fabs(a[mx][i]))mx=j;
_for(j,1,n+1)swap(a[i][j],a[mx][j]);
if(fabs(a[i][i])<eps){puts("No Solution");return;}
_for(j,i+1,n+1)a[i][j]/=a[i][i];
a[i][i]=1;
_for(j,1,n){
if(i==j)continue;
_for(k,i+1,n+1)a[j][k]-=a[j][i]*a[i][k];
}
}
_for(i,1,n)printf("%.2lf\n",a[i][n+1]);
}
例題:[SDOI2006]線性方程組
思路
本題主要難點在于如何區分無解和無限組解,
舉個例子:
\[A= \left[\begin{matrix} 1&0&1\\ 0&2&1\\ 1&2&2\\ \end{matrix}\middle|\begin{matrix} 2 \\ 3 \\ 5 \\ \end{matrix}\right], B= \left[\begin{matrix} 1&0&1\\ 0&2&1\\ 1&2&2\\ \end{matrix}\middle|\begin{matrix} 2 \\ 3 \\ 7 \\ \end{matrix}\right], C= \left[\begin{matrix} 1&0&1\\ 0&2&1\\ 1&1&2\\ \end{matrix}\middle|\begin{matrix} 2 \\ 3 \\ 5 \\ \end{matrix}\right] \]-
在 \(A\) 中,第三行的前 \(3\) 項是前兩行相加的結果,第 \(4\) 項也是,所以有用的行只有兩個(這個數量稱為秩),有一個數是不確定的,所以有不唯一解,
-
在 \(B\) 中,第三行的前 \(3\) 項是前兩行相加的結果,但第 \(4\) 項不是,會出現 \(0=2\) 的情況,所以無解,
-
在 \(C\) 中,它的秩是 \(3\),所以有唯一解:
歸納一下:
-
若秩不等于未知數的數量且消元后出現 \(\sum_{j=1}^{n}a_{i,j}=0\) 且 \(a_{i,n+1}\neq0\) 的情況,則無解,
-
若秩不等于未知數的數量且消元后出現 \(\sum_{j=1}^{n+1}a_{i,j}=0\) 的情況,則有無陣列解,
-
否則有唯一解,
代碼
點擊查看代碼
const ll N=110,inf=1ll<<40;
const double eps=1e-9;
ll n,l=1;db a[N][N];
namespace SOLVE{
inline void Gauss(){
_for(i,1,n){
ll mx=l;
_for(j,l+1,n)if(fabs(a[j][i])>fabs(a[mx][i]))mx=j;
if(fabs(a[mx][i])<eps)continue;
swap(a[mx],a[l]);
_for(j,i+1,n+1)a[l][j]/=a[l][i];
a[l][i]=1;
_for(j,1,n){
if(j==l)continue;
_for(k,i+1,n+1)a[j][k]-=a[l][k]*a[j][i];
}
++l;
}
}
inline void Jie(){
if(l<=n){
_for(i,l,n)if(fabs(a[i][n+1])>eps){puts("-1");return;}
puts("0");return;
}
_for(i,1,n){
printf("x%lld=%.2lf\n",i,a[i][n+1]);
}
}
inline void In(){
scanf("%lld",&n);
_for(i,1,n)_for(j,1,n+1)scanf("%lf",&a[i][j]);
Gauss();
Jie();
return;
}
}
練習題
[JSOI2008]球形空間產生器
思路
設球心為 \(p(p_1,p_2,\cdots,p_n)\),可列式子:
\[\begin{aligned} r=&(a_{1,1}-p_1)^2+(a_{1,2}-p_2)^2+\cdots+(a_{1,n}-p_n)^2\\ =&(a_{2,1}-p_1)^2+(a_{2,2}-p_2)^2+\cdots+(a_{2,n}-p_n)^2\\ =&\cdots\\ =&(a_{n+1,1}-p_1)^2+(a_{n+1,2}-p_2)^2+\cdots+(a_{n+1,n}-p_n)^2\\ \end{aligned} \]拆一下得:
\[\begin{aligned} r=&(a_{1,1}^2-2a_{1,1}p_1+p_1^2)+(a_{1,2}^2-2a_{1,2}p_2+p_2^2)+\cdots+(a_{1,n}^2-2a_{1,n}p_n+p_n^2)\\ =&(a_{2,1}^2-2a_{2,1}p_1+p_1^2)+(a_{2,2}^2-2a_{2,2}p_2+p_2^2)+\cdots+(a_{2,n}^2-2a_{2,n}p_n+p_n^2)\\ =&\cdots\\ =&(a_{n+1,1}^2-2a_{n+1,1}p_1+p_1^2)+(a_{n+1,2}^2-2a_{n+1,2}p_2+p_2^2)+\cdots+(a_{n+1,n}^2-2a_{n+1,n}p_n+p_n^2)\\ \end{aligned} \]化簡得:
\[\begin{aligned} r=&\sum_{i=1}^{n}a_{1,i}^2-\sum_{i=1}^{n}2a_{1,i}p_i\\ =&\sum_{i=1}^{n}a_{2,i}^2-\sum_{i=1}^{n}2a_{2,i}p_i\\ =&\cdots\\ =&\sum_{i=1}^{n}a_{n+1,i}^2-\sum_{i=1}^{n}2a_{n+1,i}p_i\\ \end{aligned} \]即:
\[\begin{aligned} \sum_{i=1}^{n}a_{1,i}^2=&r+\sum_{i=1}^{n}2a_{1,i}p_i\\ \sum_{i=1}^{n}a_{2,i}^2=&r+\sum_{i=1}^{n}2a_{2,i}p_i\\ \cdots=&r+\cdots\\ \sum_{i=1}^{n}a_{n+1,i}^2=&r+\sum_{i=1}^{n}2a_{n+1,i}p_i\\ \end{aligned} \]發現就是個方程,但是這個 \(r\) 混在里面很難受,那么我們把它當成一個系數為 \(1\) 的未知數,則原式可以化成如下增廣矩陣:
\[\left[ \begin{matrix} 2a_{1,1}&2a_{1,2}&\cdots&2a_{1,n}&1\\\\ 2a_{2,1}&2a_{2,2}&\cdots&2a_{2,n}&1\\\\ \cdots&\cdots&\cdots&\cdots&1\\\\ 2a_{n+1,1}&2a_{n+1,2}&\cdots&2a_{n+1,n}&1 \end{matrix} \middle| \begin{matrix} \sum_{i=1}^{n}a_{1,i}^2\\\\ \sum_{i=1}^{n}a_{2,i}^2\\\\ \cdots\\\\ \sum_{i=1}^{n}a_{n+1,i}^2\\ \end{matrix} \right] \]就可以直接高斯消元了,
Latex 當算草紙真好用,
代碼
點擊查看代碼
const ll N=20,inf=1ll<<40;
const double eps=1e-6;
ll n;ldb a[N][N];
namespace SOLVE{
inline ll rnt(){
ll x=0,w=1;char c=getchar();
while(!isdigit(c)){if(c=='-')w=-1;c=getchar();}
while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x*w;
}
inline void Gauss(){
_for(i,1,n){
ll mx=i;
_for(j,i+1,n)if(fabs(a[j][i])>fabs(a[mx][i]))mx=j;
swap(a[mx],a[i]);
_for(j,i+1,n+1)a[i][j]/=a[i][i];
a[i][i]=1.0;
_for(j,1,n){
if(i==j)continue;
for_(k,n+1,i)a[j][k]-=a[j][i]*a[i][k];
}
}
}
inline void In(){
scanf("%lld",&n),++n;
_for(i,1,n){
_for(j,1,n-1){
scanf("%Lf",&a[i][j]);
a[i][n+1]+=a[i][j]*a[i][j];
a[i][j]*=2.0;
}
a[i][n]=1;
}
Gauss();
_for(i,1,n-1)printf("%.3Lf ",a[i][n+1]);
return;
}
}
[USACO10HOL]Driving Out the Piggies G
思路
好題,還對高斯消元能出什么題沒什么認知就 重繪了我對高斯消元能出什么題的認知,
我們設到節點 \(i\) 的期望步數為 \(x_i\),節點 \(i\) 在圖上的出度為 \(d_i\),
那么顯然:
\[x_i=[i=1]+\sum_{(i,j)\in\mathbb{E}}x_j\times(1-\frac{p}{q})\times d_j \]顯然 \(x_i\times\frac{p}{q}\) 為實際爆炸的概率,那么考慮如何求 \(x_i\),我們進行一個移項:
\[x_i-\sum_{(i,j)\in\mathbb{E}}x_j\times(1-\frac{p}{q})\times d_j=[i=1] \]直接沖一個高斯消元即可,
代碼
點擊查看代碼
const ll N=310,inf=1ll<<40;
const double eps=1e-13;
ll n,m,tu[N][N];ldb p,q,d[N],a[N][N];
namespace SOLVE{
inline ll rnt(){
ll x=0,w=1;char c=getchar();
while(!isdigit(c)){if(c=='-')w=-1;c=getchar();}
while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x*w;
}
inline void Pre(){
a[1][n+1]=1;
_for(i,1,n){
_for(j,1,n){
if(i==j)a[i][j]=1.0;
else if(tu[i][j])a[i][j]=-(1.0-p/q)/d[j];
}
}
}
inline void Gauss(){
_for(i,1,n){
ll mx=i;
_for(j,i+1,n)if(a[mx][i]<a[j][i])mx=j;
swap(a[mx],a[i]);
_for(j,i+1,n+1)a[i][j]/=a[i][i];
a[i][i]=1;
_for(j,1,n){
if(i==j)continue;
for_(k,n+1,i)a[j][k]-=a[i][k]*a[j][i];
}
}
}
inline void In(){
scanf("%lld%lld%Lf%Lf",&n,&m,&p,&q);
_for(i,1,m){
ll x=rnt(),y=rnt();
d[x]+=1.0,d[y]+=1.0;
tu[x][y]=tu[y][x]=1;
}
Pre(),Gauss();
_for(i,1,n)printf("%.9Lf\n",a[i][n+1]*p/q);
return;
}
}
本文來自博客園,作者:Keven-He,轉載請注明原文鏈接:https://www.cnblogs.com/Keven-He/p/GaussElimination.html
轉載請註明出處,本文鏈接:https://www.uj5u.com/qita/500769.html
標籤:其他
上一篇:leetcode 448. Find All Numbers Disappeared in an Array 找到所有陣列中消失的數字(簡單)
下一篇:圖論-虛擬節點分層建圖
