高斯消元求多元一次方程(講解+板子)


一、我們解多元一次方程需要什么?

因為未知數有多個,所以我們需要方程的數量也不同

如果你要求n元一次方程,那么你至少需要給出n個方程才可能會求出來所有未知量的大小

 

二、增廣矩陣是什么?

我們求多元一次方程的方法就是按照大學線性代數課程中的方法。

我們首先需要構造一個增廣矩陣,然后我們把這個增廣矩陣化為上三角矩陣,如果順利的話就可以求出來所有未知數的解,但也有可能只能求出來一些未知數的解。無解情況也會處理。

 

1、增廣矩陣:就是在系數矩陣的右邊添上一列,這一列是線性方程組的等號右邊的值。(下面給出n個未知量,m個方程構成的增廣矩陣)

a11x1+a12x2+a13x3+...+a1nxn=b1

a21x1+a22x2+a23x3+...+a2nxn=b2

......

 

am1x1+am2x2+am3x3+...+amnxn=bm

 

 

 

2、上三角矩陣:由增廣矩陣經過一系列初等變化而成。(最終形式如下,以4個未知數,4個方程為例)

a11   a12     a13    a14

     a22    a23   a24

0      0        a33    a34

     0        0        a44

 

就是第一行第一個元素下面的值都要為0、第二行第二個元素下面的值都要為0、第三行第三個元素下面的值都要為0......

 

簡單說一下初等變化中比較實用的兩個規則(使用初等變化之后矩陣所有性質都是不改變的):

1、矩陣中任意兩行或者兩列可以交換

2、矩陣中一行中每一個元素,可以加上另一行對應的位置k倍的元素,也就是

a11+k*a21    a12+k*a22   ...    a1n+k*a2n  可以替換掉下面這一行(k可以取任意值)

a11               a12                      a1n

 

三、矩陣的秩

就是將矩陣轉化為上三角矩陣之后不是全0行的行數。(全0行,也就是矩陣這一行所有的值都是0)

示例(R就代表矩陣的秩):

 

 

 

 

 

 

 上圖截取自百度文庫:點我跳轉

 

 

四、什么是列主元消去法


列主元消去法是在系數矩陣中按列選取元素絕對值得最大值作為主元素(可以通過交換行的方式來滿足這個需求)

這樣做的原因就是可以減小誤差,因為求解方程的時候不能避免出現相除的情況,而計算機內部不能保存分數,也就導致了結果和預想的不一樣。

例如,你讓變量a=1/3,b=3.你輸出a*b的結果就不是1

 

代碼:

    col=0; // 當前處理的列
    for(k = 0;k < equ && col < var;k++,col++)
    {// 枚舉當前處理的行.
     // 找到該col列元素絕對值最大的那行與第k行交換.(為了在除法時減小誤差)
        max_r=k;
        for(i=k+1;i<equ;i++)
        {
            if(abs(a[i][col])>abs(a[max_r][col])) max_r=i;
        }
        if(max_r!=k)
        {// 與第k行交換.
            for(j=k;j<var+1;j++) swap(a[k][j],a[max_r][j]);
        }
    }

 

 

五、怎么把系數矩陣轉換成上三角矩陣

首先,按照列主元消去法把這一列中最大元素那一行和第一行交換。之后讓其他行的第一個元素變成0,那么就會滿足公式

a11*k+ay1=0  或者   a11+k*ay1=0     (y>1)

k=-ay1/a11     或者   k=-a11/ay1

這樣的話就讓第y行每一個元素加上a11*k  或者   ay1*k   (這就看你選擇那個公式了,其實都差不多)

 

代碼:

for(i=k+1;i<equ;i++)
        {// 枚舉要刪去的行.
            if(a[i][col]!=0)
            {
                LCM = lcm(abs(a[i][col]),abs(a[k][col]));
                ta = LCM/abs(a[i][col]);
                tb = LCM/abs(a[k][col]);
                if(a[i][col]*a[k][col]<0)tb=-tb;//異號的情況是相加
                for(j=col;j<var+1;j++)
                {
                    a[i][j] = a[i][j]*ta-a[k][j]*tb;
                }
            }
        }

 

六、判斷矩陣的解

如果你是n元一次方程,而且你給了n個方程,且這n個方程化成的上三角矩陣每一行都至少有一個元素不為0(要排除0 0 ... 0 0 a,a!=0,因為着一種情況就代表着每一個系數都是0,但是結果卻不是0,這樣的情況是無解的)

那么這樣的方程n個未知數都可以求出來

 

如果你有n個未知數,但是你只給出來m各方程(m<n),那么這個注定是求不出來所有未知數的解的(可能會求出來一部分未知數)

