高斯消元法、LU分解法与克莱姆法则解方程组的C++实现
數值分析老師布置了這樣的實習作業:
隨機生成一個 𝑛𝑛n 階方陣與 𝑛𝑛n 階向量構成 𝐴𝑥=𝑏𝐴𝑥 = 𝑏Ax=b
- 至少生成 555 個方程并取不同的 𝑛𝑛n
- 編程實現克萊姆法則求 𝑥𝑥x
- 編程實現高斯消去法基本方法求 𝑥𝑥x
- 編程實現高斯消去法的 LULULU 分解求 𝑥𝑥x
- 編程實現列主元高斯消去法求 𝑥𝑥x
- 列表比較并分析以上方法的運算時間與效率
看到這個作業后,我內心其實是拒絕的。要是用 C++ 的話,估計得寫好久。
但作業也不能不交是吧,于是就寫了下面的代碼。
#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); }當 nnn 為 100100100 時,四種方法的運行時間如下表:
| 3ms | 4ms | 3ms | 632ms |
克萊姆法則的復雜度為 O(n4)O(n^4)O(n4),明顯慢于其他 O(n3)O(n^3)O(n3) 的算法。
剛開始很不理解,LULULU 分解法的動作明顯比高斯消元更冗余。然后搜了搜它們的優缺點,如下:
LULULU 和 GaussGaussGauss 都是 O(n3)O(n^3)O(n3) 。更精確地講,乘除法大概都是 n33\frac{n^3}{3}3n3? 次,時間上差別不大,不過:
LULULU 具有承襲性,這是 LULULU 的優點。
LULULU 只適用于解所有順序主子式都大于 000 的,通用性欠缺,這是 LULULU 的缺點。
LULULU 法不保證具有數值穩定性,這是 LULULU 的缺點。( GaussGaussGauss 法可以用選取列主元技巧保證數值穩定性)
集合 LULULU 與 GaussGaussGauss 優點,同時規避掉這些缺點的,是 LUPLUPLUP 分解法。
參考鏈接 LU分解法與Gauss消元法兩者的比較。
所謂承襲性,就是說當計算 Ax=yiAx=y_iAx=yi? (1≤i≤n1 \le i \le n1≤i≤n) 時,若使用 GaussGaussGauss ,那么只能依次調用 nnn 次,總復雜度為 O(k×n3)O(k \times n^3)O(k×n3)。而使用 LULULU 分解法時,先用 O(n3)O(n^3)O(n3) 的復雜度將原式轉化為 LUx=yiLUx=y_iLUx=yi? ,之后進行 kkk 次回代就能求出所有的 xix_ixi? ,總復雜度為 O(n3+k×n2)O(n^3+k \times n^2)O(n3+k×n2)。
另外,在調試代碼的時候,我最初將數組開的很大,試了試 200×200200 \times 200200×200 的矩陣,結果其他方法都能正確求解,唯獨克萊姆法則輸出全為 nannannan 。而對于 100×100100 \times 100100×100 的矩陣,四種方法都能正確求解。
后來發現,使用克萊姆法則,需要求解行列式的值。而我隨機生成的矩陣,這個行列式的值會隨著矩陣的規模增大而急劇增大,溢出 doubledoubledouble。C++黨已經輸麻了
總結
以上是生活随笔為你收集整理的高斯消元法、LU分解法与克莱姆法则解方程组的C++实现的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: w ndows10电脑配置看哪里,Win
- 下一篇: bootstrap-table.js增加