高斯消元法的C++簡單實現

高斯消元法

　　首先，我們導入幾個概念。

定義1： 一個矩陣稱為階梯形（行階梯形），若它有以下三個性質：

　　1.每一非零行在每一零行之上；

　　2.某一行的先導元素所在的列位於前一行先導元素的后面；

　　3.某一行先導元素所在列下方元素都是零。

　　比如，

定義2：若一個階梯形矩陣還滿足以下性質，稱它為簡化階梯形（簡化行階梯形）：

　　1.每一非零行的先導元素是1；

　　2.每一先導元素1是該元素所在列的惟一非零元素。

　　比如，

定理1：每個矩陣行等價於惟一的簡化階梯形矩陣。即簡化階梯形矩陣是惟一的。

---

 　　下面，我們用一個具體例子來說明高斯消元法的主要步驟。

　　原矩陣：

　　第一步，由最左的非零列開始，這是一個主元列。主元位置在該列頂端。

　　第二步，在主元列中選取一個非零元作為主元。若有必要的話，對換兩行使這個元素移到主元位置上。

　　第三步，用倍加行變換將主元下面的元素變成0.

　　第四步，暫時不管包含主元位置的行以及它上面的各行，對剩下的子矩陣使用上述的三個步驟直到沒有非零行需要處理為止。

　　對每一行重復上述步驟。

　　第五步，由最右面的主元開始，把每個主元上方的各元素變成0.若某個主元不是1，用倍乘變換將它變成1.

　　最后，我們就得到了原矩陣的簡化階梯形。

　　其中，第1~4步稱為行化簡算法的向前步驟，產生唯一的簡化階梯形的第5步，稱為向后步驟。

 

---

 

C++實現

　　我們嘗試用C++來實現以上步驟。這里只是簡單的實現，也就是用代碼描述了上述步驟，沒有考慮過多的問題。歡迎大家在評論里指出問題，提出更好的建議，以便於日后改進。

　　大概的實現思路就是先實現向前步驟：

　　首先，我們對於每一行找到第一個不為零的元素，並且將這一行置為1 * * * *的形式，用這一行乘上倍數加到之后的每一行。

　　再實現向后步驟：

　　然后，我們從最后一行開始，選擇主元，加到之前的每一行上，使得該列的元素都為零。

　　最后，我們就完成了化簡，得到了簡化階梯形。

　　以上算法只是一個粗略實現，主要體現在：

　　1.對於主元的選定不夠最優；

　　2.會出現精度問題；

　　3.對於某些情況無法處理。

　　先暫時貼上代碼，之后有時間再進行優化。

#include <iostream>
#include <cstdio>

using namespace std;

int main()
{
    double martix[100][100];
    int n, m; // n行m列

    scanf("%d %d", &n, &m);

    // 輸入
    for(int i = 0; i < n; i++)
        for(int j = 0; j < m; j++)
            scanf("%lf", &martix[i][j]);

    // 向前步驟
    for(int i = 0; i < n - 1; i++)
    {
        // 找主元
        int pos = 0;
        for(int j = 0; j < m; j++)
            if(martix[i][j])
            {
                pos = j;
                break;
            }

        if(martix[i][pos] != 1 && martix[i][pos] != 0)
        {
            double tmp = martix[i][pos];
            for(int j = pos; j < m; j++)
            {
                martix[i][j] = martix[i][j] / tmp;
            }
        }
        for(int j = i + 1; j < n; j++)
        {
            if(!martix[j][pos])
                continue;
            double tmp = martix[j][pos];
            for(int k = pos; k < m; k++)
            {
                martix[j][k] = martix[j][k] - martix[i][k] * tmp;
            }
        }
    }

    // 向后步驟
    for(int i = n - 1; i > 0; i--)
    {
        int pos = 0;
        for(int j = 0; j < m; j++)
            if(martix[i][j])
            {
                pos = j;
                break;
            }

        if(martix[i][pos] != 1 && martix[i][pos] != 0)
        {
            double tmp = martix[i][pos];
            for(int j = pos; j < m; j++)
            {
                martix[i][j] = martix[i][j] / tmp;
            }
        }

        for(int j = 0; j < i; j++)
        {
            if(!martix[j][pos])
                continue;
            double tmp = martix[j][pos];
            for(int k = pos; k < m; k++)
            {
                martix[j][k] = martix[j][k] - martix[i][k] * tmp;
            }
        }
    }

    // 輸出
    for(int i = 0; i < n; i++)
    {
        for(int j = 0; j < m; j++)
            printf("%-10.2f", martix[i][j]);
        printf("\n");
    }
    return 0;
}