高斯消元1

高斯消元

高斯消元是线性代数规划中的一个算法，可用来为线性方程组求解。所谓线性方程组，是由$M$个$N$元一次方程共同构成的。线性方程组的所有系数可以写成一个$M$行$N$列的“系数矩阵”，再加上每个方程等号右侧的常数，可以写成一个$M$行$N+1$列 的“增广矩阵”
对“增广矩阵”进行的操作有三种

1. 用一个非零的数乘其一行
2. 把其中一行的若干倍加到另一行上
3. 交换两行的位置

称为“初等行变换”，得到“上三角矩阵”或“下三角矩阵”，进一步化简得到“简化阶梯形矩阵”。通过初等行变换把增广矩阵变为简化阶梯形矩阵的线性方程组求解算法就是高斯消元算法。
方程的解写作形如$\{\begin{array}{l}{x}_1=4-2x_2\\ x_3=1\end{array}$的情况下，我们把$x_1,x_3$这样的未知量称为“主元” ，把$x_2$这样的未知量称为“自由元”。
系数和常数构成$n\ast (n+1)$ 的增广矩阵,有三种特殊情况：

1. 最后有$n$个主元，那么有唯一解。
2. 有$k$个主元，自由元$n-k$个，有无穷多个解。
3. 存在系数全为零、常数不为零的行，无解。

以下为两道板子题

P3389 【模板】高斯消元法

题目链接

P3389 【模板】高斯消元法

Code

#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;

int n ;
double a[110][110] ;

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

void gauss()
{
	for ( int i = 1 ; i <= n; ++i )
	{
		for ( int j = i + 1 ; j <= n ; ++j )
			if ( fabs(a[j][i]) > fabs(a[i][i]) )
				for ( int k = 1 ; k <= n + 1 ; ++k )
					swap( a[i][k] , a[j][k] ) ;
		if ( fabs(a[i][i]) < 1e-8 ) continue ;
		for ( int j = 1 ; j <= n ; ++j )
		{
			if ( i == j ) continue ;
			if ( fabs(a[j][i]) < 1e-8 ) continue ;
			double r = a[i][i] / a[j][i] ;
			for ( int k = 1 ; k <= n + 1 ; ++k )
				a[j][k] = a[i][k] - a[j][k] * r ;
		}
	}
}

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() ;
	int ok = 0 ;
	for ( int i = 1 ; i <= n ; ++i )
		if ( fabs(a[i][i]) < 1e-8 && fabs(a[i][n+1]) < 1e-8 )
		{
			ok = 1 ; break ;
		}
	if ( ok ) printf("No Solution\n") ;
	else for ( int i = 1 ; i <= n ; ++i ) printf("%.2lf\n",a[i][n+1] / a[i][i]) ; ;
}

P2455 [SDOI2006]线性方程组

题目链接

P2455 [SDOI2006]线性方程组

题目分析

本题不同于上一题，要判断是无解还是无穷多解，形如

2
0 2 3
0 0 0

一类数据理应为无穷多个解，程序会判为无解，应遵循以下规则进行变换

对于第$i$列，设$j\in [i+1,n]$，如果$b[j][i]>b[i][i]$，直接交换。但是如果相等的话，则按照后面位的系数，将绝对值较小的进行交换。
判断无解的情况要放在判断无限解的情况之前。

Code

#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;

int n ;
double a[105][105] ;

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

bool cmp(int i, int j)
{
	if(fabs(a[j][i]) > fabs(a[i][i])) return 1;
	else if (fabs(a[j][i]) < fabs(a[i][i])) return 0;
	// 如果b[j][i]和b[i][i]绝对值相等,取后面的位数系数绝对值小的
	for (int k = i + 1; k <= n; ++k)
		if(fabs(a[j][k]) < fabs(a[i][k])) return 1;
		else if (fabs(a[j][k]) > fabs(a[i][k])) return 0;
	return 0;
}

void gauss()
{
	for (int i = 1; i <= n; ++i)
	{
		for (int j = i + 1; j <= n; ++j)
			if (fabs(a[j][i]) > fabs(a[i][i]))
				for (int k = 1; k <= n + 1; ++k)
					swap(a[i][k], a[j][k]);
		if (fabs(a[i][i]) < 1e-8) continue;
		for (int j = 1; j <= n; ++j)
		{
			if (i == j) continue;
			if (fabs(a[j][i]) < 1e-8) continue;
			double r = a[i][i] / a[j][i];
			for (int k = 1; k <= n + 1; ++k)
				a[j][k] = a[i][k] - a[j][k] * r;
		}
	}
}

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();
	int ok = 0;
	for (int i = 1; i <= n; ++i)
		if (fabs(a[i][i]) < 1e-8 && fabs(a[i][n+1]) > 1e-8)
		{
			ok = 1; break;
		}
	if ( ok != 1 )
	{
		for (int i = 1; i <= n; ++i)
		{
			if (fabs(a[i][i]) < 1e-8 && fabs(a[i][n+1]) < 1e-8)
			{
				ok = 2; break;
			}
		}
	}
	if (ok == 1) printf("-1\n");
	else if (ok == 2) printf("0\n");
	else for (int i = 1; i <= n; ++i)
	{
		double x = a[i][n+1] / a[i][i];
		if (fabs(x) < 1e-8) x = -x;
		printf("x%d=%.2lf\n",i,x);
	}
}

总结与反思

1. 实型运算会出现-0.00的情况，需要特判
2. 实型运算中极易出现精度丢失，比较关系不能直接使用$==$等运算符，要用绝对值做差与一极小值比较判断大小关系，防止精度丢失
3. 处理细节问题，判断是无解或无穷多个解