高斯消元提供了一种简单的方法解决方程组联立的问题。

基础知识

高斯消元是线性代数中的一类方法,也是一种算法,用来求解线性方程组,并可以求出矩阵的秩,以及求解出可逆方阵的逆矩阵。

根据初等变换原理,将增广矩阵转换为行阶梯矩阵,然后回代求出方程的解。

在程序设计中,常用的算法是:
(我们设方程组中方程的个数为equ,变元的个数为var)

  1. 把方程组转换成增广矩阵。
  2. 利用初等行变换来把增广矩阵转换成行阶梯阵。 枚举k从0到equ – 1,当前处理的列为col(初始为0) ,每次找第k行以下(包括第k行),col列中元素绝对值最大的列与第k行交换。如果col列中的元素全为0,那么则处理col + 1列,k不变。
  3. 转换为行阶梯阵,判断解的情况。

无解
当方程中出现(0, 0, …, 0, a)的形式,且a != 0时,说明是无解的。
唯一解
条件是k = equ,即行阶梯阵形成了严格的上三角阵。利用回代逐一求出解集。
无穷解
条件是k < equ,即不能形成严格的上三角形,自由变元的个数即为equ – k,但有些题目要求判断哪些变元是不缺定的。

这里单独介绍下这种解法:

首先,自由变元有var - k个,即不确定的变元至少有var - k个。我们先把所有的变元视为不确定的。在每个方程中判断不确定变元的个数,如果大于1个,则该方程无法求解。如果只有1个变元,那么该变元即可求出,即为确定变元。

以上介绍的是求解整数线性方程组的求法,复杂度是O(n^3)。浮点数线性方程组的求法类似,但是要在判断是否为0时,加入EPS,以消除精度问题。

贴代码

代码(雍哥修改的不知道谁的)——整型参数型

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
#define maxn 405

int equ,var;
//有equ个方程,val个变元 equ行(0~equ-1) val+1列(增广矩阵)(0~val)
int mat[maxn][maxn];//行列式矩阵
int x[maxn];//解
bool free_x[maxn];//是否是不确定变元

void Debug()
{
for(int i=0;i<equ;++i)
{
for(int j=0;j<=var;++j)
printf("%01d ",mat[i][j]);
puts("");
}puts("");
}
inline int gcd(int a,int b){return b==0?a:gcd(b,a%b);}
inline int lcm(int a,int b){return a*b/gcd(a,b);}

//-2 有浮点数解 -1 无解 0有一个解 >0无穷解,返回自由变元个数
int gauss()
{
int k,max_r,col;//当前行,当前列绝对值最大的行;当前列
//初始化
memset(x,0,sizeof(x));
memset(free_x,true,sizeof(free_x));
//转换为阶梯矩阵
for(k=0,col=0;k<equ&&col<var;++k,++col)
{
max_r=k;
for(int i=k+1;i<equ;++i)
if(abs(mat[i][col])>abs(mat[max_r][col]))max_r=i;
if(max_r!=k)
for(int i=k;i<=var;++i)swap(mat[k][i],mat[max_r][i]);
if(!mat[k][col]){--k;continue;}
for(int i=k+1;i<equ;++i)if(mat[i][col])
{
int l=lcm(abs(mat[i][col]),abs(mat[k][col]));
int ta=l/abs(mat[i][col]),tb=l/abs(mat[k][col]);
if(mat[i][col]*mat[k][col]<0)tb=-tb;
for(int j=col;j<=var;++j)
mat[i][j]=mat[i][j]*ta-mat[k][j]*tb;
}
}//Debug();
//1.无解的情况 存在行(0,0,...,a) a!=0
for(int i=k;i<equ;++i)
if(mat[i][col])return -1;
//2.无穷解的情况 秩不是equ
if(k<var)
{
for(int i=k-1;i>=0;--i)
{
int free_x_num=0,free_index;
for(int j=0;j<var;++j)
if(mat[i][j]&&free_x[j])++free_x_num,free_index=j;
if(free_x_num>1)continue;
int tmp=mat[i][var];
for(int j=0;j<var;++j)
if(mat[i][j]&&j!=free_index)tmp-=mat[i][j]*x[j];
x[free_index]=tmp/mat[i][free_index];
free_x[free_index]=false;
}
return var-k;
}
//3.唯一解 严格三角阵 满秩
for(int i=var-1;i>=0;--i)
{
int tmp=mat[i][var];
for(int j=i+1;j<var;++j)
if(mat[i][j])tmp-=mat[i][j]*x[j];
if(tmp%mat[i][i])return -2;//有浮点数
x[i]=tmp/mat[i][i];
}
return 0;
}

对于浮点参数型,基本实现形式差不多,只不过求解的时候注意点浮点数的运算就行了。

代码(来自这里):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
#define MAXN 405
const double EPS=1e-9;

double A[MAXN][MAXN];
double X[MAXN];
bool Free[MAXN];
int equation, variable;

void debug()
{
printf("The Matrix is :\n");
for(int i=0; i<=variable+1; i++)
{
printf("%8d",i);
if(i==0) printf("|");
}
printf("\n");
for(int i=0; i<=variable+1; i++)
{
printf("-------");
if(i==0) printf("-+");
else printf("-");
}
printf("\n");

for(int i=1; i<=equation; i++)
{
printf("%8d|",i);
for(int j=1; j<=variable+1; j++)
{
printf("%8.2lf",A[i][j]);
}
printf("\n");
}
printf("\n");
}