如果有k個未知數沒有求出來,那就以為着有k-1個變元,就比如x+y=1,如果y確定了,那么x也就確定了。所以變元就只有一個(從這里我們也看出來了,如果存在變元,那么這個方程的解也就有無數個)

 

代碼:

// 高斯消元法解方程組(Gauss-Jordan elimination).(-2表示有浮點數解,但無整數解,
//-1表示無解,0表示唯一解,大於0表示無窮解,並返回自由變元的個數)
//有equ個方程,var個變元。增廣矩陣行數為equ,分別為0到equ-1,列數為var+1,分別為0到var.

// 1. 無解的情況: 化簡的增廣陣中存在(0, 0, ..., a)這樣的行(a != 0).
    for (i = k; i < equ; i++)
    {
        // 對於無窮解來說,如果要判斷哪些是自由變元,那么初等行變換中的交換就會影響,則要記錄交換.
        if (a[i][col] != 0) return -1;
    }
    // 2. 無窮解的情況: 在var * (var + 1)的增廣陣中出現(0, 0, ..., 0)這樣的行,即說明沒有形成嚴格的上三角陣.
    // 且出現的行數即為自由變元的個數.
    if (k < var)
    {
        // 首先,自由變元有var - k個,即不確定的變元至少有var - k個.
        for (i = k - 1; i >= 0; i--)
        {
            // 第i行一定不會是(0, 0, ..., 0)的情況,因為這樣的行是在第k行到第equ行.
            // 同樣,第i行一定不會是(0, 0, ..., a), a != 0的情況,這樣的無解的.
            free_x_num = 0; // 用於判斷該行中的不確定的變元的個數,如果超過1個,則無法求解,它們仍然為不確定的變元.
            for (j = 0; j < var; j++)
            {
                if (a[i][j] != 0 && free_x[j]) free_x_num++, free_index = j;
            }
            if (free_x_num > 1) continue; // 無法求解出確定的變元.
            // 說明就只有一個不確定的變元free_index,那么可以求解出該變元,且該變元是確定的.
            temp = a[i][var];
            for (j = 0; j < var; j++)
            {
                if (a[i][j] != 0 && j != free_index) temp -= a[i][j] * x[j];
            }
            x[free_index] = temp / a[i][free_index]; // 求出該變元.
            free_x[free_index] = 0; // 該變元是確定的.
        }
        return var - k; // 自由變元有var - k個.
    }
    // 3. 唯一解的情況: 在var * (var + 1)的增廣陣中形成嚴格的上三角陣.
    // 計算出Xn-1, Xn-2 ... X0.
    for (i = var - 1; i >= 0; i--)
    {
        temp = a[i][var];
        for (j = i + 1; j < var; j++)
        {
            if (a[i][j] != 0) temp -= a[i][j] * x[j];
        }
        if (temp % a[i][i] != 0) return -2; // 說明有浮點數解,但無整數解.
        x[i] = temp / a[i][i];
    }
    return 0;

 

七、全部代碼

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

const int MAXN=50;



int a[MAXN][MAXN];//增廣矩陣
int x[MAXN];//解集
bool free_x[MAXN];//標記是否是不確定的變元



/*
void Debug(void)
{
    int i, j;
    for (i = 0; i < equ; i++)
    {
        for (j = 0; j < var + 1; j++)
        {
            cout << a[i][j] << " ";
        }
        cout << endl;
    }
    cout << endl;
}
*/

inline int gcd(int a,int b)
{
    int t;
    while(b!=0)
    {
        t=b;
        b=a%b;
        a=t;
    }
    return a;
}
inline int lcm(int a,int b)
{
    return a/gcd(a,b)*b;//先除后乘防溢出
}

