高斯消元法 1

求解线性方程组的直接解法

5.1 Gauss 消去法

① 三角方程组

先举一个简单的例子来说明消去法的基本思想.

例1. 用消去法解方程组

(1)⎧x 1+x 2+x 3=6,

⎪

(2) ⎨4x 2-x 3=5,

⎪2x -2x +x =1. (3)

23⎩1

解 第一步. 将方程(1)乘上-2加到方程(3)上去, 消去(3)中的未知数x 1, 得到 -4x 2-x 3=-11. ( 4 )

第二步. 将方程(2)加到方程(4)上去, 消去方程(4)中的未知数x 2, 得到与原方程组等

价的三角形方程组

⎧x 1+x 2+x 3=6, ⎪

(5 ) ⎨4x 2-x 3=5,

⎪-2x =-6.

3⎩

显然, 方程组(5)是容易求解的, 解为x *=(1, 2, 3) T .

上述过程相当于

⎛1116⎫⎛1116⎫⎛1116⎫ ⎪ ⎪ ⎪

(A |b ) = 04-1 5⎪→ 04-1 5⎪→ 04-1 5⎪

2-211⎪ 0-4-1-11⎪ 00-2-6⎪⎝⎭⎝⎭⎝⎭

(-2)⨯r 1+r 3→r 3 r 2+r 3→r 3其中用r i 表示矩阵的第i 行.

下面我们讨论求解一般线性方程组的高斯消去法. 一般地

⎧a 11x 1+a 12x 2+ +a 1n x n =b 1⎪a x + +a x =b ⎪2222n n 2

⎨

⎪ ⎪⎩a nn x n =b n

当a 11a 22…a nn ≠0时, 可解出 x n =b n /ann for k=n-1:1 x k =(b 1- ak,k+1x k +1-…- akn x n )/ akk end

注: x k , b k 可用同一组单元. 并可解出一个未知数即代入其它方程消去该未知数

Gauss 消元法的流程图为: 流程图中，a ij , b i k 是循环次数。

(i , j 1,2,..., n ) 分别为线性方程组的系数矩阵和常数向量；

② 顺序消去法

一般地, k =1对n 阶方程组消去第k 个元(a kk ≠0):

⎡a kk a k , k +1⎢

⎢a k +1, k a k +1, k +1⎢ ⎢

a n , k +1⎢⎣a nk

a kn

a k +1, n

a nn

a k , k +1b k ⎤⎡a kk

⎥⎢

b k +1⎥⎢(m k +1, k ) a k +1, k +1

→

⎥⎢

⎥⎢b n ⎥⎦⎢⎣(m n 1) a n , k +1

a kn

a k +1, n

a nn

b k ⎤⎥b k +1⎥ ⎥⎥b n ⎥⎦

这里各行变换:i 行-k 行×m ik , 其中m ik =a ik /a kk , i =k +1,…,n 而后, k=2对n -1阶方程组

消第2个元…我们有如下顺序消元算法:

for k=1:n -1 if a kk ≠0 for i =k +1:n m ik =aik /a kk i 行=i行-k 行×m ik end else stop end end

（每行包括右端项! ）可细化, 也可存储m ik 于a ik 得:

for k=1:n -1 if a kk ≠0 for i =k +1:n a ik =aik /akk for j=k+1:n a ik =a ik -a ik ×a kj end b i =bi -a ik ×b k end else stop end end

顺序消元过程和回代过程连起来就可得精确解. 顺序消元算法也可将系数矩阵和右端项分开:

for k=1:n -1 if a kk ≠0 for i =k +1:n a ik =aik /akk for j=k+1:n a ik =a ik -a ik ×a kj end end else stop end end (注意m ik 在a ik ) for k=1:n -1 for i =k +1:n b i =bi -a ik ×b k end

end

Gauss 消去法运算量

消去第k 个元素时，对矩阵作加法和乘法运算各(n-k ) ×(n-k ) 次, 除法(n-k ) 次. 对右端作加法和乘法运算各(n-k ) 次. 分别共12+22+…+(n -1) 2=n (n -1/2)(n -1)/3和1+2+…+(n-1)=n (n -1)/2次加法乘法, 消元时还有1+2+…+(n-1)= n(n -1)/2次除法. 另外回代过程中加法和乘法运算各n (n -1)/2次, 除法n 次. 运算量主要是消元的贡献, 加法和乘法运算各约n 3/3. 定理1. 设Ax=b, 其中A ∈R n ⨯n .

(k )

(1) 如果a kk ≠0 (k=1,2, , n ), 则可通过Gauss 消去法将Ax=b约化等价的三角形(1) (1) )

⎡a 11⎤x ⎡b 1(1) ⎤a 12 a 1(1n ⎡1⎤⎢⎥⎢(2) ⎥(2) (2) ⎥⎢x a a 2222n ⎥⎢⎥=⎢b 2⎥, 方程组⎢⎢ ⎥⎢ ⎥⎢ ⎥⎢⎥⎢(n ) ⎥(n ) ⎥⎢x a ⎢⎥⎢nn ⎦⎣m ⎦⎣⎣b n ⎥⎦

且计算公式为:

(a) 消元过程

(k ) 设a kk ≠0, 对k =1, 2, , n -1计算

(k ) (k )

⎧m ik =a ik /a kk ⎪(k +1)

(k ) (k )

⎨a ij =a ij -m ik a kj

⎪(k +1) (k ) (k )

=b i -m ik b k ⎩b i

i , j =k +1, k +2, , n

(b) 回代过程

(n ) (n )

⎧x n =b n /a nn ⎪n ⎨(i ) (i ) (i )

x =(b -a x ) /a ∑i i ij j ii ⎪

j =i +1⎩

i =n -1, , 2, 1

(2) 如果A 为非奇异矩阵, 则可通过Gauss 消去法(及交换两行的初等变换) 将方程组Ax=b约化为

(1) ⎡a 11⎢⎢⎢⎢⎢⎣

(1) )

⎤x ⎡b 1(1) ⎤a 12 a 1(1n ⎡1⎤

⎥⎢(2) ⎥(2) (2) ⎥⎢x a 22 a 2b n ⎥⎢2⎥=⎢2⎥.

⎥⎢ ⎥⎢ ⎥

⎥⎢(n ) ⎥(n ) ⎥⎢x a nn ⎥⎦⎣m ⎦⎢⎣b n ⎥⎦

行列式和逆矩阵

易知顺序消元过程中，行列式不变. 因此det (A )=det (U )=u 11u 22…u nn , 这里U 是顺序消元过程结束时的上三角矩阵. A 的顺序主子式。也是同样情况, 即有det (A k )=det (U k )=u 11u 22…u kk . 。顺便指出, 消第k 个元素时，左上角的元素(称之为

主元素) 非零时才能进行. 所有主元素非零等价于det (A k ) ≠0, k =1,2,…, n -1.

求A 的逆矩阵相当于解n 个方程组, 它们的系数矩阵是A , n 个右端项是n 个单位向量. 用增广矩阵表示即(A ∣I ) ，其中I 是单位矩阵. 消元过程可对所有右端项同时进行, 消元过程结束再分别回代. 也可巧妙安排就地求逆. ③ 主元素法

如果a kk =0顺序消元无法进行, 如果a kk 很小，可以进行。但将引出很大误差. 例2. 用顺序消去法解方程组

⎡0. 0001⎢1⎣

解 用精确运算:

1⎤⎡x 1⎤⎡1⎤

=⎢⎥ ⎢⎥⎥1⎦⎣x 2⎦⎣2⎦

⎡0. 00011⎤⎡0. 000111⎤

→⎢⎥⎢⎥

120--9998⎣⎦⎣⎦

解得(10000/9999,9998/9999)≈(1.[1**********]100,0.[1**********]900)

若在十进三位尾数舍入的浮点计算机系统中运算, 第二行将是

(0 -10000 ∣ -10000)

从而得到解x 1=0, x 2=1, 与实解相去甚远. 究其原因, 小主元作分母扩大误差所致. 如果两个方程(两行) 交换次序再消元:

⎡0. 000111⎤⎡112⎤⎡12⎤⎢⎥→⎢⎥→⎢⎥1120. 00011101⎣⎦⎣⎦⎣⎦

得到解x 1=1, x 2=1, 与实解很近.

下面引进的列主元素法在消元前，先选该列中绝对值最大的做主元(交换两

行, 每行包括右端项!): for k=1:n -1 找p :

a pk =max a kk , a k +1, k , , a nk

)