/***************************************************************************/
/// 功能说明:
/// 1、不需带入参数,参数已经在 build() 函数中设置好了
/// 2、函数返回 -1 :无解
/// 3、函数返回 0 :有且仅有一组解
/// 4、函数返回 其他正数 :有多解,返回的值是不确定变元的个数
/**************************************************************************/
int gauss()
{
int col=1, row=1;
for(; row<=equation && col<=variable; row++, col++)
{
int max_row=row;

//找到 col 那列元素绝对值最大的那行与当前行交换,减小精度误差
for(int i=row+1; i<=equation; i++)
if(abs(A[i][col])>abs(A[max_row][col]))
max_row=i;

//如果 col 那列元素最大是 0,表明这一列全部是 0,处理下一列
if(abs(A[max_row][col])<=EPS)
{
row--;
continue;
}

//如果不是同一行,交换元素
if(max_row!=row)
{
for(int i=col; i<=variable+1; i++)
swap(A[max_row][i], A[row][i]);
}

//枚举要删去的行
//for(int i=row+1; i<=equation; i++)
for(int i=1; i<=equation; i++)
{
if(i==row) continue;
if(abs(A[i][col])<=EPS) continue;

double ta=A[row][col];
double tb=A[i][col];

for(int j=1; j<=variable+1; j++)
A[i][j]=A[i][j]*ta-A[row][j]*tb;
}
// printf("row=%d col=%d\n", row, col);
// debug();
}

// debug();

// 1、没有解的情况:
for(int i=row; i<=equation; i++)
if(A[i][variable+1]!=0) return -1; //表明无解

// 2、无穷解的情况:
if(row<variable+1)
{
for(int i=row-1; i>=1; i--)
{
int free_x_num=0, free_index=-1;
for(int j=1; j<=variable; j++)
if(Free[j] && A[i][j]!=0)
{
free_x_num++;
free_index=j;
}
if(free_x_num>1) continue;

Free[free_index]=false;

double temp=A[i][variable+1];
for(int j=1; j<=variable; j++)
if(!(abs(A[i][j])<=EPS) && j!=free_index)
temp-=A[i][j]*X[j];
X[free_index]=temp/A[i][free_index];
}
return variable-row+1;
}

// 3、唯一解的情况:
for(int i=variable; i>=1; i--)
{
double temp=A[i][variable+1];
for(int j=i+1; j<=variable; j++)
if(!(abs(A[i][j])<=EPS))
temp-=A[i][j]*X[j];
X[i]=temp/A[i][i];
}
return 0;
}

void build(int n, int m)
{
memset(Free, true, sizeof(Free));
memset(X, 0, sizeof(X));
variable=n, equation=m;
for(int i=1; i<=equation; i++)
for(int j=1; j<=variable+1; j++)
scanf("%lf",&A[i][j]);
debug();
}

void PRIMAXNT(int flag)
{
if(flag==-1)
{
printf("MAXNO Solution !\n");
return;
}
if(flag==0)
{
for(int i=1; i<=variable; i++)
printf("x%d : %4.4lf\n", i, X[i]);
return;
}
printf("There are %d variables who are indeterminate !\n",flag);
for(int i=1; i<=variable; i++)
{
printf("x%d : ", i);
if(Free[i]) printf("indeterminate !\n");
else printf("%4.4lf\n",X[i]);
}
}

主要解决的问题

高斯消元主要解决这几类问题:

普通的消元问题

给出几个方程组,求出每个解。一般套用模板就可以了。

染色问题

一个方格只有01两种状态,改变其中一个方格会影响其它周围方格的状态。

对于这种问题,首先要建立方程组。关于方程组的建立,首先要把方格转换成一行,然后构造一个 n * n 的矩阵。那么对于矩阵中的( i,j )就表示j可以影响i的状态。

那么如果 i 能影响 j 点,就应该把( j,i )设为1。

用高斯消元的时候,注意把中间过程 mod 2,因为只有 01 两种状态。

mod p问题

方程组对p取余,相当于一个线性同余型的问题吧。

对于这种问题,分为两种情况:

  1. 如果p为 质数 ,这样就可以参考染色问题,根据性质,直接将中间过程 mod p,最后得出的解保证小于 p 即可。
  2. 如果p为 非质数 ,这样就不能对中间过程 mod p,因为这样会形成循环节而扩大解集。我也没有什么好的解决方法,暂时可以中间过程不要进行任何取余操作。求解集的时候可以枚举 0~p 的数字,满足同余关系就可以作为解的一部分。这个也要注意解不要超过 p 。

找出逆矩阵

高斯消元法可以用来找出一个可逆矩阵的逆矩阵。设 A 为一个 N * N 的矩阵,其逆矩阵可被两个分块矩阵表示出来。将一个 N * N 单位矩阵放在 A 的右手边,形成一个 N * 2N 的分块矩阵 B = [ A,I ] 。经过高斯消元法的计算程序后,矩阵 B 的左手边会变成一个单位矩阵 I ,而逆矩阵 A - 1 会出现在B 的右手边。

假如高斯消元法不能将A 化为三角形的格式,那就代表A 是一个不可逆的矩阵。

应用上,高斯消元法极少被用来求出逆矩阵。高斯消元法通常只为线性方程组求解。

参考文献