第3 章线性方程组直接求解
万里长征,是一步一步走出来的。线性方程与线性方程组的直接求解
3.1 引 言
工程技术、科学研究中的很多实际问题,如电学中的网络问题、最小二
乘法的曲线拟合问题、三次样条插值问题、自由振动问题等,都需要求解线
性方程组。n 阶线性方程组的一般形式为
a11x1 +a12x2 + … +a1nxn =b1 
a21x1 +a22x2 + … +a2nxn =b2 
. . . . . 
an1x1 +an2x2 + … +annxn =bn 
ì
.
í
... 
..
. 
式中,aij和bi 为常数;xi 为变元,这里i,j=1,2,…,n。
n 阶线性方程组一般形式的矩阵形式为
Ax =b 
式中, 
A = 
a11 a12 … a1n 
a21 a22 … a2n 
. . . . 
an1 an2 … ann 
é
.
êêêêê 
ù
.
úúúúú
, x = 
x1 
x2 
. 
xn 
é
.
êêêêê 
ù
.
úúúúú
, b = 
b1 
b2. 
bn 
é
.
êêêêê 
ù
.
úúúúú 
其中,A 为系数矩阵,b 为右端向量,x 为解向量。常把右端向量b 并入系
数矩阵A,作为A 的第n+1列,记为
..A = [A b] = 
a11 a12 … a1n b1 
a21 a22 … a2n b2 
. . . . . 
an1 an2 … ann bn 
é
.
êêêêê 
ù
.
úúúúú 
称..A 为增广矩阵。
线性方程组分类的一种方法是按阶数高低分类,阶数在150以上的线

0 
　计算方法(Python 
版)

性方程组为高阶线性方程组;阶数在100 以下的线性方程组为低阶线性方程组。线性方
程组分类的另一种方法是按系数矩阵
A 
中零元素的多少来分类,系数矩阵
A 
的大部分元
素(大约在80% 以上)是零元素的线性方程组为稀疏线性方程组;系数矩阵
A 
的大部分元
素(大约在80% 以上)是非零元素的线性方程组为稠密线性方程组。实际应用中,经常见
到高阶稀疏线性方程组、低阶稠密线性方程组、对称正定线性方程组和带状线性方程
组等。

由克拉默法则可知,若系数行列式det(A)≠0,则线性方程组Ax=
b 
有唯一解。但
是克拉默法则并不实用。当用克拉默法则求解
n 
阶线性方程组时,总乘除次数MD= 
(n+1)个行列式×[(1)n!次乘法+
n 
次除法]。随着方程组阶数的增加,求解的

n2-· 
总计算量以不低于指数级别的速度增长。例如,当阶数n=20 时,总乘除次数MD≈ 
9.7×104×1035×10

20;当阶数n=30 时,MD≈7.36;当阶数n=40 时,MD≈5.52,这时就

算用每秒做1千万亿次乘除法的计算机求解,耗费的年数也是个天文数字。

线性方程组常用的数值解法主要分为两类:直接求解方法和迭代求解方法。

直接求解方法是指经过有限次四则运算,可以求出线性方程组精确解的方法。虽然
运算方法在理论上没有误差,但实际计算过程一般会产生舍入误差。有时舍入误差大得
无法忽略,因此在设计算法时必须考虑舍入误差的增长,估计运算结果误差的大小。本章
将介绍直接求解方法中的消元法和三角分解法,这些方法求解
n 
阶线性方程组的时间复
杂度为O(空间复杂度为O(求解低阶稠密矩阵时效率较高。

n3), n2), 

迭代求解方法是指构造一种迭代方法,由某个(套)迭代初值(粗略解), 得到近似解序
列,用序列极限逐步逼近线性方程组精确解的方法。与第2章非线性方程迭代求根相比
较,线性方程组迭代求解需要一套数据才能够启动迭代过程,一次迭代的结果也包含多个
数值,却仍为单点迭代法。迭代求解方法存在收敛性和收敛速度的问题。如果收敛,那么
随着运算量的增大,误差会趋向无穷小。迭代求解的程序设计简单,内存需求小。求解高
阶稀疏矩阵时,迭代求解方法效率较高。第4章将介绍线性方程组的迭代求解方法。

3.顺序高斯消元法
2 

用顺序高斯消元法求解线性方程组的过程分为消元过程和回代过程。消元过程是自
上而下,把原线性方程组化为上三角方程组的过程;回代过程是自下而上,对这个上三角
方程组进行求解的过程。

3.2.1 
消元过程
消元的过程是对线性方程组做同解变换,把原线性方程组变为上三角方程组的过程。
消元的次序,是先用第1个方程消去它之下所有方程左边第1个变元,并更新其他变元的
系数,完成第1次循环;然后不再用到第1个方程(最上边1行)和它下边所有方程的最左
边1列(变元系数全为0), 对右下剩余的矩形区域,重新消去最上边方程之下所有方程的
最左边1列,更新其他变元的系数,完成第2次循环;如此反复,直到把这个线性方程组化
为上三角方程组。具体地说,消去某一个方程的第一个变元,就是让矩形区域最上边方程


