高斯消元
*原题:题目背景
输入一个包含 n 个方程 n 个未知数的线性方程组。
方程组中的系数为实数。
求解这个方程组。
下图为一个包含 m 个方程 n 个未知数的线性方程组示例:
输入格式
第一行包含整数 n。
接下来 n 行,每行包含 n+1 个实数,表示一个方程的 n 个系数以及等号右侧的常数。
输出格式
如果给定线性方程组存在唯一解,则输出共 n 行,其中第 i 行输出第 i 个未知数的解,结果保留两位小数。
如果给定线性方程组存在无数解,则输出 Infinite group solutions。
如果给定线性方程组无解,则输出 No solution。
数据范围
1≤n≤100,
所有输入系数以及常数均保留两位小数,绝对值均不超过 100。
输入样例:
3
1.00 2.00 -1.00 -6.00
2.00 1.00 -3.00 -9.00
-1.00 -1.00 2.00 7.00
输出样例:
1.00
-2.00
3.00
解题思路:利用高斯消元法,将系数矩阵转换为阶梯矩阵,时间复杂度为O(n3)。阶梯矩阵就是每行的第一个非0元素为1,或者该行全部为0(系数矩阵该行全部为0)。为方便计算,需要将答案作为单独一列加入到系数矩阵
转换为阶梯行列式的步骤
1-> 为方便计算,找到某一列绝对值最大的元素的那一个行,将该行作为操作行,将其置顶,如果操作行数为i,那么需要将该行第一个非0元素转换为1,并将该元素下面的所有元素清0,通过方程式变形和两两之间的运算来修改未知数的系数(详情可参考线性代数矩阵变为阶梯矩阵)
以下为操作代码
int c,r; //表示行和列
//化解行列式,将行列式化为最简阶梯行列式
for(c=0,r=0;c<n;c++) //遍历每一列,将该列绝对值最大的数所在行置顶,方便计算
{
int t=r; //记录最大绝对值元素所在行
for(int i=r;i<n;i++) //每次操作还未确立的行,找到此列绝对值最大的元素所在行,将其记录
if(fabs(g[i][c])>fabs(g[t][c])) t=i; //寻找最大的数值是因为可以避免系数变得太大,精度较高.
if(fabs(g[t][c])<eps) continue; //绝对值最大等于(接近)0,那么所有该列数都是0,无需进一步操作
for(int i=c;i<n+1;i++) //否则从第c列开始交换行,从第c列开始是操作完成的数据就无需再操作了
swap(g[t][i],g[r][i]); //r为最顶行,t为最大绝对值所在行,按列交换
for(int i=n;i>=c;i--) //将最大绝对值元素所在行转换为最简行列式(阶梯行列式),第一个非0元素,变为1,其余元素除以相同的值
g[r][i]/=g[r][c]; //到过来算是因为没有存储使其变为1的约数,顺序算会导致本身变为1,则后面的元素都会除以1
for(int i=r+1;i<n;i++) //在把该列该行下面的所有元素消为0
if(fabs(g[r][c])>eps) //该元素不为0,则操作,否则无需操作
for(int j=n;j>=c;j--) //方程两边进行相同的操作
g[i][j]-=g[r][j]*g[i][c]; //每次操作是减掉此行第一个非0元素与上一行第一个非0元素的积,这样就能保住此行第一个非0
//元素变为0
r++; //操作完一轮之后,r++,原来的最上面一行就操作完成了
}
当我们推出 r<n时,矩阵则有无解和无穷多解
因为r<n代表有一行系数全为0,那么等式左边为0,如果等式右边为0,那么不管未知数怎么取,等式都成立,即有无穷多解。若右边不等于0,则等式矛盾,不管未知数怎么取值,等式都不成立,即无解
//所有行和列转换完成后,判断是否有解
if(r<n) //说明剩下的方程个数是小于n的,解不为一,判断是无解还是无穷多解
{
for(int i=r;i<n;i++) //r及一下的面的行,由于已经转换为阶梯行列式,所以所有系数都为0,即等式左边全为0,若等式右边
{ //不为0,那么左边不等于右边,则无解,否则左边等于右边,有无穷多组解
if(fabs(g[i][n])>eps)
return 0; //无解
else
return 2; //无穷多解
}
}
如果r=n,代表n个未知数有n个一次方程,那么有唯一解
关于如何求解每个未知数,我们只需消掉后面的未知数即可,如果是第1行,则消掉第二个到最后一个元素,第二行则消掉第3个到最后一个元素,依次类推,且从倒数第二行向上推,最后一个无需推理,在先前计算阶梯矩阵时,答案已经被计算出来。
每一行为了消去应该消去的未知数,建立循环,每次消去下一行和下一个元素
else //r==n 有唯一解,从下往上推出每个未知数
{
for(int i=n-1;i>=0;i--) //从下往上遍历每一行,当前这一行减去下一行与当前这一行的相应系数,便于消除一个未知数
for(int j=i+1;j<n+1;j++)
g[i][n]-=g[j][n]*g[i][j]; //等式左右两边进行相同操作
//如两个方程 x1+3*x2=d1 x2=d2,那么 将1式减去2式*3则可以得出x1答案,因为是阶梯行列式
//所以需要依次向上推,最后一行无需推,最后一行由于先前计算阶梯矩阵已经计算出结果了
return 1; //唯一解
}
每次为了消去相应元素,如第一行消元,第一次循环应该是消去第2个元素,即行数为i,减j=i+1行第j=i+1个元素,直到循环到j=n,而i行的答案矩阵的元素减去下一行的答案矩阵的元素和当前行第j个元素的成绩(如x1+3x2=d1,x2=d2,那么x1=d2-d13)(自行根据代码推导更容易理解)
完整代码
#include <iostream>
#include <cmath>
#include <algorithm>
using namespace std;
const int N=200;
double g[N][N]; //系数矩阵,从00开始,g[0][0]表示第一行x1的系数
int n;
const double eps = 1e-6; //误差不超过eps,当a<eps时,可以将a认为为0,不直接写等于0是防止计算机在运算时精度丢失
int gauss() //返回0表示无解,返回1表示唯一解,返回2表示无穷多解
{
int c,r; //表示行和列
//化解行列式,将行列式化为最简阶梯行列式
for(c=0,r=0;c<n;c++) //遍历每一列,将该列绝对值最大的数所在行置顶,方便计算
{
int t=r; //记录最大绝对值元素所在行
for(int i=r;i<n;i++) //每次操作还未确立的行,找到此列绝对值最大的元素所在行,将其记录
if(fabs(g[i][c])>fabs(g[t][c])) t=i; //寻找最大的数值是因为可以避免系数变得太大,精度较高.
if(fabs(g[t][c])<eps) continue; //绝对值最大等于(接近)0,那么所有该列数都是0,无需进一步操作
for(int i=c;i<n+1;i++) //否则从第c列开始交换行,从第c列开始是操作完成的数据就无需再操作了
swap(g[t][i],g[r][i]); //r为最顶行,t为最大绝对值所在行,按列交换
for(int i=n;i>=c;i--) //将最大绝对值元素所在行转换为最简行列式(阶梯行列式),第一个非0元素,变为1,其余元素除以相同的值
g[r][i]/=g[r][c]; //到过来算是因为没有存储使其变为1的约数,顺序算会导致本身变为1,则后面的元素都会除以1
for(int i=r+1;i<n;i++) //在把该列该行下面的所有元素消为0
if(fabs(g[r][c])>eps) //该元素不为0,则操作,否则无需操作
for(int j=n;j>=c;j--) //方程两边进行相同的操作
g[i][j]-=g[r][j]*g[i][c]; //每次操作是减掉此行第一个非0元素与上一行第一个非0元素的积,这样就能保住此行第一个非0
//元素变为0
r++; //操作完一轮之后,r++,原来的最上面一行就操作完成了
}
//所有行和列转换完成后,判断是否有解
if(r<n) //说明剩下的方程个数是小于n的,解不为一,判断是无解还是无穷多解
{
for(int i=r;i<n;i++) //r及一下的面的行,由于已经转换为阶梯行列式,所以所有系数都为0,即等式左边全为0,若等式右边
{ //不为0,那么左边不等于右边,则无解,否则左边等于右边,有无穷多组解
if(fabs(g[i][n])>eps)
return 0; //无解
else
return 2; //无穷多解
}
}
else //r==n 有唯一解,从下往上推出每个未知数
{
for(int i=n-1;i>=0;i--) //从下往上遍历每一行,当前这一行减去下一行与当前这一行的相应系数,便于消除一个未知数
for(int j=i+1;j<n+1;j++)
g[i][n]-=g[j][n]*g[i][j]; //等式左右两边进行相同操作
//如两个方程 x1+3*x2=d1 x2=d2,那么 将1式减去2式*3则可以得出x1答案,因为是阶梯行列式
//所以需要依次向上推,最后一行无需推,最后一行由于先前计算阶梯矩阵已经计算出结果了
return 1; //唯一解
}
}
int main(int argc, char** argv) {
cin>>n;
for(int i=0;i<n;i++)
for(int j=0;j<n+1;j++)
scanf("%lf",&g[i][j]);
int t=gauss();
if(t==1)
{
for(int i=0;i<n;i++)
{
printf("%.2lf\n",g[i][n]);
}
}
else if(t==2)
puts("Infinite group solutions");
else puts("No solution");
return 0;
}