浅谈高斯消元的实现和简单应用
一、高斯消元的原理
對于n元的m個線性方程組成的方程組,我們將其以矩陣的形式記錄下來:
a11 a12 a13 ...... a1n b1
a21 a22 a23 ...... a2n b2
...
...
...
an1 an2 an3 ...... ann bn
然后進行初等行列變換,嘗試構造出一個上三角矩陣,逐步使系數不為零的項減少;
等最后只剩下一個系數不為零時,進行回代,逐步求出已知解。(詳解過程咨詢小學老師)
二、高斯消元的實現
老老實實的回代代碼參見其他人的博客,這里介紹一種比較毒瘤的不回代暴力消元法:
1.Process
對于每個方程,按照一定的規則(后話)挑選一個主元,記錄該主元對應第幾個方程,然后用初等行列變換消去其他所有該元的系數;
最后我們得到的是一個每行只有一個數不為零,每列只有一個數不為零的鬼畜矩陣(自己腦補)
此時令ans向量對應的數字出去該行的非零系數,即為對應該元的解。
2.Code
設a數組為已經記錄系數的數組(格式見上方),那么a應該是n行n+1列的,最后一列代表方程的常數項;
設w數組記錄每個方程的主元是第幾項,v數組記錄答案;
當多解時輸出“Multiple solutions”,無解時輸出”No solution”;
double a[max_n][max_n+1],v[max_n];int w[max_n]; void gauss(){double eps=1e-6;for(int i=1;i<=n;++i){ //Enumerate the equation;int p=0; //Record the position of the largest number; double mx=0; //Recording the largest number;for(int j=1;j<=n;++j)if(fabs(a[i][j])-eps>mx){mx=fabs(a[i][j]);p=j; //fabs() returns the absolute value of float; }if(!p){if(fabs(a[i][n+1]<eps))printf("Multiple solutions");else printf("No solution");return; }w[i]=p;for(int j=1;j<=n;++j)if(i!=j){ //other equationsdouble t=a[j][p]/a[i][p];for(int k=1;k<=n+1;++k) //n+1 is importanta[j][k]-=a[i][k]*t;}}for(int i=1;i<=n;++i) v[w[i]]=a[i][n+1]/a[i][w[i]]; }3.notice
(1)精度的設置
眾所周知浮點數是有精度丟失的,在高斯消元中,精度的選擇要依題目而定,精度過低會導致系數較小的數被誤認為系數為零,而精度過高也有可能會導致誤差大于精度,導致本應該系數消為0的項誤認為系數不為零,所以精度的選擇是很哲學的問題,要注意。
推薦范圍:1e-4到1e-10
(2)主元的選取原則
接著(1)說,我們知道,用浮點數是有精度丟失的,那么讓一個較大的數除以一個極小的數誤差自然大的可想而知,所以我們想得到在精度允許的條件下系數最大的主元,所以對于每個方程,我們都應該選擇最大系數的元作為主元。
(3)在模2意義下的高斯消元
使用bitset優化運行時間,詳見相關應用中第三個例題的講解;
三、相關應用
這里給出高斯消元的幾道基礎題目,難度適合初學者。
1.[Luogu P3389]【模板】高斯消元
Description
給定一個線性方程組,對其求解
輸入格式:
第一行,一個正整數 n
第二至 n+1行,每行 n+1個整數,為 a1,a2?an和 b,代表一組方程。
輸出格式:
共n行,每行一個數,第 i行為 xi(保留2位小數)
如果不存在唯一解,在第一行輸出"No Solution".
Solution
如上所述。
Code
#include<iostream> #include<cstdio> #include<cmath> #include<cstring> #include<algorithm> using namespace std;const int max_n=110; double a[max_n][max_n+1],v[max_n]; int n,w[max_n]; inline int rd(){int x=0;bool f=0;char c=getchar();while(!isdigit(c)){if(c=='-') f=1;c=getchar();}while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=getchar();}return f?-x:x; }void gauss(){double eps=1e-6;for(int i=1;i<=n;++i){//enumerate the equation;int p=0; //Record the position of the largest number; double mx=0; //Recording the largest number;for(int j=1;j<=n;++j)if(fabs(a[i][j])-eps>mx){mx=fabs(a[i][j]);p=j;//fabs() returns the absolute value of float; }if(!p){printf("No Solution");return; }w[i]=p;for(int j=1;j<=n;++j)if(i!=j){ //other equationsdouble t=a[j][p]/a[i][p];for(int k=1;k<=n+1;++k)//n+1 is importanta[j][k]-=a[i][k]*t;}}for(int i=1;i<=n;++i) v[w[i]]=a[i][n+1]/a[i][w[i]];for(int i=1;i<=n;++i) printf("%.2lf\n",v[i]); }int main(){n=rd();for(int i=1;i<=n;++i)for(int j=1;j<=n+1;++j)a[i][j]=rd();gauss();return 0; }2.[BZOJ 1013][JSOI 2008] 球形空間產生器sphere
詳解參考我的隨筆:http://www.cnblogs.com/COLIN-LIGHTNING/p/8982341.html
轉載于:https://www.cnblogs.com/COLIN-LIGHTNING/p/8981923.html
總結
以上是生活随笔為你收集整理的浅谈高斯消元的实现和简单应用的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: Unity制作自适应透明背景(PC端)
- 下一篇: SCARA机器人与 DELTA机器人