// 高斯消元法解方程組(Gauss-Jordan elimination).(-2表示有浮點數解,但無整數解,
//-1表示無解,0表示唯一解,大於0表示無窮解,並返回自由變元的個數)
//有equ個方程,var個變元。增廣矩陣行數為equ,分別為0到equ-1,列數為var+1,分別為0到var.
int Gauss(int equ,int var)
{
    int i,j,k;
    int max_r;// 當前這列絕對值最大的行.
    int col;//當前處理的列
    int ta,tb;
    int LCM;
    int temp;
    int free_x_num;
    int free_index;

    for(int i=0; i<=var; i++)
    {
        x[i]=0;
        free_x[i]=true;
    }

    //轉換為階梯陣.
    col=0; // 當前處理的列
    for(k = 0; k < equ && col < var; k++,col++)
    {
        // 枚舉當前處理的行.
// 找到該col列元素絕對值最大的那行與第k行交換.(為了在除法時減小誤差)
        max_r=k;
        for(i=k+1; i<equ; i++)
        {
            if(abs(a[i][col])>abs(a[max_r][col])) max_r=i;
        }
        if(max_r!=k)
        {
            // 與第k行交換.
            for(j=k; j<var+1; j++) swap(a[k][j],a[max_r][j]);
        }
        if(a[k][col]==0)
        {
            // 說明該col列第k行以下全是0了,則處理當前行的下一列.
            k--;
            continue;
        }
        for(i=k+1; i<equ; i++)
        {
            // 枚舉要刪去的行.
            if(a[i][col]!=0)
            {
                LCM = lcm(abs(a[i][col]),abs(a[k][col]));
                ta = LCM/abs(a[i][col]);
                tb = LCM/abs(a[k][col]);
                if(a[i][col]*a[k][col]<0)tb=-tb;//異號的情況是相加
                for(j=col; j<var+1; j++)
                {
                    a[i][j] = a[i][j]*ta-a[k][j]*tb;
                }
            }
        }
    }

    //  Debug();

    // 1. 無解的情況: 化簡的增廣陣中存在(0, 0, ..., a)這樣的行(a != 0).
    for (i = k; i < equ; i++)
    {
        // 對於無窮解來說,如果要判斷哪些是自由變元,那么初等行變換中的交換就會影響,則要記錄交換.
        if (a[i][col] != 0) return -1;
    }
    // 2. 無窮解的情況: 在var * (var + 1)的增廣陣中出現(0, 0, ..., 0)這樣的行,即說明沒有形成嚴格的上三角陣.
    // 且出現的行數即為自由變元的個數.
    if (k < var)
    {
        // 首先,自由變元有var - k個,即不確定的變元至少有var - k個.
        for (i = k - 1; i >= 0; i--)
        {
            // 第i行一定不會是(0, 0, ..., 0)的情況,因為這樣的行是在第k行到第equ行.
            // 同樣,第i行一定不會是(0, 0, ..., a), a != 0的情況,這樣的無解的.
            free_x_num = 0; // 用於判斷該行中的不確定的變元的個數,如果超過1個,則無法求解,它們仍然為不確定的變元.
            for (j = 0; j < var; j++)
            {
                if (a[i][j] != 0 && free_x[j]) free_x_num++, free_index = j;
            }
            if (free_x_num > 1) continue; // 無法求解出確定的變元.
            // 說明就只有一個不確定的變元free_index,那么可以求解出該變元,且該變元是確定的.
            temp = a[i][var];
            for (j = 0; j < var; j++)
            {
                if (a[i][j] != 0 && j != free_index) temp -= a[i][j] * x[j];
            }
            x[free_index] = temp / a[i][free_index]; // 求出該變元.
            free_x[free_index] = 0; // 該變元是確定的.
        }
        return var - k; // 自由變元有var - k個.
    }
    // 3. 唯一解的情況: 在var * (var + 1)的增廣陣中形成嚴格的上三角陣.
    // 計算出Xn-1, Xn-2 ... X0.
    for (i = var - 1; i >= 0; i--)
    {
        temp = a[i][var];
        for (j = i + 1; j < var; j++)
        {
            if (a[i][j] != 0) temp -= a[i][j] * x[j];
        }
        if (temp % a[i][i] != 0) return -2; // 說明有浮點數解,但無整數解.
        x[i] = temp / a[i][i];
    }
    return 0;
}
int main(void)
{
//    freopen("in.txt", "r", stdin);
//    freopen("out.txt","w",stdout);
    int i, j;
    int equ,var;
    while (scanf("%d %d", &equ, &var) != EOF)
    {
        memset(a, 0, sizeof(a));
        for (i = 0; i < equ; i++)
        {
            for (j = 0; j < var + 1; j++)
            {
                scanf("%d", &a[i][j]);
            }
        }
//        Debug();
        int free_num = Gauss(equ,var);
        if (free_num == -1) printf("無解!\n");
        else if (free_num == -2) printf("有浮點數解,無整數解!\n");
        else if (free_num > 0)
        {
            printf("無窮多解! 自由變元個數為%d\n", free_num);
            for (i = 0; i < var; i++)
            {
                if (free_x[i]) printf("x%d 是不確定的\n", i + 1);
                else printf("x%d: %d\n", i + 1, x[i]);
            }
        }
        else
        {
            for (i = 0; i < var; i++)
            {
                printf("x%d: %d\n", i + 1, x[i]);
            }
        }
        printf("\n");
    }
    return 0;
}

 

 

八、樣例輸入

示例:求方程x+y=1和x-y=1的解

2方程 2未知數

系數矩陣:

1 1 1

1 -1 1

 

 

 

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM