天天看點

高斯消元(and牛頓疊代法)詳解及模闆

高斯消元的最主要作用為求解線性方程組,說白了,就是解方程……

解方程還有一種為牛頓疊代法,我們每次枚舉一個值 X0 X 0 ,代入方程看是否為根,不是的話則将 X0 X 0 的值變為:

X0=X0−F(X0)/F'(X0) X 0 = X 0 − F ( X 0 ) / F ′ ( X 0 ) (其中 F'(x) F ′ ( x ) 為 F(x) F ( x ) 的導數)

其證明過程借用了高數中的基本知識泰勒級數,此處不再贅述。

有趣的是,在

雷神之錘Ⅲ

中,卡馬克用牛頓疊代法,令 X0=0x5f3759df X 0 = 0 x 5 f 3759 d f 後求出數字開根後的倒數的速度居然比庫函數快4倍!而這個數字也成為了神奇數字

消元的過程就是模拟人手算解方程的過程

其算法過程就2個

  1. 代入消元
  2. 将算出的值代回,最終求得所有未知數的解

    通俗點,比如說對于這個方程組

    高斯消元(and牛頓疊代法)詳解及模闆
    第一個過程是将其中的一個變量消元化為

    x=num || y=num || z=num

    的形式,其中

    num

    為常數

    第二個過程就是将這個已經求出來的

    x || y || z

    帶回到原方程組中求出剩餘兩個未知數點值;

這裡我們用矩陣模拟這個消元回代的過程

還是用上面的那個方程做例子

高斯消元(and牛頓疊代法)詳解及模闆

轉化為矩陣後為:

高斯消元(and牛頓疊代法)詳解及模闆

其中最後一列為等号右邊的值,設從上到下每行依次為

R1,R2,R3

,從左到右每一列依次為

x,y,z

的系數;

需要注意的是這裡我們要把x最大的放在第一行,這裡是為了在進行除法點時候減小誤差

現在開始化簡:

3*R2-2*R1 ->

高斯消元(and牛頓疊代法)詳解及模闆

3*R3-R1->

高斯消元(and牛頓疊代法)詳解及模闆

R3*(-1)-R2->

高斯消元(and牛頓疊代法)詳解及模闆

到這裡已經可消元完畢了,對應式子如下:

高斯消元(and牛頓疊代法)詳解及模闆

然後就是回代求解了,這裡就不細說了,解方程大家都會

以下模闆來自kuangbin

#include<stdio.h>
#include<algorithm>
#include<iostream>
#include<string.h>
#include<math.h>
using namespace std;
const int MAXN=;
int a[MAXN][MAXN];//增廣矩陣
int x[MAXN];//解集
bool free_x[MAXN];//标記是否是不确定的變元
inline int gcd(int a,int b)
{
    int t;
    while(b!=)
    {
        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=; i<=var; i++)
    {
        x[i]=;
        free_x[i]=true;
    }

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

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

轉載注明出處,謝謝合作