浅谈高斯消元的实现和简单应用
一、高斯消元的原理
對(duì)于n元的m個(gè)線性方程組成的方程組,我們將其以矩陣的形式記錄下來:
a11 a12 a13 ...... a1n b1
a21 a22 a23 ...... a2n b2
...
...
...
an1 an2 an3 ...... ann bn
然后進(jìn)行初等行列變換,嘗試構(gòu)造出一個(gè)上三角矩陣,逐步使系數(shù)不為零的項(xiàng)減少;
等最后只剩下一個(gè)系數(shù)不為零時(shí),進(jìn)行回代,逐步求出已知解。(詳解過程咨詢小學(xué)老師)
二、高斯消元的實(shí)現(xiàn)
老老實(shí)實(shí)的回代代碼參見其他人的博客,這里介紹一種比較毒瘤的不回代暴力消元法:
1.Process
對(duì)于每個(gè)方程,按照一定的規(guī)則(后話)挑選一個(gè)主元,記錄該主元對(duì)應(yīng)第幾個(gè)方程,然后用初等行列變換消去其他所有該元的系數(shù);
最后我們得到的是一個(gè)每行只有一個(gè)數(shù)不為零,每列只有一個(gè)數(shù)不為零的鬼畜矩陣(自己腦補(bǔ))
此時(shí)令ans向量對(duì)應(yīng)的數(shù)字出去該行的非零系數(shù),即為對(duì)應(yīng)該元的解。
2.Code
設(shè)a數(shù)組為已經(jīng)記錄系數(shù)的數(shù)組(格式見上方),那么a應(yīng)該是n行n+1列的,最后一列代表方程的常數(shù)項(xiàng);
設(shè)w數(shù)組記錄每個(gè)方程的主元是第幾項(xiàng),v數(shù)組記錄答案;
當(dāng)多解時(shí)輸出“Multiple solutions”,無(wú)解時(shí)輸出”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)精度的設(shè)置
眾所周知浮點(diǎn)數(shù)是有精度丟失的,在高斯消元中,精度的選擇要依題目而定,精度過低會(huì)導(dǎo)致系數(shù)較小的數(shù)被誤認(rèn)為系數(shù)為零,而精度過高也有可能會(huì)導(dǎo)致誤差大于精度,導(dǎo)致本應(yīng)該系數(shù)消為0的項(xiàng)誤認(rèn)為系數(shù)不為零,所以精度的選擇是很哲學(xué)的問題,要注意。
推薦范圍:1e-4到1e-10
(2)主元的選取原則
接著(1)說,我們知道,用浮點(diǎn)數(shù)是有精度丟失的,那么讓一個(gè)較大的數(shù)除以一個(gè)極小的數(shù)誤差自然大的可想而知,所以我們想得到在精度允許的條件下系數(shù)最大的主元,所以對(duì)于每個(gè)方程,我們都應(yīng)該選擇最大系數(shù)的元作為主元。
(3)在模2意義下的高斯消元
使用bitset優(yōu)化運(yùn)行時(shí)間,詳見相關(guān)應(yīng)用中第三個(gè)例題的講解;
三、相關(guān)應(yīng)用
這里給出高斯消元的幾道基礎(chǔ)題目,難度適合初學(xué)者。
1.[Luogu P3389]【模板】高斯消元
Description
給定一個(gè)線性方程組,對(duì)其求解
輸入格式:
第一行,一個(gè)正整數(shù) n
第二至 n+1行,每行 n+1個(gè)整數(shù),為 a1,a2?an和 b,代表一組方程。
輸出格式:
共n行,每行一個(gè)數(shù),第 i行為 xi(保留2位小數(shù))
如果不存在唯一解,在第一行輸出"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] 球形空間產(chǎn)生器sphere
詳解參考我的隨筆:http://www.cnblogs.com/COLIN-LIGHTNING/p/8982341.html
轉(zhuǎn)載于:https://www.cnblogs.com/COLIN-LIGHTNING/p/8981923.html
總結(jié)
以上是生活随笔為你收集整理的浅谈高斯消元的实现和简单应用的全部?jī)?nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: Unity制作自适应透明背景(PC端)
- 下一篇: SCARA机器人与 DELTA机器人