如果给定一个由equ个方程,var个变元组成的方程组,然后问解集的话,就要用到高斯消元法了。我们把系数放入一个矩阵,在其最后一列放原来每个方程等号右边的值,形成了一个equ*(var+1)的矩阵,每一次去枚举行(当然列是先确定的,从左向右),找到本行下面当前列里绝对值最大的那个数,找到后,将这两行交换(为了在除法时减小误差)。
接下来,要把矩阵转成阶梯矩阵,比如
X1 X2 X3 X4 X5
0 X2 X3 X4 X5
0 0 X3 X4 X5
0 0 0 X4 X5
0 0 0 0 X5
从上到下,变元个数依次减少,方法就是每次用当前行的方程去对下面行的方程进行消元处理
形成之后,判断解。
1.如果有 0 0 0 ......0 0 T 这种行出现,且T!=0 那么就是无解
2.如果有 0 0 0 ......0 0 0 出现,就是无穷多个解。
3.形成严格的上三角矩阵,逐次回代,就能得到唯一的解。
详参模板代码:
接下来,要把矩阵转成阶梯矩阵,比如
X1 X2 X3 X4 X5
0 X2 X3 X4 X5
0 0 X3 X4 X5
0 0 0 X4 X5
0 0 0 0 X5
从上到下,变元个数依次减少,方法就是每次用当前行的方程去对下面行的方程进行消元处理
形成之后,判断解。
1.如果有 0 0 0 ......0 0 T 这种行出现,且T!=0 那么就是无解
2.如果有 0 0 0 ......0 0 0 出现,就是无穷多个解。
3.形成严格的上三角矩阵,逐次回代,就能得到唯一的解。
详参模板代码:
inline int gcd(int a,int b)//消元的时候为计算扩大的倍数而要用到。
{
int t;
while(b!=0)
{
t=b;
b=a%b;
a=t;
}
return a;
}
inline int lcm(int a,int b)
{
return a/gcd(a,b)*b;
}
int Gauss(int equ,int var)//equ个方程,对应equ行,var个变元,对应var行
{
int i,j,k;
int max_r;//当前列里绝对值最大的那行
int col;//当前处理列,从0开始
int ta,tb;
int LCM;
int temp,free_index,free_num;
for (int i = 0; i <= var; i++)
{
x[i] = 0;
free_x[i] = 1;//初始化他们都不是自由变元
}
col=0;//从第0列开始处理
for(k=0;k<equ && col<var;k++,col++)
{
max_r=k;
for(i=k+1;i<equ;i++)//找到第k+1行到第equ-1行里col列上绝对值最大的数
{
if(abs(mat[i][col]) > abs(mat[max_r][col]))
max_r=i;
}
if(max_r!=k)
{
for(j=k;j<var+1;j++)//交换
swap(mat[max_r][j],mat[k][j]);
}
if(mat[k][col]==0)//如果第k行的第col列为0,那么k+1行开始col列一定全为0,此时就去处理下一列
{
k--;
continue;
}
for(i=k+1;i<equ;i++)//消元
{
if(mat[i][col]!=0)
{
LCM=lcm(abs(mat[k][col]),abs(mat[i][col]));
ta=LCM/abs(mat[i][col]);//倍数
tb=LCM/abs(mat[k][col]);//倍数
if(mat[k][col] * mat[i][col]<0)//异号的话是相加
tb=-tb;
for(j=col;j<var+1;j++)
mat[i][j]=mat[i][j]*ta-mat[k][j]*tb;//把k行下面的每个方程都消掉一个元
}
}
}
for(i=k;i<equ;i++)//系数全是0的行一定只在k-equ-1,如果有增广列上不为0,则无解
if(mat[i][col] != 0)
return -1;
if(k<var)//如果没有形成严格的上三角矩阵的话,有无穷多个解,而自由元个数有var-k个,那么不确定的变元至少有var-k个
{
for(i=k-1;i>=0;i--)
{
free_num=0;//不确定变元个数
for(j=0;j<var;j++)//寻找每一行上的不确定变元个数
{
if(mat[i][j]!=0 && free_x[j])
{
free_num++;
free_index=j;//记录这个不确定变元的列序号
}
}
if(free_num>1)//超过1个,无法求出
continue;
temp=mat[i][var];
for(j=0;j<var;j++)//只有一个的话,他是确定的
{
if(mat[i][j]!=0 && j!=free_index)
temp-=x[j]*mat[i][j];
}
x[free_index]=temp/mat[i][free_index];
free_x[free_index]=0;//标记他是确定的变元
}
return var-k;//自由元个数
}
for(i=var-1;i>=0;i--)//如果是有唯一解,那么逐次回代即可求出解
{
temp=mat[i][var];
for(j=i+1;j<var;j++)
if(mat[i][j]!=0)
temp-=mat[i][j]*x[j];
if(temp%mat[i][i]!=0)//说明没有整数解,但是有浮点数解
return -2;
x[i]=temp/mat[i][i];
}
return 0;
}