p 行k 行 i k =p

if a kk ≠0

for i =k +1:n m ik =aik /a kk i 行=i 行-k 行×m ik end else stop end

end

i k =p记录每次选主元(两行交换) 的指标. 跟前面一样, 算法还可细化, 把矩阵和右端分开算(这时就要用指标i k ) 。列主元素法乘数m ik 绝对值不大于1, 不会增加误差。列主元素法用来求行列式时，要注意两行交换行列式变号。还有一种全主元素法, 它是在整个右下(n-k ) ×(n-k ) 矩阵找绝对值最大的做主元(交换行及列). 这对误差控制有利, 但搜索太费时. 通常列主元素法误差控制就已可以了. 以下是列主元素法的流程图:

流程图中，a ij , b i

(i , j 1,2,..., n ) 分别为线性方程组的系数矩阵和常熟

向量；k 是循环次数；d 是主元素；

二 实验部分

本章实验内容：

实验题目：Gauss 消元法, 追赶法, 范数。 实验内容：①编制用Gauss 消元法求解线性方程组Ax=f的程序。

②编制用追赶法求解线性方程组Ax=f的程序。

③编制向量和矩阵的范数程序。

实验目的：①了解Gauss 消元法原理及实现条件, 熟练掌握Gauss

消元法解方程组的算法, 并能计算行列式的值。

②掌握追赶法，能利用追赶法求解线性方程组。

③理解向量和矩阵范数定义, 性质并掌握其计算方法.

编程要求：利用Gauss 消元法, 追赶法解线性方程组。分析误差。 计算算法:①Gauss 消元法:

1. 消元过程

(k )

设a kk ≠0, 对k =1, 2, , n -1计算

(k ) (k ) ⎧m ik =a ik /a kk ⎪(k +1)

(k ) (k )

⎨a ij =a ij -m ik a kj ⎪(k +1)

=b i (k ) -m ik b k (k ) ⎩b i

i , j =k +1, k +2, , n

⒉回代过程

(n ) (n )

⎧x n =b n /a nn ⎪n ⎨(i ) (i ) (i )

x =(b -a x ) /a ∑i i ij j ii ⎪

j =i +1⎩

i =n -1, , 2, 1

②追赶法：

1. 分解Ax=f: β1=c 1/b 1, βi =c i /(b i -a i βi -1)

( i =2, 3, , n -1)

2. 解Lg=f,求g:

g 1=f 1/b 1, g i =(f i -a i g i -1) /(b i -a i βi -1) (i =2, 3, , n )

3．解Ux=g,求x:x n =y n , x i =y i -βi x i +1 （i =n -1, n -2, , 2, 1） ③范数：

常用向量范数有：（令x =( x1, x 2, …, x n ) T ）

1-范数： ║x ║1=│x 1│+│x 2│+…+│x n │

2-范数： ║x ║2=(│x 1│+│x 2│+…+│x n │)

2

2

2

1/2

∞-范数： ║x ║∞=max(│x 1│,│x 2│,…,│x n │)

常用的三种向量范数导出的矩阵范数是：

1-范数:║A ║1= max{║Ax ║1/║x ║1=1}=max ∑a jj

1≤j ≤n

i =1n

2-范数:║A ║2=max{║Ax ║2/║x ║2=1}=1, λ1是A T A 的

最大特征值.

∞-范数:║A ║∞=max{║Ax ║∞/║x ║∞=1}=max ∑a ij

1≤i ≤n

j =1n

实验例题⑴：用Gauss 消元法解方程组

⎛2-11⎫⎛x 1⎫⎛4⎫ ⎪ ⎪ ⎪-1-23 ⎪ x 2⎪= 5⎪ 1 ⎪ ⎪31⎪⎝⎭⎝x 3⎭⎝6⎭

实验例题⑵: 用追赶法解三对角方程组Ax=b, 其中

00⎤⎡2-10⎡1⎤⎢-12-10⎥⎢0⎥0⎢⎥⎢⎥A =⎢0-12-10⎥, b =⎢0⎥ ⎢⎥⎢⎥00-12-1⎢⎥⎢0⎥⎢⎢00-12⎥⎣0⎦⎣0⎥⎦

实验例题⑶：设

⎛0. 60. 5⎫A = 0. 10. 3⎪⎪， ⎝⎭

