【數值計算】數值解析--n元一次聯立方程組:直接解法
加減法(中學所學)是我們平常用的解法之一。 例如,現有如下所示的二元一次方程組。
將等式兩邊同乘以一個實數,上下係數合併,消去其中一元未知數的方法便是熟知的加減法。
之後,把帶入式1,解得。
把上式用行列式表示如下,
之後,第2行乘以,上下相減,得到
的形式。隨後,從最下面的式子中解出未知數,代入,得到:
前面的操作是指,把係數行列的左下部分(不包括對角元素)全部變成0,求解如下所示的三角行列(upper triangle matrix)。
這樣的處理稱作前進消去(forward elimination)。 根據前進消去,未知數可通過解出, 把
高斯消去法的程式碼實現
首先,把常數向量新增到的最後一列組成的行列。
其中第i行j列的元素表示為
由於要在C,C++上實現,這裡的元素序號從0開始計數。
前進消去可以表示成下式
這裡的
根據上式,把用更新得到上三角行列。這一步的C++程式碼如下。
// 前進消去(forward elimination) // - 將左下角元素變為0 for(int k = 0; k < n-1; ++k){ double akk = A[k][k]; for(int i = k+1; i < n; ++i){ double aik = A[i][k]; for(int j = k+1; j < n+1; ++j){ A[i][j] = A[i][j]-aik*(A[k][j]/akk); } } }
這裡的二維陣列A[n][n+1]中儲存著係數行列和常數向量。
後退代入的公式如下。
C++程式碼:
// 後代入(back substitution) // - x_n的解為b_n/a_nn,x_n,代入n-1行求出x_(n-1) A[n-1][n] = A[n-1][n]/A[n-1][n-1]; for(int i = n-2; i >= 0; --i){ double ax = 0.0; for(int j = i+1; j < n; ++j){ ax += A[i][j]*A[j][n]; } A[i][n] = (A[i][n]-ax)/A[i][i]; }
計算結果儲存於A[i][n]中。
高斯消去法的完整程式碼如下。
/*!
* 高斯消去法(無主元交換)
* @param[inout] A 為n×n的係數項與n×1的常數項(b)合併後n×(n+1)的行列
* @param[in] n n元一次方程組
* @return 1:成功
*/
int GaussElimination(vector< vector<double> > &A, int n)
{
// 前進消去(forward elimination)
// - 左下角元素為0
for(int k = 0; k < n-1; ++k){
double akk = A[k][k];
for(int i = k+1; i < n; ++i){
double aik = A[i][k];
for(int j = k+1; j < n+1; ++j){
A[i][j] = A[i][j]-aik*(A[k][j]/akk);
}
}
}
// 後退代入(back substitution)
A[n-1][n] = A[n-1][n]/A[n-1][n-1];
for(int i = n-2; i >= 0; --i){
double ax = 0.0;
for(int j = i+1; j < n; ++j){
ax += A[i][j]*A[j][n];
}
A[i][n] = (A[i][n]-ax)/A[i][i];
}
return 1;
}
上面使用STL的vector代替2維陣列。
主元選擇
我們先重新看一遍高斯消去法中的前進消去公式。
會發現右側的第二項中包含除以的計算。當前進消去過程中對角元素中有0時,則會停止計算。即使不為0,當有絕對值很小的數出現時也會有誤差出現。
針對上述問題,這裡我們對主元消元法(pivotal elimination method)進行說明。首先,我們知道聯立方程的公式的順序是可以交換的。也就是說,處理第行時,把其對角元素與同列還未處理的元素 的值進行比較,把絕對值最大的那一行與行互換後,應用到前進消去公式中。上述操作便是主元消元。 像這樣的只交換行的操作稱為部分主元消元(partial pivoting), 列方向也交換的操作稱為完全主元消元(complete pivoting)。值得注意的是,部分主元消元中未知數的順序不發生變化,完全主元消元中,未知數的順序會發生變化。
包含主元消元的高斯消去法程式碼如下
/*!
* 主元選擇(Pivoting)
* - 只交換行的部分主元消去法
* @param[inout] A
* @param[in] n n元一次方程組
* @param[in] k 物件行
* @return 1:成功
*/
inline void Pivoting(vector< vector<double> > &A, int n, int k)
{
// 檢索k行後k列的絕對值最大元素所在行
int p = k; // 絕對值最大的行
double am = fabs(A[k][k]);// 最大値
for(int i = k+1; i < n; ++i){
if(fabs(A[i][k]) > am){
p = i;
am = fabs(A[i][k]);
}
}
// 若k != p,行交換(主元選擇)
if(k != p) swap(A[k], A[p]);
}
/*!
* 高斯消去法(pivoting)
* @param[inout] A
* @param[in] n n元一次方程組
* @return 1:成功
*/
int GaussEliminationWithPivoting(vector< vector<double> > &A, int n)
{
// 前進消去(forward elimination)
for(int k = 0; k < n-1; ++k){
// 主元消去
Pivoting(A, n, k);
double akk = A[k][k];
for(int i = k+1; i < n; ++i){
double aik = A[i][k];
for(int j = k; j < n+1; ++j){
A[i][j] = A[i][j]-aik*(A[k][j]/akk);
}
}
}
// 後退代入(back substitution)
A[n-1][n] = A[n-1][n]/A[n-1][n-1];
for(int i = n-2; i >= 0; --i){
double ax = 0.0;
for(int j = i+1; j < n; ++j){
ax += A[i][j]*A[j][n];
}
A[i][n] = (A[i][n]-ax)/A[i][i];
}
return 1;
}
高斯・約當法
將高斯消去法的前進消去擴增,可以求出正方行列的逆行列,這種方法稱作高斯約當法(Gauss-Jordan method)。首先我們假設是一個的正方行列。 若 ,則是一個正則行列(regular matrix)(反之為奇異行列(singular matrix))。如果把的單位行列表示成,且非奇異,則只存在一個行列滿足。 這個可以用逆矩陣表示。