高斯消元的最主要作用為求解線性方程組,說白了,就是解方程……
解方程還有一種為牛頓疊代法,我們每次枚舉一個值 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個
- 代入消元
-
将算出的值代回,最終求得所有未知數的解
通俗點,比如說對于這個方程組
第一個過程是将其中的一個變量消元化為高斯消元(and牛頓疊代法)詳解及模闆
的形式,其中x=num || y=num || z=num
num
為常數
第二個過程就是将這個已經求出來的
帶回到原方程組中求出剩餘兩個未知數點值;x || y || z
這裡我們用矩陣模拟這個消元回代的過程
還是用上面的那個方程做例子

轉化為矩陣後為:
其中最後一列為等号右邊的值,設從上到下每行依次為
R1,R2,R3
,從左到右每一列依次為
x,y,z
的系數;
需要注意的是這裡我們要把x最大的放在第一行,這裡是為了在進行除法點時候減小誤差
現在開始化簡:
3*R2-2*R1 ->
3*R3-R1->
R3*(-1)-R2->
到這裡已經可消元完畢了,對應式子如下:
然後就是回代求解了,這裡就不細說了,解方程大家都會
以下模闆來自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 ;
}
轉載注明出處,謝謝合作