计算A 的行列范数.

程序①:Gauss消元法

function x=Gauss(A,b)

%A是线性方程组的系数矩阵, b 为自由项. n=length(A) for k=1:n-1

m(k+1:n,k)=A(k+1:n,k)/A(k,k);

A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m(k+1:n,k)*A(k,k+1:n); b(k+1:n)=b(k+1:n)-m(k+1:n,k)*b(k); end

x=zeros(n,1); x(n)=b(n)/A(n,n); for k=n-1:-1:1

x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);

end x=x';

disp(sprintf('k x(k)')); for i=0:n

disp(sprintf('%d %f ',i,x(i+1))); end

数值结果：x=Gauss(A,b)

n =3

k x(k)

0 1.111111 1 0.777778 2 2.555556

程序②:追赶法

function [x,y,beta]=zhuiganfa(a,b,c,f)

%a,b,c是三对角阵的对角线上的元素, f 是自由项. n=length(b);

beta(1)=c(1)/b(1); for i=2:n

beta(i)=c(i)/(b(i)-a(i)*beta(i-1)); end

y(1)=f(1)/b(1); for i=2:n

y(i)=(f(i)-a(i)*y(i-1))/(b(i)-a(i)*beta(i-1)); end

x(n)=y(n); for i=n-1:-1:1

x(i)=y(i)-beta(i)*x(i+1); end

disp(sprintf('k x(k) y(k) beta(k)')); for i=0:n

disp(sprintf('%d %f ',i,x(i+1),y(i+1),beta(i+1))); end

数值结果:

a=[0 -1 -1 -1 -1]'; b=[2 2 2 2 2]'; c=[-1 -1 -1 -1 0]'; f=[1 0 0 0 0]';

[x,y,beta]=zhuiganfa(a,b,c,f)

k x(k) y(k) beta(k) 0 0.833333 5.000000e-001 -0.500000 1 0.666667 3.333333e-001 -0.666667 2 0.500000 2.500000e-001 -0.750000 3 0.333333 2.000000e-001 -0.800000 4 0.166667 1.666667e-001 0.000000

程序③： 1.列范数:

function fan=lie(A) %A为已知矩阵 n=length(A) for j=1:n x(j)=0 for i=1:n

x(j)=x(j)+abs(A(i,j)); end end

fan=max(x)

disp(sprintf('n x(n)')); for i=0:n

disp(sprintf(' %d %f',i,x(i+1))); end 数值结果: fan=lie(A) fan =0.8000

n x(n) 0 0.700000 1 0.800000

2. 行范数:

function fan=hang(A) %A为已知矩阵 n=length(A) for i=1:n x(i)=0 for j=1:n

x(i)=x(i)+abs(A(i,j)); end end

fan=max(x)

disp(sprintf('n x(n)')); for i=0:n

disp(sprintf(' %d %f',i,x(i+1))); end

数值结果: fan=hang(A) fan =1.1000 n x(n) 0 1.100000 1 0.400000

总结:从代数上看, 直接分解法和Gauss 消去法本质上一样, 但如果我们

采用”双精度累加”计算 a i b i , 那么直接三角分解法的精度要比Gauss 消去法为高.

求线性方程组的直接法, 其算式繁杂, 给人以枯燥沉闷的感觉. 为了改善教学效果, 本章着重介绍了三对角方程组的追赶法. 三对

角方程组以及其拓广形式的带状方程组有着广泛的实际应用.

追赶法是解三对角线方程组(对角元占优势) 的有效方法, 它

具有计算量少, 方法简单, 算法稳定等优点, 具有鲜明的对称美.

复习思考题

1. 用消去法解线性方程组为什么最好选主元？怎样的方程组可以不用选主元？

2. 用高斯—约当消去法求矩阵的逆，其理论根据是什么？

3. 求矩阵A 的LU 分解的紧凑格式有何规律？写出LU 分解法解Ax = b 的算法。

4. 乔累斯基分解法适用于哪类线性方程组？有何优点？

5. 追赶法有何优点？哪类方程组可用追赶法求解？

6. 向量范数与矩阵范数是如何定义的？各有何性质？常用范数有哪些？

7. 什么叫病态方程组？什么叫矩阵的条件数？条件数刻画了矩阵的什么性质？能否用选主元的方法求得病态方程组的高精度的解？