C語言實現高斯消元解線性方程組

目錄

• 前述
• 高斯消元
• 示例
• C代碼
• 最后

前述

在上初中的時候，我們學習了方程組，那個時候只知道消元依次求解。

在大學的時候，我們學習了線性代數，我們可以經過初等變換一次性求出所有解。

這篇文章主要是用C語言實現高斯列主元消去法求解多元一次方程。

高斯列主元消去法

由于涉及到的數學公式太麻煩了，所以從網上找了一張圖片，介紹高斯消去法的，如下圖:

Gauss列主元素法

示例

假如現在有一個三元一次方程組，如下圖：

求解多元一次方程組可以分成三個步驟:

1. 首先根據方程組構建增廣矩陣
2. 其次對增廣矩陣經過行列式的初等變化變成上三角矩陣
3. 最后從后往前回代求解。

構造增廣矩陣

系數矩陣就是將方程組的系數組成矩陣。

而增廣矩陣就是在系數矩陣的右邊添上一列，這一列是線性方程組的等號右邊的值。

下圖即為行列式的增廣矩陣：

組上三角矩陣

這里說的組上三角矩陣是指經過若干步初等變換，將矩陣左上角和右下角連線組成的對角線左下方的元素全部清零。

這個步驟主要涉及到主元以及初等變換兩個概念。

主元指在消去過程中起主導作用的元素，主元通常選擇絕對值最大的元素，用它做除法能夠減小舍入誤差的擴散，使得數值解比較可靠。

以下為行列式的初等變換：

1. 換行變換:交換兩行(列)
2. 倍法變換:將行列式的某一行(列)的所有元素同乘以數k
3. 消法變換:把行列式的某一行(列)的所有元素乘以一個數k并加到另一行(列)的對應元素上

而下面的圖則是經過若干步初等變化組成的上三角矩陣:

迭代求解

在組成上三角矩陣之后，就可以從下往上依次回代求出方程的解了

C代碼

#include
#include

#define?MAX_MATRIX?10

/**
?*?@brief?SwapRow?進行行交換
?*?@param?m?待計算的矩陣
?*????????row?待交行的行
?*????????max_row?待交換的另一行
?*????????n?矩陣行數
?*/
static?void?SwapRow(double?m[][MAX_MATRIX],?int?row,?int?max_row,?int?n)?{
??double?swap;
??for?(int?k?=?row;?k?<=?n;?k++)?{
????swap?=?m[row][k];
????m[row][k]?=?m[max_row][k];
????m[max_row][k]?=?swap;
??}
}

/**
?*?@brief?組上三角矩陣
?*?@param?m?待計算的矩陣
?*????????n?矩陣行數
?*/
static?void?SelectColE(double?m[][MAX_MATRIX],?int?n)?{
??int?max_row_e?=?0;??//主元所在行
??double?ratio?=?0;???//消元因數
??for?(int?j?=?0;?j?????max_row_e?=?j;
????for?(int?i?=?j;?i???????if?(fabs(m[i][j])?>?fabs(m[max_row_e][j]))?{
????????max_row_e?=?i;
??????}
????}
????if?(max_row_e?!=?j)?{
??????SwapRow(m,?j,?max_row_e,?n);???//與最大主元所在行交換
????}
????//消元
????for?(int?i?=?j?+?1;?i???????ratio?=?m[i][j]?/?m[j][j];
??????for?(int?k?=?j;?k?1;?k++)?{
????????m[i][k]?-=?m[j][k]?*?ratio;
??????}
????}
??}
}

/**
?*?@brief:?Gauss?高斯列主元消元法求解線性方程(A*X?=?B)
?*?@param:?m?由于A|B組成的增廣矩陣,X為待求的解
?*?????????n?求解的元數,n要小于MAX_MATRIX
?*?@result:所求結果存放在m[][n]中
?*/
void?Gauss(double?m[][MAX_MATRIX],?int?n)?{
??SelectColE(m,?n);???//?列選主元并消元成上三角
??//?回代求解,結果存在m[][n]中
??for(int?i?=?n?-?1;?i?>=?0;?i--)?{
????for(int?j?=?i?+?1;?j???????m[i][n]?-=?m[i][j]?*?m[j][n];
????}
????m[i][n]?/=?m[i][i];
??}
}

double?a[3][MAX_MATRIX]?=?{
??{3,-1,?1,?4},??//A|B
??{1,?1,?1,?6},
??{2,?3,-1,?12}
};

int?main(int?argc?,char?**argv)?{
??Gauss(a,?3);
??printf("%f,%f,%f\r\n",a[0][3],?a[1][3],?a[2][3]);
??return?0;
}

上述程序運行完成之后，終端輸出：2.000000,3.000000,1.000000

代碼參考鏈接:https://cloud.tencent.com/developer/article/1087352

這里代碼還不太完善，里面沒有處理無解的情況，感興趣的可以繼續完善。

最后

結果很重要，過程同樣重要...