题解 P2455 【[SDOI2006]线性方程组】

update:对一处表述错误进行了修正,求过qwq

常规的思路:高斯消元

高斯消元

这是一种用来求解线性方程组的算法,在方程数较多的时候尤其省时。

对于一个 (n) 元线性方程组:

[egin{cases} a_{1_1}x_1+a_{1_2}x_2+ cdots +a_{1_n}x_n=b_1\ a_{2_1}x_1+a_{2_2}x_2+ cdots +a_{2_n}x_n=b_2\ vdots\ a_{n_1}x_1+a_{n_2}x_2+ cdots +a_{n_n}x_n=b_n\ end{cases} ]

要解它,我们的第一反应显然是消元。

这里给出对于方程组中方程变换最常用的三个法则:

  1. 两个方程位置互换,解不变。

  2. 方程两边同时乘一个非零实数 (k) ,解不变

  3. 方程乘一非零实数 (k) ,再加上另一个方程,解不变

对应到系数矩阵中去,就是矩阵的初等行变换。

我们还知道,对于任意矩阵都能通过初等行变换转化为一个上三角矩阵

所以这个线性方程组的系数矩阵:

[egin{bmatrix} a_{1_1}&a_{1_2}&dots&a_{1_n}\ a_{2_1}&a_{2_2}&dots&a_{2_n}\ vdots&vdots&&vdots&\ a_{n_1}&a_{n_2}&dots&a_{n_n}\ end{bmatrix} ]

可以被转化成上三角形式:

[egin{bmatrix} a_{1_1}&a_{1_2}&dots&a_{1_n}\ &a_{2_2}&dots&a_{2_n}\ &&&vdots\ & & &a_{n_n}\ end{bmatrix} ]

对应到方程就是:

[egin{cases} a_{1_1}x_1+a_{1_2}x_2+ cdots +a_{1_n}x_n=b_1\ a_{2_2}x_2+ cdots +a_{2_n}x_n=b_2\ vdots\ a_{n_n}x_n=b_n\ end{cases} ]

显然最后方程能被直接解出来,解出的根又可以向后代入上一个方程里面。

具体步骤:

枚举每一列c:(相当于把 (x_c) 作为当前处理的列主元)

  1. 找到绝对值最大的一行,将这一行与顶行交换位置,这行不再移动。

这一步是为了将绝对值最大的 (a_{r,c}) 做除数,减少计算时的误差

  1. 等式两边同除,将该行第一个数变为 (1)

系数化一,方便后面求解

  1. 通过等式之间加减,将当前列下面的所有数消为(0)

要构造上三角矩阵,这是必要

  1. 重复至最后一列,然后回代求解。

如果求解过程中发现有方程出现 ( ext{非零常数}=0) 的情况,则方程组无解。

出现了(0=0)的情况,则有无数多组解。

code:

#include <bits/stdc++.h>
using namespace std;

const int N=110;
const double eps=1e-6;//double精度

int n;
double a[N][N];

int gauss()
{
	int c,r;//c表示在枚举哪一列,r表示在枚举哪一行
	for(c=0,r=0;c<n;c++)
	{
		int t=r;
		for(int i=r;i<n;i++)
		{
			if(fabs(a[i][c])>fabs(a[t][c]))//定位到当前列主元的最大系数所在行 
				t=i;
		}
		
		if(fabs(a[t][c])<eps) continue;
		
		for(int i=c;i<=n;i++)	swap(a[t][i],a[r][i]);//交换这两行 
		for(int i=n;i>=c;i--)	a[r][i]/=a[r][c];//同除,系数化一 
		for(int i=r+1;i<n;i++)
		{
			if(fabs(a[i][c])>eps)
			{
				for(int j=n;j>=c;j--)
					a[i][j]-=a[r][j]*a[i][c];//相加减,消去列下面的系数 
			}
		}
		r++;
	}
	
	for(int i=n-1;i>=0;i--)//向后回代求解
	{
		for(int j=i+1;j<n;j++)
		{
			a[i][n]-=a[i][j]*a[j][n];
		}
	}
	
	if(r<n)
	{
		for(int i=r;i<n;i++)
		{
			if(fabs(a[i][n])>eps)
				return 2;//无解
		}
		return 1;//有无穷多组解
	}
	
	return 0;//有唯一解
}

int main()
{
	cin>>n;
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n+1;j++)
		{
			cin>>a[i][j];
		}
	}
	int t=gauss();
	if(t==0)//有唯一解
	{
		for(int i=0;i<n;i++) printf("x%d=%.2lf
",i+1,a[i][n]);
	}
	else if(t==1)puts("0");
	else puts("-1");
	
	return 0; 
}
原文地址:https://www.cnblogs.com/IzayoiMiku/p/13909516.html