線性方程組的直接解法——Gauss消去法

WeShiko睡大覺發表於2022-12-29

考慮線性方程組

\[\mathrm{A}x=\mathrm{b} \]

其中，\(\mathrm{A}=(a_{ij})_{n\times n}\)，\(\mathrm{b}=[b_1,b_2,\cdots,b_n]^{\mathrm{T}}\)。線上性代數的課程中，我們已經學習過Gauss消元法，具體操作是將矩陣A轉化為“階梯型”矩陣。為方便起見，本文僅僅討論係數矩陣非奇異的方程組，此時，目標是將矩陣A轉化為上三角矩陣，再執行回代過程，即可給出方程組的解。本文將給出在計算機上的具體操作及例項程式碼。

一、基本Gauss消去法

我們僅僅討論對矩陣第一列的操作，剩餘的操作可以以此類推，因而不再贅述。
在執行Gauss消去法時，我們將第一列對角元以下的元素全部變為零。記第一列消元操作後的增廣矩陣為\([\mathrm{A}^{(1)},\mathrm{b}^{(1)}]\)，容易知道

\[[\mathrm{A}^{(1)},\mathrm{b}^{(1)}]= \begin{bmatrix} a_{11} & a_{22} & \cdots &a_{1n} & b_1 \\ 0 & a_{22}^{(1)} &\cdots &a_{2n}^{(1)} & b_2^{(1)}\\ \vdots &\vdots & & \vdots &\vdots\\ 0 & a_{n2}^{(1)} & &a_{nn}^{(1)} & b_n^{(1)} \end{bmatrix}\]

其中

\[a_{ij}^{(1)}=a_{ij}-\frac{a_{i1}}{a_{11}}a_{1j}，j=2,\cdots ,n \]

\[a_{i1}^{(1)}=0 \]

\[b_i^{(1)}=b_i-\frac{a_{i1}}{a_{11}}b_1 \]

觀察到重複出現的結構\(\frac{a_{_{i1}}}{a_{_{11}}}\)，我們記它為\(l_{i1}\)，稱為消元因子，並將它儲存在原來\(a_{i1}\)的位置。在計算的過程中，先計算消元因子並儲存在相應位置，再執行後續的演算法。
對於後續部分的運算，在第k步，只要對矩陣\(A^{(k-1)}(k:n,k:n)\)執行相同操作即可。

二、列主元Gauss消去法

在執行Gauss消元法的過程中，如果\(a_{kk}^{(k-1)}\)相對於其他元素絕對值較小，則會產生較大的舍入誤差，影響計算精度，為此，我們引入了列主元Gauss消去法，基於交換矩陣的行不影響線性方程組的解。
記執行完k-1步消元后的增廣矩陣為\([\mathrm{A}^{(k-1)},\mathrm{b}^{(k-1)}]\)。考慮第k列對角元及其以下的部分。選擇絕對值最大的元所在行，與當前行執行行交換，再進行Gauss消元法。

三、計算例項

用列主元Gauss消去法解以下線性方程組：

\[\left\{ \begin{array}{} 0.5x_1+1.1x_2+3.1x_3=6,\\ 2x_1+4.5x_2+3.6x_3=0.02,\\ 5x_1+0.96x_2+6.5x_3=0.96. \end{array} \right.\]

#include <iostream>
#include <math.h>
using namespace std;

int main()
{
    double A_Extended[3][4]={0.5,1.1,3.1,6,2,4.5,3.6,0.02,5,0.96,6.5,0.96};
    double X_solution[3];
    for (int i=0;i<=2;i++)
    {
        int n=i;
        for (int p=i+1;p<=2;p++)
        {
            if (fabs(A_Extended[p][i])>fabs(A_Extended[n][i]))
            {
                n=p;
            }
        }

        for (int p=i;p<=2+1;p++)
        {
            double k=A_Extended[n][p];
            A_Extended[n][p]=A_Extended[i][p];
            A_Extended[i][p]=k;
        }

        for (int p=i+1;p<=2;p++)
        {
            A_Extended[p][i]=-A_Extended[p][i]/A_Extended[i][i];
            for (int pco=i+1;pco<=2+1;pco++)
            {
                A_Extended[p][pco]=A_Extended[p][pco]+A_Extended[p][i]*A_Extended[i][pco];
            }
        }
    }
    X_solution[2]=A_Extended[2][3]/A_Extended[2][2];
    for (int i=1;i>=0;i--)
    {
        double sum=0;
        for (int k=2;k>i;k--)
        {
            sum=sum+A_Extended[i][k]*X_solution[k];
        }
        X_solution[i]=(A_Extended[i][3]-sum)/A_Extended[i][i];
    }

    cout<<X_solution[0]<<" "<<X_solution[1]<<" "<<X_solution[2]<<endl;
    return 0; 
}

相關文章