數(shù)值分析老師布置了這樣的實習作業(yè)：

隨機生成一個 $n$ 階方陣與 $n$ 階向量構(gòu)成 $Ax=b$

• 至少生成 $5$ 個方程并取不同的 $n$
• 編程實現(xiàn)克萊姆法則求 $x$
• 編程實現(xiàn)高斯消去法基本方法求 $x$
• 編程實現(xiàn)高斯消去法的 $LU$ 分解求 $x$
• 編程實現(xiàn)列主元高斯消去法求 $x$
• 列表比較并分析以上方法的運算時間與效率

看到這個作業(yè)后，我內(nèi)心其實是拒絕的。要是用 C++ 的話，估計得寫好久。

但作業(yè)也不能不交是吧，于是就寫了下面的代碼。

#include<bits/stdc++.h> using namespace std; const int N=110; const double eps=1e-6; int n,flag,bg,ed; double d[N],a[N][N],acp[N][N],l[N][N],u[N][N];void randinit(){ for(int i=1;i<=n;i++) for(int j=1;j<=n+1;j++) acp[i][j]=(rand()%200)/10.0-10; } void init(){ for(int i=1;i<=n;i++) for(int j=1;j<=n+1;j++) a[i][j]=acp[i][j]; } void out(){ printf("x = [ "); for(int i=1;i<=n;i++) printf("%lf ",a[i][n+1]); printf("]\n"); } void pt(double a[N][N]){ for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) printf("%lf%c",a[i][j]," \n"[j==n]); }void calc_U(){ //求解上三角方程組for(int i=n;i>=1;i--){for(int j=n;j>i;j--) a[i][n+1]-=a[i][j]*a[j][n+1];a[i][n+1]/=a[i][i];} }void calc_D(){ //求解下三角方程組for(int i=1;i<=n;i++){for(int j=1;j<i;j++) a[i][n+1]-=a[i][j]*a[j][n+1];a[i][n+1]/=a[i][i];} }void divide(){ //將A分解為L和U;for(int i=1;i<=n;i++){if(abs(a[i][i])<eps) return (void)(flag=true);for(int k=i+1;k<=n;k++){l[k][i]=a[k][i]/a[i][i];if(abs(a[k][i])>eps) for(int j=n+1;j>=i;j--)a[k][j]-=a[i][j]*(a[k][i]/a[i][i]);}}for(int i=1;i<=n;i++) for(int j=i;j<=n;j++) l[i][j]=i==j;for(int i=1;i<=n;i++) for(int j=i;j<=n;j++) u[i][j]=a[i][j]; }void gauss_base(){ //高斯消去法（基本）divide(); //化解為上三角calc_U(); //回代 }void gauss(){ //高斯消去法（列主元、無回代）for(int i=1;i<=n;i++){int mx=i;for(int j=i+1;j<=n;j++) if(a[j][i]>a[mx][i]) mx=j;if(abs(a[mx][i])<eps) return (void)(flag=true);for(int j=i;j<=n+1;j++) swap(a[i][j],a[mx][j]);for(int j=n+1;j>=i;j--) a[i][j]/=a[i][i];for(int k=1;k<=n;k++) if(i!=k)for(int j=n+1;j>=i;j--) a[k][j]-=a[k][i]*a[i][j];} }double det(){ //計算detdouble ret=1;for(int i=1;i<=n;i++){int mx=i;for(int j=i+1;j<=n;j++) if(a[j][i]>a[mx][i]) mx=j;if(abs(a[mx][i])<eps){ flag=true; return 0; }if(mx!=i) ret=-ret;ret*=a[mx][i];for(int j=i;j<=n+1;j++) swap(a[i][j],a[mx][j]);for(int j=n+1;j>=i;j--) a[i][j]/=a[i][i];for(int k=1;k<=n;k++) if(i!=k)for(int j=n+1;j>=i;j--) a[k][j]-=a[k][i]*a[i][j];}return ret; }void cramer(){ //克萊姆法則求解for(int i=0;i<=n;i++){init();if(i!=0) for(int j=1;j<=n;j++) swap(a[j][i],a[j][n+1]);d[i]=det();}for(int i=1;i<=n;i++) a[i][n+1]=d[i]/d[0]; }void LU(){ //LU分解法divide();for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) a[i][j]=l[i][j];for(int i=1;i<=n;i++) a[i][n+1]=acp[i][n+1];calc_D();for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) a[i][j]=u[i][j];calc_U(); }int main(){srand(time(0));printf("input n (2<=100) :\n");scanf("%d",&n);do{flag=false;randinit();init();divide();}while(flag); //生成有唯一解的方程組init();bg=clock(); gauss_base(); ed=clock();printf("gauss_base ans is : \n"); out();printf("gauss_base run time is %d ms\n\n",ed-bg);init();bg=clock(); gauss(); ed=clock();printf("gauss ans is : \n"); out();printf("gauss run time is %d ms\n\n",ed-bg);init();bg=clock(); LU(); ed=clock();if(n<=10){printf("L is :\n"); pt(l);printf("u is :\n"); pt(u);}printf("LU ans is : \n"); out();printf("LU run time is %d ms\n\n",ed-bg);bg=clock(); cramer(); ed=clock();printf("cramer ans is : \n"); out();printf("cramer run time is %d ms\n\n",ed-bg); }