第3 章　线性方程组直接求解　 51 
等号的左右两端都乘以1个因子,加到待消元的方程上,使待消元方程第一个变元的系数
为0,从而实现这一步消元。
为了用程序实现消元和回代,需要用符号表示方程组的一般形式。n 阶线性方程组
在消元过程一开始,可以表示为
a1(11)x1 +a1(12)x2 + … +a1(1,n)xn =a1(1,n)+1 
a2(11)x1 +a2(12)x2 + … +a2(1,n)xn =a2(1,n)+1 
. . 
a(1) n,1x1 +a(1) n,2x2 + … +a(1) n,nxn =a(1) n,n+1 
ì
.
í
... 
... 
下面说明其中符号的含义: 
(1)xj 为变元。下标j 为变元序号,且j 为列号。
(2)a(k) ij 为变元xj 的系数。下标i 为行号,下标j 为列号。a(k) ij 每更新(改变)一次, 
上标k 增1。k 初值为1,因此k 为a(k) ij 的更新次数+1。
外层循环每循环一轮,就用涉及的矩形区域中上边的第1行消去它下边各行左边的
第1列,此矩形区域最上1行和最左1列在之后的消元过程中不再涉及(上标不再改变), 
除此之外的右下矩形区域中所有的系数被更新一遍(上标k 增加1)。第k 次消元的初始
状态为
a1(11)x1 +a1(12)x2 + … +a1(1,k)xk + … +a1(1,n)xn =a1(1,n)+1 
a2(22)x2 + … +a2(2,k)xk + … +a2(2,n)xn =a2(2,n)+1 
. . . . 
a(k) k,kxk + … +a(k) k,nxn =a(k) k,n+1 
. . . 
a(k) n,kxk + … +a(k) n,nxn =a(k) n,n+1 
ì
.
í
.... 
.... 
在第k 次更新的初始,涉及的待更新矩形区域中,系数的上标都为k,此矩形区域中
最左上角的系数的行标、列标也为k。更新的过程,就是用第k 行消去第k+1行至第n 
行的第k 列,并且让第k+1行至第n 行、第k+1列至第n+1列的矩形区域中的系数更
新,上标变为k+1。在此后的消元过程中,第k 行及之上和第k 列及其左边的系数不再
变化。若
以上标k 为循环变量,则每循环1次,消去1列。消元过程中,k 共循环n-1次, 
消去n-1列,把原方程组变为上三角方程组。上标k 每循环1次,需要更新第k+1行
至第n 行、第k+1列至第n+1列的矩形区域中的系数,这一过程需要2重循环,因此消
元过程总共需要3重循环,外层循环变量为上标k,中层循环变量为行标i,内层循环变量
为列标j。
消元结束之后得到的上三角方程组中,每向下一行,上标k 增加1,所以上标k 等于
行标i,如下所示: 
a1(11)x1 +a1(12)x2 + … +a1(1,n)xn =a1(1,n)+1 
a2(22)x2 + … +a2(2,n)xn =a2(2,n)+1 
. . . 
a(n) n,nxn =a(n) n,n+1 
ì
.
í
... 
...

