1. 程式人生 > 實用技巧 >淺談高斯消元

淺談高斯消元

作用

高斯消元可以用來解線性方程組,也有一些別的用途。

演算法流程

我們來看這樣的一個方程組:

$ x+y+z=3$

\(2x+3y+z= 6\)

\(4x+y+2z=7\)

然後我們把每一項的係數提取出來,就得到如下的矩形 :

$1 \ 1 \ 1 \ 3 $

$ 2 \ 3 \ 1 \ 6 $

$ 4 \ 1 \ 2 \ 7 $

然後,我們就可以用矩陣的第一行乘\((a[i][1]\div a[1][1])\)加到矩陣的第2行、第3行,使矩陣的2~3行的第1項變為0。如圖:

$ 1 \ \ \ \ 1 \ \ \ \ 1 \ \ \ \ 3 $

$ 0 \ \ \ \ 1 \ -1 \ \ \ \ 0 ​$

$ 0 \ -3 \ -2 \ -5 ​$

然後我們再用第2行乘3並加到第三行,使第3行的第2項邊為0。 如圖:

$ 1 \ \ \ \ 1\ \ \ \ 1 \ \ \ \ 3 $

$ 0 \ \ \ \ 1 \ -1 \ \ \ \ 0 $

$ 0 \ \ \ \ 0 \ -5 \ -5 $

然後我們發現矩陣第3行只有一個未知數\(z\),於是我們就可以直接求出\(z\) 的值,然後第2行就只有一個未知數\(y\)了。求出\(y\)後第一行就只有一個未知數\(x\)。這樣,我們就求出了這個方程組的解。

所以說,我們需要通過某些操作使得矩形變成只有一個倒三角上的數字可以不為0,然後依次求解每一個未知數的值。

注: 對於一個有n個未知數的方程組,若在矩陣處理完畢後前n列有至少一列上的數字全部為0,則該方程有多解,即這一行係數所對應的邊量可以取任意值。

程式碼

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<algorithm>
#define ll long long
#define INF 0x7fffffff
#define re register
#define qwq printf("qwq\n");

using namespace std;

int read()
{
	register int x = 0,f = 1;register char ch;
	ch = getchar();
	while(ch > '9' || ch < '0'){if(ch == '-') f = -f;ch = getchar();}
	while(ch <= '9' && ch >= '0'){x = x * 10 + ch - 48;ch = getchar();}
	return x * f;
}

double a[105][105],Max,x[105];
int n,flag;

void gauss(int n)
{
	int k;
	for(int i = 1; i < n; i++)
	{
		k = i;
		Max = 0;
		for(int j = i + 1; j <= n; j++) if(Max < fabs(a[j][i])) Max = a[j][i], k = j;
		if(k != i) for(int j = i; j <= n + 1; j++) swap(a[k][j],a[i][j]);
		for(int j = i + 1; j <= n; j++)
		{
			double b = (a[j][i] / a[i][i]);
			for(int k = i; k <= n + 1; k++) a[j][k] = a[j][k] - a[i][k] * b;
		}
	}
	for(int i = n; i >= 1; i--)
	{
		if(a[i][i] == 0) {printf("No Solution\n"); flag = 1; return ;}
		for(int j = i + 1; j <= n; j++) a[i][n + 1] = a[i][n + 1] - a[i][j] * x[j];
		x[i] = a[i][n + 1] / a[i][i];
	}
} 

int main()
{
	n = read();
	for(int i = 1; i <= n; i++)
		for(int j = 1; j <= n + 1; j++)
			scanf("%lf",&a[i][j]);
	gauss(n);
	if(flag) return 0;
	for(int i = 1; i <= n; i++) printf("%.2f\n",x[i]);
    return 0;
}