當 $n$ 為 $100$ 時，四種方法的運行時間如下表：

基本Gauss列主元GaussLU分解法克萊姆法則

3ms	4ms	3ms	632ms

克萊姆法則的復雜度為 $O(n^4)$，明顯慢于其他 $O(n^3)$ 的算法。

剛開始很不理解，$LU$ 分解法的動作明顯比高斯消元更冗余。然后搜了搜它們的優(yōu)缺點，如下：

$LU$ 和 $Gauss$ 都是 $O(n^3)$ 。更精確地講，乘除法大概都是 $\frac{n^3}{3}$ 次，時間上差別不大，不過：
$LU$ 具有承襲性，這是 $LU$ 的優(yōu)點。
$LU$ 只適用于解所有順序主子式都大于 $0$ 的，通用性欠缺，這是 $LU$ 的缺點。
$LU$ 法不保證具有數(shù)值穩(wěn)定性，這是 $LU$ 的缺點。（ $Gauss$ 法可以用選取列主元技巧保證數(shù)值穩(wěn)定性）
集合 $LU$ 與 $Gauss$ 優(yōu)點，同時規(guī)避掉這些缺點的，是 $LUP$ 分解法。

參考鏈接 LU分解法與Gauss消元法兩者的比較。

所謂承襲性，就是說當計算 $Ax=y_i$ ($1\le i\le n$) 時，若使用 $Gauss$ ，那么只能依次調(diào)用 $n$ 次，總復雜度為 $O(k\times n^3)$。而使用 $LU$ 分解法時，先用 $O(n^3)$ 的復雜度將原式轉(zhuǎn)化為 $LUx=y_i$ ，之后進行 $k$ 次回代就能求出所有的 $x_i$ ，總復雜度為 $O(n^3+k\times n^2)$。

另外，在調(diào)試代碼的時候，我最初將數(shù)組開的很大，試了試 $200\times 200$ 的矩陣，結(jié)果其他方法都能正確求解，唯獨克萊姆法則輸出全為 $nan$ 。而對于 $100\times 100$ 的矩陣，四種方法都能正確求解。
后來發(fā)現(xiàn)，使用克萊姆法則，需要求解行列式的值。而我隨機生成的矩陣，這個行列式的值會隨著矩陣的規(guī)模增大而急劇增大，溢出 $double$。C++黨已經(jīng)輸麻了

總結(jié)

以上是生活随笔為你收集整理的高斯消元法、LU分解法与克莱姆法则解方程组的C++实现的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。