52 　计算方法(Python 版) 
下面举例说明消元的具体过程,如表3.1所示。
表3.1 消元的具体过程
实 例符号表示
x+2y-z=2 ① 
3x-y+z=4 ② 
3x+2y-2z=1 ③ { a1(11)x1+a1(12)x2+a1(13)x3=a1(14) ① 
a2(11)x1+a2(12)x2+a2(13)x3=a2(14) ② 
a3(11)x1+a3(12)x2+a3(13)x3=a3(14) ③ { 
①× -31
( ) +②→② 
①× -31
( ) +③→③ 
令l21=-a2(11) 
a1(11),则①×l21+②→② 
令l31=-a3(11) 
a1(11),则①×l31+③→③ 
x+2y-z=2 ① 
-7y+4z=-2 ② 
-4y+z= -5 ③ { a1(11)x1+a1(12)x2+a1(13)x3=a1(14) ① 
a2(22)x2+a2(23)x3=a2(24) ② 
a3(22)x2+a3(23)x3=a3(24) ③ { 
②× --4 
( -7) +③→③ 令l32=-a3(22) 
a2(22),则②×l32+③→③ 
x+2y-z=2 ① 
-7y+4z=-2 ② 
-97
z=-27 
7 ③ 
ì
.
í
.. 
.. 
a1(11)x1+a1(12)x2+a1(13)x3=a1(14) ① 
a2(22)x2+a2(23)x3=a2(24) ② 
a3(33)x3=a3(34) ③ { 
表3.1中的li,k 称为行乘子,其计算公式为
li,k =-
a(k) i,k 
a(k) k,k i ∈ [k +1,n] 
当外层循环为第k 次循环时,用方程k 消去方程i 的变元xk ,实现对第i 行的系数
a(k) i,j 的1次更新的过程为:将方程k(第k 行)的对应系数a(k) k,j(j∈[k,n+1])乘以行乘子
li,k ,加上a(k) i,j ,得到a(k+1) i,j 。这一过程的计算公式为
a(k+1) i,j =a(k) k,j ×li,k +a(k) i,j i ∈ [k +1,n],j ∈ [k,n +1] (3.1) 
式(3.1)的主要目的是消去方程i 的变元xk ,即要求a(k+1) i,k =0,那么a(k+1) i,k =a(k) k,k × 
li,k +a(k) i,k =0,对此等式变形可得li,k =-a(k) i,k 
a(k) k,k ,与行乘子li,k 的计算公式一致。
当a(k) k,k =0时,上述计算过程无法进行;当a(k) k,k 的绝对值很小时,误差会被放大。
定义3.1 a(k) k,k(k=1,2,…,n)对求解起突出作用,称a(k) k,k 为主元素(主元)。
线性方程组的消元过程可以用矩阵来表示: 
..A = [A b] 
= 
a1(11) a1(12) … a1(1,n) a1(1,n)+1 
a2(11) a2(12) … a2(1,n) a2(1,n)+1 
. . . . . 
a(1) n,1 a(1) n,2 … a(1) n,n a(1) n,n+1 
é
.
êêêêê
ê 
ù
.
úúúúú
ú 
一系列倍加变换→ 
a1(11) a1(12) … a1(1,n) a1(1,n)+1 
a2(22) … a2(2,n) a2(2,n)+1 
. . . 
a(n) n,n a(n) n,n+1 
é
.
êêêêê
ê 
ù
.
úúúúú
ú

第3 章　线性方程组直接求解　 53 
消元过程对应于增广矩阵..A = [A b] 的一系列倍加变换,最终把系数矩阵A 化为
上三角矩阵。用方程k 消去方程i 的变元xk 的过程,对应于增广矩阵..A 的1次倍加变
换,即用行乘子li,k 乘以..A 的第k 行,加到第i 行上去。
3.2.2 回代过程
在消元过程中把原线性方程组变形为上三角方程组,然后在回代过程中,自下而上地
求解这个上三角方程组。当回代到某一行时,把已求出的变元代入,就变成求解只有1个
未知变元的线性方程。
下面举例说明回代的具体过程,如表3.2所示。
表3.2 回代的具体过程
实 例符号表示
x+2y-z=2 ① 
-7y+4z=-2 ② 
-97
z=-27 
7 ③ 
ì
.
í
.. 
.. 
a1(11)x1+a1(12)x2+a1(13)x3=a1(14) ① 
a2(22)x2+a2(23)x3=a2(24) ② 
a3(33)x3=a3(34) ③ { 
由③得
z= 
-27 
7 
-97
=3 
由③得
x3=a3(34) 
a3(33) 
把z=3代入②得
y=-2-4z 
-7 =2 
把x3 代入②得
x2=a2(24)-a2(23)x3 
a2(22) 
把z=3,y=2代入①得
x=2-(2y-z)=1 
把x3 和x2 代入①得
x1 =
a1(14)- Σ3 
j=2(a1(1,j)xj) 
a1(11) 
若行标i、列标j 和上标k 都从1开始,则对n 阶线性方程组回代的一般公式为
xk =
a(k) k,(n+1)- Σn 
j=k+1(a(k) k,jxj) 
a(k) k,k , k =n,n -1,n -2,…,1 
可以依次求出xn ,xn-1,xn-2,…,x1。
3.2.3 顺序高斯消元法的算法和程序
顺序高斯消元法的几点说明如下。
(1)当用第k 行消去第i 行的变元xk 时,a[i][k]变为0,并且不参与之后的运算。
为了节省存储空间,用a[i][k]存放行乘子li,k 。
(2)与上面的叙述不同,算法和程序中列表元素的下标从0开始。
(3)程序没有对线性方程组无解和主元为零时求不出解的情况进行判断。

54 　计算方法(Python 版) 
算法3.1 顺序高斯消元法的算法。
输入原方程组的阶数n 
输入原方程组的增广矩阵a[n][n+1] 
forkinrange(n-1),循环1次消去1列 (消元过程) 
foriinrange(k+1,n),循环1次更新1行
行乘子a[i][k]/=-a[k][k] 
forjinrange(k+1,n+1) 
a[i][j]+=a[i][k]*a[k][j] 
forkinrange(n-1,-1,-1) (回代过程) 
s=0 
forjinrange(k+1,n) 
s+=a[k][j]*x[j] 
x[k]=(a[k][n]-s)/a[k][k] 
输出原方程组的解向量x 
程序3.1 顺序高斯消元法对应的程序。 
def inputdata(a,n): 
print("请输入原方程组的增广矩阵: ") 
for i in range(n): 
for j in range(n+1): 
a[i][j]=eval(input('a[{}][{}]='.format(i,j))) 
def outputdata(x,n): 
print("原方程组的解为: ") 
for i in range(n): 
print('x[{}]={}'.format(i,x[i]),end=' ') 
n=int(input("请输入原方程组的阶数: ")) 
a=[[0]*(n+1) for i in range(n)] 
inputdata(a,n) 
for k in range(n-1): 
for i in range(k+1,n): 
a[i][k]/=-a[k][k] 
for j in range(k+1,n+1): 
a[i][j]+=a[i][k]*a[k][j] 
x=[0]*n 
for k in range(n-1,-1,-1): 
s=0 
for j in range(k+1,n): 
s+=a[k][j]*x[j]

第3 章　线性方程组直接求解　 55 
x[k]=(a[k][n]-s)/a[k][k] 
outputdata(x,n) 
上述程序和算法的时间复杂度为O(n3),空间复杂度为O(n2)。具体地说,顺序高斯
消元法的乘除次数MD=n3 
3+n2-n3
≈n3 
3。
证明 (1)消元过程乘除次数=2n3+3n2-5n 
6 。
消去xk(k=1,2,…,n-1)一列时,乘法次数=方程个数(n-k)×含右端项时的系
数个数(n+1-k);除法次数=方程个数(n-k)。
又因为
Σn-1 
k=1(n -k)=Σn-1 
k=1
k,Σn 
k=1
k =n(n +1) 
2 ,Σn 
k=1
k2 =n(n +1)(2n +1) 
6 
所以消元过程乘除次数=Σn-1 
k=1((n -k)(n +1-k))+ Σn-1 
k=1(n -k) 
=Σn-1 
k=1(k(k +1))+ Σn-1 
k=1
k =Σn-1 
k=1
k2 +2Σn-1 
k=1
k 
= (n -1)n(2n -1) 
6 +2×n(n -1) 
2 =2n3 +3n2 -5n 
6 
(2)回代过程乘除次数=n(n+1) 
2 。
当回代求xk(k=1,2,…,n)时,乘法次数=n-k;除法次数=1。
所以回代过程乘除次数=Σn 
k=1(n -k)+n =Σn 
k=1(k -1)+n =Σn 
k=1
k =n(n +1) 
2 。
(3)顺序高斯消元法的乘除次数(MD)=消元过程乘除次数+回代过程乘除次数
=2n3+3n2-5n 
6 +n(n+1) 
2 
=n3 
3+n2-n3
≈n3 
3 
证毕。
当方程组的阶数增长时,顺序高斯消元法的运算量增长速度,比克拉默法则求解线性
方程组的运算量增长速度要慢得多。如果线性方程组的阶数=20,那么克拉默法则需要
的乘除次数≈1021,而顺序高斯消元法需要的乘除次数≈3000。
上述程序和算法没有对线性方程组无解和主元为零时求不出解的情况进行判断。
由克拉默法则可知,如果线性方程组的系数行列式det(A)≠0,则线性方程组Ax=b 有
唯一解,否则线性方程组Ax=b 可能无解(如0x +0y=1),也可能有无穷多解(如
0x+0y=0)。
性质3.1 当线性方程组Ax=b 有唯一解时,使用顺序高斯消元法能够求出这个解
的条件是所有的主元素a(k) k,k(k=1,2,…,n)都不为0。

56 　计算方法(Python 版) 
若a(k) k,k =0,则行乘子li,k =-a(k) i,k 
a(k) k,k 
无意义,无法完成消元。
定义3.2 若系数矩阵A = 
a11 a12 … a1n 
a21 a22 … a2n 
. . . . 
an1 an2 … ann 
é
.
êêêêêê 
ù
.
úúúúúú
,则A 的r 阶顺序主子式Δr = 
a11 a12 … a1r 
a21 a22 … a2r 
. . . . 
ar1 ar2 … arr 
(r=1,2,…,n)。
定理3.1 对于n 阶线性方程组Ax=b,顺序高斯消元法能够求出解的充要条件是: 
系数矩阵A 的所有顺序主子式Δi(i=1,2,…,n)均不为0。
证明 (1)因为det(A)=Δn ≠0,所以线性方程组Ax=b 有唯一解。
(2)用数学归纳法证明:Δi(i=1,2,…,n)均不为0.a(k) k,k(k=1,2,…,n)都不为0。
① 当r=1时,Δr=Δ1=a1(11)=a(r) r,r,所以当r=1时,Δi(i=1,2,…,r)均不为0.a(k) k,k 
(k=1,2,…,r)都不为0。
② 假设当r=1,2,…,k 时,Δi (i=1,2,…,r)均不为0.a(k) k,k (k=1,2,…,r)都不
为0。
因为消元过程对应于一系列倍加变换,而倍加变换不改变行列式的值,所以
Δk = 
a11 … a1k 
. . . 
ak1 … akk 
倍加变换(消元)
a1(11) … a1(k1) 
. . . 
0 … a(k) kk 
并且
Δk+1 = 
a11 … a1k a1,k+1 
. . . . 
ak1 … akk ak,k+1 
ak+1,1 … ak+1,k ak+1,k+1 
倍加变换(消元)
a1(11) … a1(k1) a1(1,k)+1 
. . . . 
0 … a(k) kk a(k) k,k+1 
0 … 0 a(k+1) k+1,k+1 
............ 
...................... 
= 
a1(11) … a1(k1) 
. . . 
0 … a(k) kk 
·a(k+1) k+1,k+1 - 
a1(1,k)+1 
. 
a(k) k,k+1 
·0=Δk ·a(k+1) k+1,k+1 
又因为Δk ≠0,所以Δk+1不为0.a(k+1) k+1,k+1不为0,那么r=1,2,…,k,k+1时,Δi(i= 
1,2,…,r)均不为0.a(k) k,k(k=1,2,…,r)都不为0。
③ 由①、②可得,Δi(i=1,2,…,n)均不为0.a(k) k,k(k=1,2,…,n)都不为0。
(3)按性质3.1,对于n 阶线性方程组Ax=b,顺序高斯消元法能够求出解的充要条
件是:系数矩阵A 的所有顺序主子式Δi(i=1,2,…,n)均不为0。
证毕。

第 
3 
章　线性方程组直接求解　

3.列主元高斯消元法
3 

3.3.1 
列主元高斯消元法的主要思想
顺序高斯消元法的缺点如下。

57
(1)当线性方程组Ax=
b 
的系数行列式deA)≠0,且存在一个主元素a(kk)=0时,
t(k,
线性方程组Ax=
b 
有唯一解,但顺序高斯消元法不能求出解。

(2)小主元(主元素的绝对值远小于其他元素的绝对值)时,能够求出解,但求出的解
误差很大。
交换增广矩阵的2行(2个方程的位置互换)或交换2列(2个变元的位置互换), 得到
的方程组与原方程组同解。通过交换行、列来避免主元为0和小主元,然后再消元求解的
方法,称为选主元高斯消元法。列主元高斯消元法是一种常见的选主元高斯消元法,与顺

(k)
序高斯消元法相比,它的改进之处是:在用主元素ak,消去第k+1 行至第
n 
行的变元xk

k 

((

之前,从第
k 
行至第
n 
行所有xk 
的系数(即a(k) k)
k 
,…,k) 中,寻找绝对值最大者

k,aan,
k 
)

k 
,k+1,
(k)
k 
,若ma_i≠k,则交换第
k 
行和第ma_
i 
行,使最大的系数a(k)成为主元
。


amax_i,xx_i,

定理3.2 
当线性方程组Ax=
b 
有唯一解时,列主元高斯消元法必(ma) 能(k) 求出解。
证明用反证法。假设线性方程组Ax=
b 
有唯一解,但列主元高斯消元法求不出

(k)(k)(k)
解。这就是说,在消元过程中,当消去某一个变元(设为xk 
)时,ak,ak 
,…,ak 
全部
为0,即使在其中取绝对值最大者为主元,主元仍然为0。
k 
,k+1,n,
那么系数行列式

(1)(1)

a11 … a1,

k

(1) (1)
a1,
k 
… a1,

n 

...

..................

1 

.. .

a11 a12 … a1n 

(k-1) (k-1)1 a-1,… a-1,

kk 
kn

.................................. 

(k-1)

… ak-k

1,

0

a21 a22 … a2n

det(A)
= 

消去x1~xk-1

.... 

(k)(k)
ak,
k 
… ak,

n

0… 0 

.. . 

... 

((ank,) … ank,)

kn 

an1 an2… ann 

0… 0 

a(1)
11 … a(1)
1,k-1 a(k)k,
k 
… a(k)k,
n 
a(1)1,
k 
… a(1)1,
n 
. . . · . . . -. . . 
0 … a(k-1)k-1,k-1 a(k)n,
k 
… a(k)n,
n 
a(k-1)k-1,
k 
… a(k-1)k-1,
n 
a(1)
11 … a(1)
1,k-1 a(k)k,
k 
… a(k)k,
n 
. . . · . . . 
0 … a(k-1)k-1,k-1 a(k)n,
k 
… a(k)n,
n 

·0

= 

= 

ak(k,) … a(k)kk,

n 

... 

an(k,) … an(k,)kn 

=0。

因为a() ()
k 
,…,()全部为0,所以

k,aan,

kk 
,kk+1,kk 


58 　计算方法(Python 版) 
由上可得,系数行列式det(A)= 
a(1) 
11 … a(1) 
1,k-1 
. . . 
0 … a(k-1) k-1,k-1 
·0=0。
但是线性方程组Ax=b 有唯一解.系数行列式det(A)≠0,这与det(A)=0矛盾, 
所以假设错误,也就是说,当线性方程组Ax=b 有唯一解时,列主元高斯消元法必能求
出解。证
毕。
列主元高斯消元法能够避免小主元时舍入误差被放大的情况,提高解的精度。
例3.1 求解线性方程组0.003x1+3x2=2.001 ① 
x1+x2 =1 ②,设存储精度为4位有效数字。
(精确解:x1=13
,x2=23
,比较顺序高斯消元法和列主元高斯消元法的误差大小) 
解法1 顺序高斯消元法。
①× - 1 
0.003 
.
è .
.
. ÷ 
+②→②,得-999x2=-666,那么x2=-666 
-999≈0.6667(因为存储精
度为4位有效数字,所以x2 舍入误差的绝对值≈3.3×10-5)。
把x2≈0.6667代入①,得x1=(2.001-3×0.6667)/0.003=0.3(因为除以小主元
0.003,所以x1 舍入误差的绝对值≈3.3×10-2,误差被放大1000倍)。
解法2 列主元高斯消元法。
选列主元,得到同解线性方程组
x1+x2 =1 ① 
0.003x1+3x2=2.001 ② 
①× -0.003 
1 
.
è .
.
. ÷ 
+②→②,得2.997x2=1.998,那么x2=1.998 
2.997≈0.6667(因为存储精
度为4位有效数字,所以x2 舍入误差的绝对值≈3.3×10-5)。
把x2≈0.6667代入①,得x1≈1-0.6667=0.3333(因为避免了小主元,所以x1 的舍
入误差没有被放大)。
3.3.2 列主元高斯消元法的算法和程序
用列主元高斯消元法求解n 阶线性方程组,在消元过程中,使用主元a(k) k,k 之前,从第
k 行至第n 行xk 的系数中,寻找最大者a(mka)x_i,k 。若a(mka)x_i,k 为0,则方程组无解,否则交换
第k 行和第max_i 行。
算法3.2 列主元高斯消元法的算法。

第3 章　线性方程组直接求解　 59 
输入原方程组的阶数n 
输入原方程组的增广矩阵a[n][n+1] 
forkinrange(n-1),循环1次消去1列 (消元过程) 
暂设最大系数的值max=a[k][k],最大系数所在行的行号max_i=k 
foriinrange(k+1,n) 
Y |a[i][k]|>|max| N 
令max=a[i][k],max_i=i 
Y max为0 N 
输出"原方程组无解。" 
break 
Y max_i≠k N 
交换行k与行max_i 
foriinrange(k+1,n),循环1次更新1行
行乘子a[i][k]/=-a[k][k] 
forjinrange(k+1,n+1) 
a[i][j]+=a[i][k]*a[k][j] 
Y 直至循环结束仍未执行break N 
forkinrange(n-1,-1,-1)(回代过程) 
s=0 
forjinrange(k+1,n) 
s+=a[k][j]*x[j] 
x[k]=(a[k][n]-s)/a[k][k] 
输出原方程组的解向量x 
程序3.2 列主元高斯消元法对应的程序。 
def inputdata(a,n): 
print("请输入原方程组的增广矩阵: ") 
for i in range(n): 
for j in range(n+1): 
a[i][j]=eval(input('a[{}][{}]='.format(i,j))) 
def outputdata(x,n): 
print("原方程组的解为: ") 
for i in range(n): 
print('x[{}]={}'.format(i,x[i]),end=' ') 
n=int(input("请输入原方程组的阶数: ")) 
a=[[0]*(n+1) for i in range(n)] 
inputdata(a,n) 
for k in range(n-1):

60 　计算方法(Python 版) 
max,max_i=a[k][k],k 
for i in range(k+1,n): 
if abs(a[i][k])>abs(max): 
max,max_i=a[i][k],i 
if max==0: 
print("原方程组无解。") 
break 
if max_i! =k: 
a[k],a[max_i]=a[max_i],a[k] 
for i in range(k+1,n): 
a[i][k]/=-a[k][k] 
for j in range(k+1,n+1): 
a[i][j]+=a[i][k]*a[k][j] 
else: 
x=[0]*n 
for k in range(n-1,-1,-1): 
s=0 
for j in range(k+1,n): 
s+=a[k][j]*x[j] 
x[k]=(a[k][n]-s)/a[k][k] 
outputdata(x,n)
3.4 全主元高斯消元法
3.4.1 全主元高斯消元法的主要思想 
全主元高斯消元法是另一种选主元高斯消元法。全主元高斯消元法与顺序高斯消元
法和列主元高斯消元法相比,它的改进之处是:在用主元素a(k) k,k 消去第k+1行至第n 行
的变元xk 之前,从a(k) k,k 右下方的矩形区域所有xk 的系数中,寻找最大者a(mka)x_i,max_j (即
a(mka)x_i,max_j= max k≤i,j≤n(a(k) i,j )),若max_i≠k,则交换第k 行和第max_i 行;若max_j≠k, 
则交换第k 列和第max_j 列,使最大的系数a(mka)x_i,max_j成为主元。
在程序中不存储变元的名称,确定变元的依据是变元的位置,即变元位于哪一列。因
此当交换两列时,需要记录交换后各变元的位置(分别位于哪一列),求解完毕,再按照原
方程组各变元的次序输出各变元的值。
定理3.3 当线性方程组Ax=b 有唯一解时,全主元高斯消元法必能求出解。
证明略。
与列主元高斯消元法相比,全主元高斯消元法也能够避免小主元时舍入误差被放大
的情况,而且解的精度一般比列主元高斯消元法高。
例3.2 求解线性方程组
x1+x2=2 ① 
2x1+105x2=105 ②,设存储精度为4位数字。
(x1≈-1.0000200004,x2≈0.9999799996,比较列主元高斯消元法和全主元高斯消
元法的误差大小)。
解法1 列主元高斯消元法。
选列主元,得到同解线性方程组:

第3 章　线性方程组直接求解　 61 
2x1 +105x2 =105 
x1 + x2 =2 
①② 
①× -12
.
è .
.
. ÷ 
+②→②,得x2= 
105 
-2 
.
è .
.
. ÷ 
+2 
105 
-2 
.
è .
.
. ÷ 
+1
≈1(因为存储精度为4位有效数字,所以x2 
舍入误差的绝对值≈2.0×10-5) 
把x2=1代入①得:2x1+105×1=105,所以x1=0(因为x2 乘以系数105,所以x1 
舍入误差的绝对值≈1,误差被放大105 倍)。
解法2 全主元高斯消元法。
选全主元,得到同解线性方程组: 
105x2 +2x1 =105 
x2 +x1 =2 
①② 
①× - 1 
105 
.
è .
.
. ÷ 
+②→②,得x1= 
105 
-105 
.
è .
.
. ÷ 
+2 
2 
-105 
.
è .
.
. ÷ 
+1
≈1(因为存储精度为4位有效数字,所以
x1 舍入误差的绝对值≈2.0×10-5) 
把x1=1代入①得:105x2+2×1=105,所以x2=105-2 
105 ≈1(因为回代时避免了乘
以大系数,所以x1 的舍入误差没有被放大)。
与顺序高斯消元法相比,列主元高斯消元法增加的运算量主要如下: 
(1)为寻找最大系数而增加的比较次数≈ Σn-1 
k=1(n -k)=Σn-1 
k=1
k = (n -1)n 
2 ≈n2 
2。
(2)行交换增加的系数交换次数≈ Σn-1 
k=1(n-k+2)=Σn-1 
k=1
k+2(n-1)=(n -1)n 
2 + 
2n -2≈n2 
2。
与顺序高斯消元法相比,全主元高斯消元法增加的运算量主要如下。
(1)为寻找最大系数而增加的比较次数≈ Σn-1 
k=1(n -k +1)2 =Σn-1 
k=1(k +1)2 = 
n(n +1)(2n +1) 
6 -1≈n3 
3。
(2)行交换增加的系数交换次数≈ Σn-1 
k=1(n-k+2)=Σn-1 
k=1
k+2(n-1)=(n -1)n 
2 + 
2n -2≈n2 
2。
(3)列交换增加的系数交换次数≈ Σn-1 
k=1
n =n(n -1)≈n2。
总之,与全主元高斯消元法增加的运算量相比,列主元高斯消元法增加的运算量可以
忽略。全主元高斯消元法的程序结构比列主元高斯消元法要复杂得多。只要原方程组有

2 
　计算方法(Python 
版)

解,这两种方法都能求出解。一般而言,全主元高斯消元法比列主元高斯消元法精度高。

3.4.2 
全主元高斯消元法的算法和程序
在程序中用列表xj 
记录对增广矩阵做一系列交换的过程中各变元的位置(变元位
于哪一列)。增广矩阵的列j(列序号从0开始) xj 
[

j]=j] 
所对应变元的下标=j]。也就是说, 
当xj[
k 
时,解向量中的元素x[存放的是变元xk 
。程序开始时,原方程组各变元
的次序表示为xj 
[增广矩阵的列
j 
移动到列ma_表示为赋值xj 
[

j]=j;xj, max_j]= 
xjxj, j] maxj]

[j]; 交换列
j 
和列ma表示为交换xj 
[和xj 
[_的值。求解完毕,按照原
方程组各变元的次序输出各(_) 变元的值。
算法3.全主元高斯消元法的算法。

3 

输入原方程组的阶数n 
输入原方程组的增广矩阵a[n][n+1], 记录各变元的初始次序xj[j]=j,(j=0,1,2,…,n-1) 
forkinrange(n-1), 循环1次消去1列(消元过程) 

暂设最大系数的值max=a[k][k],max的行号max_i=k,列号max_j=k 
foriinrange(k,n)( 找绝对值最大的系数及其位置) 

forjinrange(k,n) 
Y |a[i][j]|>|max| N 
令max=a[i][j],max_i=i,max_j=j 
Y max为0 N 
输出"原方程组无解。" 
break 
Y max_i≠k N 
交换行k与行max_i 
max_j≠k Y 
交换列k与列max_j 
N 
交换xj[k]与xj[max_j] 
消元模块的2重循环,同列主元高斯消元法,此处略
Y 直至循环结束仍未执行break N 

原方程组无解

回代模块的2重循环,同列主元高斯消元法,此处略
forkinrange(n)( 按原方程组各变元的次序输出各变元的值) 
在列表xj中,元素值为k的元素下标i=xj.nek)

idx(
输出原方程组变元xk的值x[

i] 


第3 章　线性方程组直接求解　 63 
程序3.3 全主元高斯消元法对应的程序。 
def inputdata(a,n): 
print("请输入原方程组的增广矩阵: ") 
for i in range(n): 
for j in range(n+1): 
a[i][j]=eval(input('a[{}][{}]='.format(i,j))) 
def outputdata(x,n,xj): 
print("原方程组的解为: ") 
for k in range(n): 
i=xj.index(k) 
print('x[{}]={}'.format(k,x[i]),end=' ') 
n=int(input("请输入原方程组的阶数: ")) 
a=[[0]*(n+1) for i in range(n)] 
inputdata(a,n) 
xj=[i for i in range(n)] 
for k in range(n-1): #消元过程开始 
max=a[k][k] 
max_i=max_j=k 
for i in range(k,n): 
for j in range(k,n): 
if abs(a[i][j])>abs(max): 
max,max_i,max_j=a[i][j],i,j 
if max==0: 
print("原方程组无解。") 
break 
if max_i!=k: 
a[k],a[max_i]=a[max_i],a[k] 
if max_j!=k: 
for i in range(k,n): 
a[i][k],a[i][max_j]=a[i][max_j],a[i][k] 
xj[k],xj[max_j]=xj[max_j],xj[k] 
for i in range(k+1,n): #真正的消元,循环1 次更新1 行 
a[i][k]/=-a[k][k] 
for j in range(k+1,n+1): 
a[i][j]+=a[i][k]*a[k][j] 
else: 
x=[0]*n #回代过程开始 
for k in range(n-1,-1,-1): 
s=0 
for j in range(k+1,n): 
s+=a[k][j]*x[j] 
x[k]=(a[k][n]-s)/a[k][k] 
outputdata(x,n,xj)

64 　计算方法(Python 版) 
3.5 高斯约当消元法
3.5.1 高斯约当消元法的主要思想 
与顺序高斯消元法的求解过程相比,高斯约当消元法求解过程的特点如下。
(1)在消去变元xk 所在的那一列时,不仅把主元之下的xk 消去,也把主元之上的
xk 消去,也就是把系数矩阵A 除主对角线之外的元素都化为0。
(2)在消去变元xk 所在的那一列时,先把主元a(k) k,k 化为1(让方程k 所有系数都除以
a(k) k,k),再用它去消其他方程的变元xk ,也就是把系数矩阵A 主对角线上的元素都化为1。
消元完毕,原线性方程组同解变换成如下形式: 
x1 =a1(1,n)+1 
x2 =a2(2,n)+1 
. . . 
xn =a(n) n,n+1 
ì
.
í
... 
... 
这时系数矩阵A 被化为单位矩阵,右端向量就是解向量,不再需要回代过程。因此
高斯约当消元法又称为无回代过程消元法。
例3.3 用高斯约当消元法求解线性方程组
2x1- x2-3x3=-2 ① 
2x1-3x2-2x3=-3 ② 
-x1+ x2+ x3=1 ③ 
ì
.
í
.. 
..
。
解 ①×0.5,得 
x1-0.5x2-1.5x3=-1 ① 
2x1- 3x2- 2x3=-3 ② 
-x1+ x2+ x3=1 ③ 
ì
.
í
.. 
.. 
①×(-2)+②→②,①+③→③,得 
x1-0.5x2-1.5x3=-1 ① 
-2x2+ x3=-1 ② 
0.5x2-0.5x3=0 ③ 
ì
.
í
.. 
.. 
②×(-0.5),得 
x1-0.5x2-1.5x3=-1 ① 
x2-0.5x3=0.5 ② 
0.5x2-0.5x3=0 ③ 
ì
.
í
.. 
.. 
②×0.5+①→①,②×(-0.5)+③→③,得 
x1 -1.75x3=-0.75 ① 
x2- 0.5x3=0.5 ② 
-0.25x3=-0.25 ③ 
ì
.
í
.. 
.. 
③×(-4),得 
x1 -1.75x3=-0.75 ① 
x2- 0.5x3=0.5 ② 
x3=1 ③ 
ì
.
í
.. 
.. 
③×1.75+①→①,③×0.5+②→②,得
x1 =1 ① 
x2 =1 ② 
x3=1 ③ 
ì
.
í
.. 
..

第3 章　线性方程组直接求解　 65 
所以原方程组的解为:x1=1,x2=1,x3=1。
3.5.2 高斯约当消元法的算法和程序
与顺序高斯消元法相同,高斯约当消元法也存在主元为0时求不出解,小主元时误差
较大等问题。解决方法与顺序高斯消元法类似,可以采用交换列的方法,保证在有解时一
定能求出解。下面的程序没有对线性方程组无解和主元为零时求不出解的情况进行
判断。算
法3.4 高斯约当消元法的算法。
输入原方程组的阶数n 
输入原方程组的增广矩阵a[n][n+1] 
forkinrange(n),循环1次消去1列
forjinrange(k+1,n+1),把主元ak(k,k)化为1,更新本行的系数
a[k][j]/=a[k][k] 
foriinrange(n),更新除行k之外的各行
Y i≠k N 
forjinrange(k+1,n+1),本循环用来更新本行的系数
a[i][j]+=-a[i][k]*a[k][j] 
输出原方程组的解(此时右端向量就是解向量) 
程序3.4 高斯约当消元法对应的程序。 
def inputdata(a,n): 
print("请输入原方程组的增广矩阵: ") 
for i in range(n): 
for j in range(n+1): 
a[i][j]=eval(input('a[{}][{}]='.format(i,j))) 
def outputdata(a,n): 
print("原方程组的解为: ") 
for i in range(n): 
print('x[{}]={}'.format(i,a[i][n]),end=' ') 
n=int(input("请输入原方程组的阶数: ")) 
a=[[0]*(n+1) for i in range(n)] 
inputdata(a,n) 
for k in range(n): #循环1 次消去1 列 
for j in range(k+1,n+1): #把主元a(k) 
k,k 化为1,更新本行的系数 
a[k][j]/=a[k][k] 
for i in range(n): #更新除行k 之外的各行 
if i!=k: 
for j in range(k+1,n+1): #本循环用来更新本行的系数

66 　计算方法(Python 版) 
a[i][j]+=-a[i][k]*a[k][j] 
outputdata(a,n) 
3.5.3 一次求解出多个线性方程组
行乘子li,k 的计算公式li,k =-a(k) i,k 
a(k) k,k 
中,不出现右端向量b。当用高斯约当消元法求
解多个系数矩阵A 相同,而右端向量b 不同的线性方程组时,消元过程对应的一系列倍
加变换是相同的。也就是说,以什么次序,哪一行(方程)乘以什么因子(行乘子),加到哪
一行(方程)上,不受右端向量b 的影响。因此,可以把这些方程组的右端向量合并到一
个增广矩阵中,每多一个线性方程组,增广矩阵就多增加一列。当对增广矩阵做倍加变
换,把系数矩阵A 化为单位矩阵时,右端向量就被化为对应方程组的解。这样,原来求解
一个线性方程组的倍加变换,可以把这些线性方程组全部求出,减少了总的运算量。
例3.4 用高斯约当消元法,一次消元求解下列多个线性方程组: 
① 
2x+y-z=2 
- x+ 3z=2, 
-2x+y+z=0 
ì
.
í
.. 
..
② 
2x+y-z=1 
- x+ 3z=8, 
-2x+y+z=3 
ì
.
í
.. 
..
③ 
2x+y-z=7 
- x+ 3z=0 
-2x+y+z=-3 
ì
.
í
.. 
.. 
解 构造方程组①、②、③的增广矩阵,得
..A = [A b] = 
2 1 -1 2 1 7 
-1 0 3 2 8 0 
-2 1 1 0 3 -3 
........ 
é
.
êêêê 
ù
.
úúúú 
r1/2→
1 0.5 -0.5 1 0.5 3.5 
-1 0 3 2 8 0 
-2 1 1 0 3 -3 
........ 
é
.
êêêê 
ù
.
úúúú 
r1 +r2 r1 ×2+r3→
1 0.5 -0.5 1 0.5 3.5 
0 0.5 2.5 3 8.5 3.5 
0 2 0 2 4 4 
........ 
é
.
êêêê 
ù
.
úúúú 
r2 ×2→ 
1 0.5 -0.5 1 0.5 3.5 
0 1 5 6 17 7 
0 2 0 2 4 4 
........ 
é
.
êêêê 
ù
.
úúúú 
r2/(-2)+r1 r2 × (-2)+r3→
1 0 -3 -2 -8 0 
0 1 5 6 17 7 
0 0 -10 -10 -30 -10 
........ 
é
.
êêêê 
ù
.
úúúú 
r3/(-10)→
1 0 -3 -2 -8 0 
0 1 5 6 17 7 
0 0 1 1 3 1 
........ 
é
.
êêêê 
ù
.
úúúú 
r3 ×3+r1 r3 × (-5)+r2→
1 0 0 1 1 3 
0 1 0 1 2 2 
0 0 1 1 3 1 
........ é
.
êêêê 
ù
.
úúúú 
所以方程组①的解为
x=1 
y=1 
z=1 
ì
.
í
.. 
..
,方程组②的解为
x=1 
y=2 
z=3 
ì
.
í
.. 
..
,方程组③的解为
x=3 
y=2 
z=1 
ì
.
í
.. 
..
。
3.5.4 一次求解多个线性方程组的算法和程序
算法3.5 用高斯约当消元法一次消元求解多个线性方程组的算法。

第3 章　线性方程组直接求解　 67 
输入原方程组的阶数n和原方程组的个数m 
输入合并之后的增广矩阵a[n][n+m] 
forkinrange(n),循环1次消去1列
forjinrange(k+1,n+m),把主元ak(k,k)化为1,更新行k的系数
a[k][j]/=a[k][k] 
foriinrange(n),更新除行k之外的各行
Y i≠k N 
forjinrange(k+1,n+m),本循环用来更新本行的系数
a[i][j]+=-a[i][k]*a[k][j] 
输出这些方程组的解(每一列右端向量对应一个方程组的解向量) 
程序3.5 用高斯约当消元法一次消元求解多个线性方程组对应的程序。 
def inputdata(a,n,m): 
print("请输入合并之后的增广矩阵({}行、{}列): ".format(n,n+m)) 
for i in range(n): 
for j in range(n+m): 
a[i][j]=eval(input('a[{}][{}]='.format(i,j))) 
def outputdata(a,n,m): 
for k in range(m): 
print("\n 方程组{}的解为: ".format(k)) 
for i in range(n): 
print('x{}[{}]={}'.format(k,i,a[i][n+k]),end=' ') 
n=int(input("请输入原方程组的阶数: ")) 
m=int(input("请输入原方程组的个数: ")) 
a=[[0]*(n+m) for i in range(n)] 
inputdata(a,n,m) 
for k in range(n): #循环1 次消去1 列 
for j in range(k+1,n+m): #把主元a(k) 
k,k 化为1,更新行k 的系数 
a[k][j]/=a[k][k] 
for i in range(n): #更新除行k 之外的各行 
if i!=k: 
for j in range(k+1,n+m): #本循环用来更新本行的系数 
a[i][j]+=-a[i][k]*a[k][j] 
outputdata(a,n,m) 
3.6 消元形式的追赶法
3.6.1 消元形式的追赶法的主要思想 
带状对角矩阵是一种常见的高阶稀疏矩阵。3对角矩阵是带状对角矩阵的一种。在

68 　计算方法(Python 版) 
做三次样条插值、求解微分方程等问题时,经常需要解3对角方程组。
定义3.3 如果在方阵A 中,除主对角线和两条次对角线之外的元素都为0,那么A 
是3对角矩阵,即A 可以表示为下面的形式: 
A = 
a1,0 a1,1 
a2,-1 a2,0 a2,1 
. . . 
an-1,-1 an-1,0 an-1,1 
an,-1 an,0 
é
.
êêêêêêê 
ù
.
úúúúúúú 
其中,ai,j的左下标i 为行号,右下标j 为0表示ai,j位于主对角线,j 为1表示ai,j位于右
上次对角线,j 为-1表示ai,j位于左下次对角线。系数矩阵A 中各个系数下标之间的关
系是,ai,j的下方为ai+1,j-1,ai,j的右方为ai,j+1。
定义3.4 如果线性方程组Ax=b 的系数矩阵A 是3对角矩阵,那么Ax=b 是3对
角方程组。
定义3.5 如果在方阵A 中,除主对角线、左下lw 条次对角线和右上uw 条次对角
线之外的元素都为0,那么A 是带状对角矩阵,称lw 为下带宽,uw 为上带宽,bw =uw + 
lw +1为带宽,即A 可以表示为下面的形式: 
A = 
a1,1 … a1,uw+1 
. . . 
alw+1,1 . an-uw ,n 
. . . 
an,n-lw … an,n 
é
.
êêêêêêê 
ù
.
úúúúúúú 
其中,左下标为行号,右下标为列号。
定义3.6 如果线性方程组Ax=b 的系数矩阵A 是带状对角矩阵,那么Ax=b 是带
状对角方程组。
追赶法用于求解3对角方程组。在求解3对角方程组时,原来的求解过程被大大简
化。要求解的3对角方程组可以表示为
a1,0x1 + a1,1x2 =a1,2 
a2,-1x1 + a2,0x2 +a2,1x3 =a2,2 
. . . . 
an-1,-1xn-2 + an-1,0xn-1 + an-1,1xn =an-1,2 
an,-1xn-1 + an,0xn =an,2 
ì
.
í
.... 
...
. 
求解3对角方程组,不需要存储系数矩阵A 中除3条对角线之外的零元素,只需要
存储3条对角线、右端向量和解向量,使存储空间需求从O(n2)降为O(n)。
这里按照顺序高斯消元法的求解思路求解这个3对角方程组。在消元过程中,当消
去xi 一列时,主元素为ai,0。ai,0下方只有一个系数ai+1,-1非0,因此当消去xi 一列时, 
只需要更新第i+1行中变元的系数。这一行中,也只有ai+1,0和ai+1,1需要更新。