實現求解線性方程(矩陣、高斯消去法)------c++程序設計原理與實踐(進階篇)
步驟:
其中A是一個n*n的系數方陣 向量x和b分別是未知數和常量向量:
這個系統可能有0個、1個或者無窮多個解,這取決於系數矩陣A和向量b。求解線性系統的方法有很多,這裏使用一種經典的方法——高斯消去法(https://zh.wikipedia.org/wiki/高斯消去法)。首先,我們對A和b進行交換,使得A變為一個上三角矩陣。上三角矩陣就是對角線之下的所有元素均為0。即如下形式:
實現這個目標是很容易的。為了使a(i,j)變為0,我們先將它乘以一個常量,使它等於第j列上的另一個元素,比如說等於a(k,j)。然後,用第i個方程減去第k個方程,a(i,j)即變為0,矩陣第i行其他元素的值也相應發生改變。
如果這樣一個變換最終使得對角線上所有元素都非0,方程組就有唯一解,此解可以通過”回代“求得。如果存在為0的元素,則意味著方程組有0個或者無窮多個解。
我們現在用c++來表示上述計算方法。首先,定義兩個要使用的具體Matrix類型,以簡化程序:
typedef Numeric_lib::Matrix<double, 2>Matrix2; // Matrix庫下載地址 :http://www.stroustrup.com/Programming/Matrix/Matrix.h 整個庫定義在名字空間 Numeric_lib 中 typedef Numeric_lib::Matrix<double, 1> Vector;
接下來我們將高斯消去法計算過程描述為程序:
Vector classic_gaussian_elimination(Matrix2 A,Vector b){ classical_elimination(A,b); return back_substitution(A,b); }
即,先為兩個輸入A和b創建拷貝(使用賦值函數),然後調用一個函數求解方程組,最後調用回代函數計算結果並將結果返回。關鍵之處在於,我們分解問題的方式和符號表示都完全來自於原始的數學描述。下面所要做的就是實現classic_elimination()和back_substitution()了,解決方案同意完全來自於數學教科書:
void classical_elimination(Matrix2&A,Vector& b){ const Index n=A.dim1(); //從第1列一直遍歷到倒數第二列 //將對角線之下所以元素都填充0 for(Index j=0;j<n-1;++j){ const double pivot =A(j,j); if(pivot==0)cerr<<"錯誤:其中有一對角線位為0"<<‘\n‘; //將第i行在對角線之下的元素都填充為0 for(Index i=j+1;i<n;++i){ const double mult =A(i,j)/pivot; A[i].slice(j)=scale_and_add(A[j].slice(j),-mult,A[i].slice(j)); //A[i].slice(j)表示從A[i][j]到這一行的末尾。 b(i)-=mult*b(j); //對b做對應變化 } } }
“pivot”表示當前行位於對角線上的元素,它必須是非0。因為需要用它作為除數;如果它為0,我們將放棄計算,拋出一個異常:
Vector back_substitution(const Matrix2&A,const Vector&b){ const Index n=A.dim1(); Vector x(n); for(Index i=n-1;i>=0;--i){ double s=b(i)-dot_product(A[i].slice(i+1),x.slice(i+1)); if(double m=A(i,i)) x(i)=s/m; else throw Back_subst_failure(i); } return x; }
完整示例程序:
#include<iostream> #include"Matrix.h" //Matrix庫下載地址 :http://www.stroustrup.com/Programming/Matrix/Matrix.h #include"MatrixIO.h" //MatrixIO庫下載地址 :http://www.stroustrup.com/Programming/Matrix/MatrixIO.h 僅為一維二維提供非常簡單的I/O功能 using namespace Numeric_lib; //整個庫定義在名字空間 Numeric_lib 中 using namespace std; typedef Numeric_lib::Matrix<double, 2>Matrix2; typedef Numeric_lib::Matrix<double, 1> Vector; void classical_elimination(Matrix2& A, Vector& b) { const Index n = A.dim1(); //從第1列一直遍歷到倒數第二列 //將對角線之下所以元素都填充0 for (Index j = 0; j<n - 1; ++j) { const double pivot = A(j, j); if (pivot == 0)cerr<<"錯誤:其中有一對角線位為0"<<‘\n‘; //將第i行在對角線之下的元素都填充為0 for (Index i = j + 1; i<n; ++i) { const double mult = A(i, j) / pivot; A[i].slice(j) = scale_and_add(A[j].slice(j), -mult, A[i].slice(j)); b(i) -= mult*b(j); //對b做對應變化 } } } Vector back_substitution(const Matrix2&A, const Vector&b) { const Index n = A.dim1(); Vector x(n); for (Index i = n - 1; i >= 0; --i) { double s = b(i) - dot_product(A[i].slice(i + 1), x.slice(i + 1)); if (double m = A(i, i)) x(i) = s / m; else cerr<<"錯誤:其中有一對角線位為0"<<‘\n‘; } return x; } Vector classic_gaussian_elimination(Matrix2 A, Vector b) { classical_elimination(A, b); return back_substitution(A, b); } int main() { double val2[3][3] = {2,1,-1,-3,-1,2,-2,1,2 }; double val1[3] = {8,-11,-3 }; Matrix2 A(val2); Vector b(val1); cout<<classic_gaussian_elimination(A, b); }
改進:
pivot為0的問題是可以避免的,我們可以對行進行排列,從而將0和較小的值從對角線上移開,這樣就得到了一個更魯棒的方案。“更魯棒”是指對於舍入誤差不敏感。但是,隨著我們將0置於對角線之下,元素值也會發生改變。因此,我們必須反復進行重排序,以將較小的值從對角線上移開(即,不能一次重排矩陣後就直接使用經典算法):
void elim_with_partial_pivot(Matrix2& A, Vector& b) { const Index n = A.dim1(); for (Index j = 0; j < n; ++j) { Index pivot_row = j; //查找一個合適的主元: for (Index k = j + 1; k < n; ++k) if (abs(A(k, j)) > abs(A(pivot_row, j))) pivot_row = k; //如果我們找到了一個更好的主元,交換兩行: if (pivot_row != j) { A.swap_rows(j, pivot_row); std::swap(b(j), b(pivot_row)); } //消去: for (Index i = j + 1; i < n; ++i) { const double pivot = A(j, j); if (pivot == 0)error("can‘t solve:pivot==0"); const double mult = A(i, j) / pivot; A[i].slice(j) = scale_and_add(A[j].slice(j), -mult, A[i].slice(j)); b(i) -= mult*b(j); } } }
在這裏我們使用了swap_rows()和scale_and_multipy(),這樣程序更符合習慣,我們也不必顯式編寫循環代碼了。
隨機數測試:
void solve_random_system(Index n) { Matrix2 A = random_matrix(n); Vector b = random_vector(n); cout << "A=" << A << ‘\n‘; cout << "b=" << b << ‘\n‘; try { Vector x = classic_gaussian_elimination(A, b); cout << "classic elim solution is x =" << x << ‘\n‘; Vector v = A*x; cout << "A*x=" << v << ‘\n‘; } catch (const exception& e) { cerr << e.what() << ‘\n‘; } }
程序在三種情況下會進入catch語句:
- 代碼中有bug。
- 輸入內容使classic_elimination出現錯誤(elim_with_partial_pivot在很多情況下可以做得更好)。
- 舍入誤差導致問題。
為了測試我們的程序,我們輸出 A*x,其值應該與b相單。但考慮到存在舍入誤差,若其值與b足夠接近就認為結果正確,這也是為什麽測試程序中沒有采用下面語句來判斷結果是否正確的原因:
if(A*b!=b)error ("substitution failed");
在計算機中,浮點數只是實數的近似,因此我們必須接受近似的計算結果。一般來說,應該避免使用==和!=來判斷是否正確。
Matrix庫中並沒有定義矩陣與向量的乘法運算,因此我們為測試程序定義這個運算:
Vector operator*(const Matrix2&m, const Vector&u) { const Index n = m.dim1(); Vector v(n); for (Index i = 0; i < n; ++i) v(i) = dot_product(m[i], u); return v; }
random_matrix()和random_vector()是隨機數的簡單應用。Index是索引類型,它是用typedef定義的。
完整程序:
#include<iostream> #include<random> #include <time.h> #include"Matrix.h" //Matrix庫下載地址 :http://www.stroustrup.com/Programming/Matrix/Matrix.h #include"MatrixIO.h"//MatrixIO庫下載地址 :http://www.stroustrup.com/Programming/Matrix/MatrixIO.h using namespace Numeric_lib;//整個庫定義在名字空間 Numeric_lib 中 using namespace std; typedef Numeric_lib::Matrix<double, 2>Matrix2; typedef Numeric_lib::Matrix<double, 1> Vector; void classical_elimination(Matrix2& A, Vector& b) { const Index n = A.dim1(); //從第1列一直遍歷到倒數第二列 //將對角線之下所以元素都填充0 for (Index j = 0; j<n - 1; ++j) { const double pivot = A(j, j); if (pivot == 0)cerr<<"錯誤:其中有一對角線位為0"<<‘\n‘; //將第i行在對角線之下的元素都填充為0 for (Index i = j + 1; i<n; ++i) { const double mult = A(i, j) / pivot; A[i].slice(j) = scale_and_add(A[j].slice(j), -mult, A[i].slice(j)); b(i) -= mult*b(j); //對b做對應變化 } } } Vector back_substitution(const Matrix2&A, const Vector&b) { const Index n = A.dim1(); Vector x(n); for (Index i = n - 1; i >= 0; --i) { double s = b(i) - dot_product(A[i].slice(i + 1), x.slice(i + 1)); if (double m = A(i, i)) x(i) = s / m; else ; } return x; } void elim_with_partial_pivot(Matrix2& A, Vector& b) { const Index n = A.dim1(); for (Index j = 0; j < n; ++j) { Index pivot_row = j; //查找一個合適的主元: for (Index k = j + 1; k < n; ++k) if (abs(A(k, j)) > abs(A(pivot_row, j))) pivot_row = k; //如果我們找到了一個更好的主元,交換兩行: if (pivot_row != j) { A.swap_rows(j, pivot_row); std::swap(b(j), b(pivot_row)); } //消去: for (Index i = j + 1; i < n; ++i) { const double pivot = A(j, j); if (pivot == 0)error("can‘t solve:pivot==0"); const double mult = A(i, j) / pivot; A[i].slice(j) = scale_and_add(A[j].slice(j), -mult, A[i].slice(j)); b(i) -= mult*b(j); } } } Vector classic_gaussian_elimination(Matrix2 A, Vector b) { elim_with_partial_pivot(A, b); //classical_elimination(A, b); return back_substitution(A, b); } Vector operator*(const Matrix2&m, const Vector&u) { const Index n = m.dim1(); Vector v(n); for (Index i = 0; i < n; ++i) v(i) = dot_product(m[i], u); return v; } int max0 = 10; Vector random_vector(Index n) { Vector v(n); default_random_engine ran{(unsigned int)(time(0)+2)}; uniform_int_distribution<> ureal{ 0,max0 }; for (Index i = 0; i < n; ++i) { v(i) = ureal(ran); } return v; } Matrix2 random_matrix(Index n) { Matrix2 v(n,n); default_random_engine ran{ (unsigned int)time(0) }; uniform_int_distribution<> ureal{ 0,max0 }; for (Index i = 0; i < n; ++i) { for (Index j = 0; j < n; ++j) v[i][j] = ureal(ran); } return v; } void solve_random_system(Index n) { Matrix2 A = random_matrix(n); Vector b = random_vector(n); cout << "A=" << A << ‘\n‘; cout << "b=" << b << ‘\n‘; try { Vector x = classic_gaussian_elimination(A, b); cout << "classic elim solution is x =" << x << ‘\n‘; Vector v = A*x; cout << "A*x=" << v << ‘\n‘; } catch (const exception& e) { cerr << e.what() << ‘\n‘; } } int main() { /*double val2[3][3] = {2,1,-1,-3,-1,2,-2,1,2 }; double val1[3] = {8,-11,-3 }; Matrix2 A(val2); Vector b(val1); cout<<classic_gaussian_elimination(A, b); */ solve_random_system(4); }
c++程序設計原理與實踐(進階篇)
實現求解線性方程(矩陣、高斯消去法)------c++程序設計原理與實踐(進階篇)