线性方程组求解的数值方法(讲义)
第三章 线性方程组求解
由于很多数学和工程问题都可以转化为线性方程组问题,线性方程组的数值求解具有重要意义。线性方程组解包括无解、有唯一解和有无穷解的三种情况。对于无解问题,工程上一般通过求解最小二乘问题,得到其近似解;对于无穷解问题,需要通过计算该问题对应齐次线性方程的基础解析和一个特解,构造问题的通解。本章主要讨论方程解存在且唯一情况下,求解该方程解的方法。
3.1 线性方程组回顾
3.1.1 线性方程组的概念
线性方程组的标准形式如下式所示:
a11x1a12x2axax211222am1x1am2x2a1nxnb1,a2nxnb2,amnxnbm, (3.1)
其中,x1,x2,...xn为未知量,a11,...,ann为线性方程组的系数,b1,b2,...,bm为常数项。
利用线性代数中矩阵的相关知识,线性方程Amn组也可表示为矩阵的形式:
Amnxn1bm1 (3.2)
其中,为线性方程组的系数矩阵,xn1为未知向量,bm1为常数向量。如果
bm10,则称线性方程组为”齐次线性方程组”,否则称为”非齐次线性方程组”。
在实际中,许多理论、工程问题都可以转化为线性方程组的问题加以解决。例如,在电路分析中,根据Kirchhoff电流定律,电路中任意节点的输入电流和输出电流的代数和为零。
2Ω
3Ω3Ω4Ω2Ω
1Ω
图3-1 电路模型
以此可建立如下线性方程组:
55000i1500111i0200023i30 11100i40057200i5 (3.3)
在数值分析中,解决函数的数值逼近以及离散微分方程求解等问题也需要用到线性方程组求解的相关知识。
另外,有些问题虽然形式上不属于线性方程组求解问题,但通过简单的转化可以转化为线性方程求解问题,并加以解决。
例3.1:求解方程:
e2x13x22x310x1x23x31 e
e3x210x36
解:通过将等式两边取对数,得到:
2x13x22x3ln(10)x1x23x30 3x10xln(6)23
该方程为线性方程组,可很容易的利用本章所述方法求解。
例3.2:求解方程:
(2x13x2)21 x1x21e
解:该方程可等价转化为:
2x13x212x13x21和, x1x20x1x20
从而得到方程的两组解。
一般的,对于f(a1x1...anxn)g的方程,都可通过先求解外层函数的解,然后再求内层线性方程解的形式加以解决。
3.1.2 线性方程组解的情况
在讨论线性方程组求解方法之前,首先回顾线性方程组解的情况。
对于任意线性方程组,其解可能包括三种情况:无解、有唯一解和有无穷多解。线性方程组如果有解,就称它是相容的;如果无解,就称它不相容。
利用系数矩阵A和增广矩阵B(A,b)的秩,可以方便的讨论线性方程组是否有解(即是否相容)以及解是否唯一的问题。
定理3.1 n元线性方程组Axb
无解的充分必要条件是R(A)R(A,b);
有唯一解的充分必要条件是R(A)R(A,b)n;
有无限多解的充分必要条件是R(A)R(A,b)n。
由定理3.1容易得出线性方程组理论中两个最基本的定理。
定理3.2 n元齐次线性方程组有非零解的充分必要条件是R(A)n。 定理3.3 n元线性方程组有解的充分必要条件是R(A)R(A,b)。
对于齐次方程组Ax0,其解向量满足如下性质:
性质1 若x1,x2为齐次方程组Ax0的解,则x1+2也是该方程的解。
性质2 若x1为齐次方程组Ax0的解,k为实数,则xk1也是该方程的解。
根据上述两条性质,可知所有满足公式:
xk11k22kii (3.4)
的向量,都是该方程的解。齐次线性方程组的解集的最大无关组称为该齐次线性方程组的基础解系。
对于非齐次线性方程Axb,则具有如下性质:
性质3 设x1及x2都是非齐次线性方程Axb的解,则x12为对应的齐次线性方程组Ax0的解。
性质4 设x是非齐次线性方程Axb的解,x是对应的齐次线性方程组Ax0的解,则x仍是该非齐次方程的解。
于是非齐次线性方程Axb的通解为:
xk11knrnr*(k1,knr为任意实数)。
其中,1,nr是齐次线性方程组Ax0的基础解系。
Matlab提供了求解线性方程组的方法:
1、求解齐次线性方程组Ax0基础解系的“null”函数
2、求解非齐次线性方程组Axb的特解的函数“\”(左除)
利用上述两个函数,可以很方便的计算非齐次线性方程组的通解。
例3.3:求解方程
123x14
123
3x4
2
12x34
的通解。
解:很容易判断,该方程的系数矩阵与增广矩阵的秩相对,且小于n,因此方程具有无穷多解。
为了求解该方程的通解,首先在matlab中输入系数矩阵:
然后,利用null函数,得到该方程对应齐次方程组的基础解系:
其次,输入常数向量,并利用“\”,计算非齐次方程组的特解。
至此,该非齐次方程组的通解可表示为:
00.96364
xk1-0.8321k2-0.14820 0.5547-0.22240
3.2 消元法
消元法主要包括高斯消去法和对其的改进方法,如列主元消去法和全主元消去法。消去法基于如下简单的事实:三角矩阵可以很容易的求出该方程的解。
例3.4:求解方程
345x18012x2 2001x31
解:该方程可等价表示为:
3x14x25x38 x22x32
x13
通过自下而上依次代入,可很容易得到该方程的解为:
x31
x2b22x3220
x1b14x25x3805133
消元法的基本思路是通过一系列的等价变换,将原方程转化为与之等价的三角矩阵线性方程组,然后得到原方程的解。
3.2.1 高斯消去法
在消去法解方程组的过程中,需要始终保持所变换方程与原方程等价。高斯消去法通过一系列初等行变换,将原线性方程组问题转化为三角矩阵方程求解问题,从而简化求解过程,得到原方程的解。
初等行变换包含如下三种变换:
对调两行(对调i,j两行,记作rirj);
以数k0乘某一行中的所有元素(第i乘数k,记作rik);
把某一行所有元素的k倍加到另一行对应的元素上去(第j行的k倍加到
第i行上,记作rikrj)。
把定义中的“行”换成“列”,即得矩阵的初等列变换的定义(所用记号是把“r”换成“c”)。矩阵的初等行变换与初等列变换,统称初等变换。
设A是一个mn矩阵,对A施行一次初等行变换,相当于在A的左边乘以
相应的m阶初等矩阵;对A施行一次初等列变换,相当于在A的右边乘以相应的n阶初等矩阵。因此,初等行变换可表示为:PrA;初等列变换可表示为:APc。
如果矩阵A经过有限次初等行变换变成矩阵B,就称矩阵A与B行等价,记r作A如果矩阵A经过有限次初等列变换变成矩阵B,就称矩阵A与B列B;
c等价,记作A如果矩阵A经有限次初等变换变成矩阵B,就称矩阵A与B;
B等价,记作AB。
矩阵之间的等价关系具有下列性质:
(i)反身性 AA
(ii)对称性 若AB,则BA;
(iii)传递性 若AB,BC,则AC。
不难证明:对于任何矩阵Ann,总可以经过有限次初等行变换将其转化为三角矩阵的形式。
因此,高斯消去法过程的矩阵表示为:
rAxbP1AxPb1
r P2P1AxP2Pb1
......
r Uxb'r
***
其中,U表示上三角矩阵,U0...*。 00*
高斯消去法具有如下性质:
性质1:若A的所有顺序主子式均不为0,则高斯消元无需换行即可进行到底,得到唯一解。
性质2:只要 A 非奇异,即A1存在,则可通过逐次消元及行交换,将方程组化为三角形方程组,求出唯一解。
例3.5:求解方程:
1070x17326x4 2515x36
解:首先写出该方程对应的增广矩阵:
1070 326
515 74
6
通过初等行变换,将矩阵转化为:
1070 00.16
02.55 76.1
2.5
继续利用初等行变换,将矩阵进一步转化为:
0 10700.16
0155 076.1
155
最后得到该方程的解为:
x31
6.16x31 0.1
77x0xx1010x2
在高斯消去法过程中,可能出现对角线主元为0的情况,还需要通过换行、换列等操作,灵活选择选择主元。
选择策略包括:列主元消去法和全主元消去法。
列主元消去法通过初等行变换,每次仅选子矩阵中第k列绝对值最大的元素作为主元进行下一次消去运算,其主元选择公式:
(k1)(k1)akmaxa,ki,k0;k2,...,nkin (3.5)
全主元消去法每一步选子矩阵中绝对值最大的元素为主元素,进行下一次消去运算,其主元选择公式:
注意,由于全主元消去法在变换过程中进行了初等列变换,未知向量各元素间的相对关系将发生改变,需要在反向求解时加以纠正。
例3.6:求解方程 (k1)(k1)0;k2,...,nakmaxmaxa,ki,jkinkjn(3.6)
113x12331x6 2110x33
解:首先写出增广矩阵:
113 331
110 26
3
对其进行一步消元,得到:
113 00-8
02-3 20
1
由于对角线元素为0,进一步消元无法继续,需要进行换行处理,得到:
113 02-3
00-8 21
0
最后得到方程解为:
x30
13x3 0.52
x12x23x32.5x2
例3.7:求解方程
113x12331x6 2110x33
解:首先写出增广矩阵:
113 331
110 26
3
对其进行一步消元,得到:
113 00-8
02-3 20
1
由于对角线元素为0,进一步消元无法继续,子矩阵的最大元素-8,并进行列交换,得到:
131 0-80
0-32 20
1
进一步消元,得到:
131 0-80
002 20
1
需要注意的是,在反向计算各未知变量值时,由于进行了2、3列间的列交换,变量x2和x3的位置发生交换。因此,该方程的解为:
x20.5
x30
x123x3x22.5
另外,在实际计算汇总,由于计算机的精度游侠,即使主元素不为0,但是比较小,也应该换行。
例3.8,假设某计算机精度为3位小数,求解如下线性方程组,
1041x11x2 112
解:方法一:不进行换行,直接计算。
1041 11
10041104201 11 -0.9999 0.99981-1.000 x101.000x21
方法二:先进行行交换,在进行计算:
211 2100.9999 0.9998 1 2x111
0-1.000 1.000x2111 41 10
很容易得到该方程的精确解为:
x1
1
,x22x1
1104
可以看出,如果不进行换行,选择适当的主元,计算结果较实际结果将产生很大的偏差。
最后,给出高斯消去法过程的描述:
步骤 1:将方程Axb写为增广矩阵的形式:
1
(1)2 AAb
...n
(1)(1)
步骤 2:设a110,计算第i行的消去因子: mi(1)ai(1)/a111
(i2,...,n)
步骤 3:将增广矩阵A第i行减去mi(1)1,得到一步消去后增广矩阵:
(1)
1(2)
102
(2)A...
(2)0n
(1)
A
0
(1)
(2)(2)
步骤 4:设a220,计算第i行的消去因子:mi(2)ai(2)2/a22
(i3,...,n)
步骤 5:将增广矩阵A第i行减去mi(2)2,得到二步后增广矩阵:
(1)
1
1 (2)02(1)(2)
A 02
(3)...(3)00A
0n
(k)
步骤 6:重复上述步骤4~5,设akk(A的元素从k,k开始计数)0,(k)(k)
/akk计算第i行的消去因子:mi(k)aik
(2)(2)
(k)
(ik1,...,n)
(k)
并将增广矩阵A第i行减去mi(k)k,得到k步消去后增广矩阵:
(k)
(1)
1 * *********** a110(2)
02 ******* (1)(3)
A 00 ****
0...
(k)
000...A 0
(1)
a12(2)a22
00
(1)
x...a1(1)b1n1(2)(2)x...a2bn22...
......
...
(n)(n)x0annnbn
步骤 7:重复上述步骤,最终得到上三角矩阵。
步骤 8:从下至上,以此求解未知数xn,xn-1,...,x1,即得到原方程的解。
3.3 矩阵分解
矩阵分解是矩阵分析的重要手段。通过将一般的矩阵转化为三角矩阵、对角矩阵、正交矩阵等特殊形式的矩阵及其组合。可以从中发现矩阵的特点,为矩阵分析和计算提供便利。矩阵分析中常见的分解方法包括:针对实数矩阵的LU分解、特征值分解、针对复数矩阵的谱分解、奇异值分解以及针对正定矩阵的Cholesky分解等。本节主要介绍LU分解和Cholesky分解,其它分解可参考矩阵分析相关书籍。
3.3.1 LU分解法
三角矩阵是指矩阵对角线一侧元素全部为0的矩阵,包括上三角(右三角)矩阵U和下三角(左三角)矩阵L。
a**a00
U0...*,L*...0
**b 00b
三角矩阵具有如下性质:
(3.7)
1. 上(下)三角方阵的行列式的值等于对角线元素的乘积:uii;
i1
n
2. 上(下)三角方阵的转置为下(上)三角矩阵;
3 上(下)三角方阵的逆矩阵为(上)三角矩阵,且对角元是原三角矩阵对角元的倒数;
4. 两个上(下)三角方阵的乘积也是上(下)三角矩阵,且对角元是原三角矩阵对角元的乘积。
上述性质可以根据线性代数的知识加以推导,读者可行证明。另外,利用matlab也可以初步验证上述结论的正确性。
例3.9-1,要验证性质4,可采用matlab的符号计算功能,首先输入上三角矩阵
将两矩阵相乘,得到:
可以很容易看出,对于任意3×3矩阵,性质4是成立的。
另外,当符号运算难度较大时,可以通过随机矩阵的方法初步验证某假设的正确性。
例3.9-2,验证性质3,首先产生随机矩阵,
然后计算该矩阵的逆,得到:
可以看出,该矩阵仍然是上三角矩阵,重复上述步骤多次,如果每次都能得到相同的结论,则说明该结论“很可能”是正确的。
利用matlab验证一个假设方便快捷,虽然不能从严格意义上证明该假设的正确性,确可以快速发现某假设的错误,对于提高科学研究的效率具有重要意义。实际上,数学的发展也是首选提出一种假设或方法,而后随着认识的逐步深入,才一步步对其进行严格的证明。
方阵A的LU分解由如下定理给出:
定理 3.4 设A(aij)nn唯一存在矩阵L,U,使得ALU的充分必要条件是
det(Ai)0,i1,2,
1l121
其中,Ll31l32
ln1ln2
,n1.
u11u12
u22,U
1
u1n
u2n unn
理解:该定义可以通过高斯消去法的角度加以理解
由于:AA;
通过n次初等行变换,可以得到:
Pn...P2P1APn...P2P1A
由于Pi都是下三角矩阵(根据高斯消去法的过程得到,读者可自行分析),因此,
1Pn...P2P1为下三角矩阵,记作L,
另外,根据高斯消去法的原理,Pn...P2P1A为上三角矩阵,记作U。 因此有:
LUL1AALU。
如果已知矩阵A的LU分解,可很容易计算其对应线性方程组的解。 例3.10:已知
123100123123x114252210014,求解方程:252x18。 23153510024315x320解:
首先,AxbLUxb
其次,令yUx,构造过渡方程组:Lyb,得到:
Ly(14,18,20)T
y(14,10,72),
第三,求解方程组:Uxy
T
Ux(14,10,72)Tx(1,2,3).
T
给定矩阵A,可以通过如下公式计算该矩阵的LU分解(该公式可通过matlab符号计算功能得到)。
a11a21a31a411l21l31l41
aaaa
12223242
aaaa
13233343
aaaa
14
243444
01
001
ll
3242
l
43
0u11000010
12
uu
u11
l21u11l31u11
l41u11uluululululu
21
12
31
12
32
41
12
42
222222
uu
0u
00u
uluululuulululu
12
13
14
22
23
44
33
344413
21
13
23
31
13
32
23
41
13
42
23
43
uuu
(3.8)
3333
uluululuulululuu
14
21
14
24
31
14
32
24
34
41
14
42
24
43
34
44
观察上述公式可以发现:
首先,比较矩阵第一行可得:u11a11,u12a12,u13a13,u14a14; 其次,得到u11后,比较矩阵第一列,可得得到:l21,l31,l41; 第三,得到u12,u13,u14和l21后,可得到u22,u23,u24; 第四,依次迭代,可求解L、U矩阵的所有分量。 例3.11:求下矩阵的LU分解
43A 21
解:首先得到:
u114u3 12
然后计算:
l21
a212
0.5 u114
其次计算:
u22a22l21u1211.50.5
最后得到A矩阵的LU分解为:
10L
0.51
34
U
00.5
例3.12:求下矩阵的LU分解
03B 21
解:首先得到:
u110u3 12
但在求解l21时发现,由于u110,因此,无法计算l21。
因此,不是所有方阵都可写成ALU的形式,更一般的形式是:PALU,其中,P为排列矩阵,对A进行换行变换。Matlab中提供了LU分解的函数“lu”,该函数的调用方式为“[L,U,P] = lu(A)”,读者可利用matlab的帮助命令,查看lu函数的帮助文档,说明其中L, U, P 各是什么矩阵,以及如何构成了LU分解。
3.3.2 Cholesky分解
正定矩阵是指对任意非零向量x都满足xTAx0的矩阵。工程中许多矩阵(如,随机变量的协方差矩阵等)都表现为正定矩阵的形式。因此,研究正定矩阵的分解和求解方法对于提高算法效率具有重要意义。 正定矩阵的包括如下判定性质:
A1亦为对称矩阵,且aij0; A的顺序主子式Ak亦对称正定; A的特征值i0; det(Ak)0;
另外,若A为n阶满秩方阵,则有ATA为对称正定矩阵: 理解:
一方面,AAAAATA,因此,其为对称矩阵;
T
T
T
另一方面,根据对称矩阵定义:xTATAx(Ax)T(Ax)。由于A为满秩矩阵,
Ax0无非零解,故,Ax0,(Ax)T(Ax)0。
如果矩阵A为对称正定矩阵,则其LU分解形式更为简单。
Cholesky分解:若A为正定矩阵,则有:ACTC。 理解:由于det(Ak)0,则A能够LU分解,即ALU;
令,D0
0
00
0
0 则:ALU
LDD1U
0BLD*
**
*CD1U0
00
0
0
*
* 由于A对称,则BCCTBTCB1CTBTC(BT)1B1CT
由于C和BT均为上三角矩阵,且对角线元素互为倒数,故C(BT)1为上三角矩阵,且对角线元素为1。
同理可得:B1CT为下三角矩阵,且对角线元素为1。
由上三角矩阵等于下三角矩阵,且对角元素为1,可知C(BT)1B1CTI,其中I为单位阵。
因此有:BCTACTC。
上述证明过程是构造式的证明方法,给出了LU分解和Cholesky矩阵分解的关系: 如果已知LU分解,则cholesky分解可计算为:
CD1U
如果已知cholesky分解,则LU分解可计算为:
UCD,LCTD1
其中,D为U对角线元素开方,或C的对角线元素。
获得矩阵的Cholesky分解以后,可以利用例3.12的方法求解其对应的线性方程组的解。
Matlab中提供了Cholesky分解函数“chol”。下例给出了LU分解和Cholesky分解关系的matlab代码。读者可对利用chol函数计算的Cholesky分解结果和利用LU分解计算的结果进行比较。
3.4 范数与病态
除了直接法以外,还可以利用迭代法求解线性方程组,为了描述和分析迭代法的收敛过程,首先需要介绍向量范数和矩阵范数的概念。
3.4.1 向量范数
范数是定义在向量空间上的实值函数,其公理化的定义如下: 定义3.1 设x为
n
上的函数,且满足如下性质:
正定型:对xV,x0且x0x0;
齐次性:对kF,xV,有kxkx; 三角不等式:对x,yV,有xyxy 则称x是向量x的范数。 常用的向量范数包括: (1)最大范数或l范数:x
maxxj;
1jn
(2)和范数或l1范数:x1xj;
j1
2
(3)欧几里德范数或l2范数:x2xj;
j1
n
n
(4)Holder范数或lp范数:x
p
pxjj1
n
,p1;
容易看出,l1范数、l2范数和l范数为lp范数中取p1、p2和p的特例。 例3.13:判断下面函数是否为向量范数。
x
n
1x000
, xxi
i10x0
n
解: x范数。
xi x01,x0;000 axax,不满足性质3,因此不是
i1
图3-2 范数曲线
除了x0对应0向量以外,其他任意满足xa,a0的向量都有无穷多个,所有向量构成了高维空间中的连续曲面,如图3-2所示。对于l2
范数,满足x2a
的向量为n维空间中的球面。对于l2范数,满足x1a的向量为n维空间中的立方体。读者可以利用下面matlab程序绘制不同lp范数的范数曲线,观察其特点。
通过向量范数的定义,使得不同向量间存在了“大小”关系。由于向量包含多个元素,不同向量各元素间大小不同,因此,我们不能直接说:x1x2;但是,通过定义范数,我们可以说:x1x2。
范数的另一个重要用途是定义向量序列的极限,而向量极限的概念是迭代法的基础。
向量极限定义(3.2):给定向量序列{x0,x1,...,xk,...}
limxkx
k
ξ-δ语言描述:0,K,使得∀k>K时,满足:xkx。 值得注意的是,上述定义中并没有明确说明其范数的具体形式。实际上,对于向量范数,具有如下等价性定理:
向量范数的等价性定理(3.5):设xs和 xt是 Rn上的任意两种向量范数,则存在常数 c1,c20, 使得c1xsxtc2xs 该定理具有如下重要推论:
推论:向量序列在某范数下收敛,则在任意范数下收敛。 理解:xn在s收敛:
limxnx0limc1xnxlimc2xnx0
n
s
n
s
n
s
由向量范数等价定理,可得:
c1xnxxnxc2xnx
s
t
s
夹逼定理:已知 gkfkhk,如果limgklimhkA, 则limfkA
k
k
k
所以:limxnx=0 收敛到同一值。
n
t
3.4.2 矩阵范数
与向量类似,矩阵也可定义范数。 定义3.3 设A为
mn
上的函数,且满足如下性质:
(1)A0,且A0A0; (2)对k,kAkA; (3)ABAB; 则称为Amn的矩阵范数。
我们知道mn矩阵可以看成是一个mn维向量,因此可以按向量范数的办法来定义矩阵范数。例如,Frobenius范数(简称F范数)的定义AF为:
AF(aij)
i,j1
n
2其中,A*为A的共轭转置。
除此以外,还有另外一种更为常用的矩阵范数定义方法—算子范数。 算子范数的定义3.4:
A
max(
x0
Ax
)x
(3.9)
需要注意的是,上式中,等式左端的范数A为定义在矩阵上的矩阵范数,而等式右边的Ax和x为定义在向量的上向量范数。由于该矩阵范数由向量范数生成,因此算子范数也称为诱导范数。 常见的算子范数包括: 列范数:|A||1max(
x0
Axx1AxxAxx
2
)max|aij|
1jn
n
i1
行范数:||A||max(
x0
)max|aij|
1i
n
n
j1
谱范数:||A||2max(
x0
)以上定理的证明从略。
例3.14:求下矩阵的算子1、2和范数。
123A654
789
解:算子1范数为列和范数,将各列元素相加,得到其和分别为:14,15和13,因此,有A16。
算子范数为行和范数,将各行元素相加,得到其和分别为:6,15和24,因此,
有A24。
算子2范数为ATA的最大特征值的平方根,A216.7197。 Matlab提供了计算向量和矩阵范数的函数“norm”,例如,输入norm(A, inf),可计算矩阵A的范数。
“算子范数”名称来源于矩阵在系统分析中的重要作用。在系统分析与设计中,人们需要考察当信号通过系统后所产生的变化。
图3-3 线性系统模型
而任意离散有限线性算子可表示为矩阵形式。 例3.15,线性放大器可表示为:
a0
y(n)ax(n)
00000a00
xy
0...0
00a
差分运算可表示为(此处不考虑边界问题):
11
01
y(n)x(n1)
x(n)00
0000
0x(1)x(2)x(1)
...00x(2)x(3)x(2)
...10......
011x(n1)x(n)x(n1)001x(n)0x(n)0
累加算子可表示为:
11
n
y(n)x(i)1
i1
11
x(1)
0x(1)x(1)x(2)
0x(2)...
0...n1x(i)0x(n1)i1
n
1x(n)x(i)i1
1...01...0
1111
1
1
时移算子可表示为:
01
y(n)x(n1)0
00
00x(1)0
0...00x(2)x(1)
1...00......
0...00x(n1)x(n2)0010x(n)x(n1)0
一般地,常系数线性微分方程可表示为:
a0x(n)a1x(n1)a2x(n2)b0y(n)b1y(n1)b2y(n2)
(1)(2)(1)(2)
可记为:aIaDaDxbIbDbD121200y
其中,D(1),D(2)为一阶,二阶时移算子。
除了线性时不变系统外,其他的线性变换也可表示为矩阵形式。 例如,离散傅立叶变换可表示为:
22
00j01j02j2N
eNeNe
22j20j11j12
eNeNeNN1j2ki
y(k)eNx(i).........i0
.........2
22j(N1)1j(N1)2jN(N1)0
eNeNe
x(0)y(0)x(1)y(1)
......
x(N1)y(N2)x(N1)y(N1)
...............
2
j(N1)
eN
...
...
2
j(N1)(N1) eNe
j
2
0(N1)N
正是由于矩阵在线性系统表征的重要性,矩阵被称为“线性算子”,在下一节中我们将看到,算子范数在分析线性系统稳定性方面的重要意义,这也是“算子范
数”名称的由来。
与向量范数类似,矩阵范数在收敛的意义下是等价的。 除了等价性外,矩阵范数还有另外一个重要的性质—相容性。 定义3.5 设
与
nnnnn
分别是及的矩阵范数和向量范数。若对及AFFFa
xFn,都有AxaAxa,则称矩阵范数与向量范数是相容的。
类似地,可以定义矩阵范数间的相容性: 定义3.6 设数是相容的。
利用相容性,可以得到算子范数的三个重要性质: 算子范数与其对应的向量范数相容 根据算子范数的定义可知:
Ax
AmaxAxAx
x
是矩阵范数,若对A,B
nn
,都有ABAB,则称矩阵范
算子范数是相容矩阵范数
ABmax
ABxBx
AmaxAB xx
1in
对于任意算子范数有(A)||A||,其中,(A)max|i|称为矩阵A的谱半
径。
证明:由算子范数的相容性可得:
|| Ax||||A||||x||
将任意特征值对应的特征向量u代入得到:
||||u||||u||||Au||||A||||u||A
由于为任意特征值,则:
(A)max|i|||A||
1in
最后,给出各种范数之间的关系,如图3-4所示。
1
2,...,
1
2,...,
||A||F
图3-4 范数关系图
3.4.3 病态与条件数
除了定义向量的极限以外,范数的另一个重要应用是用于分析线性系统的稳定性。 对于线性系统,当其输入信号x(t)时,得到输出信号y(t),记作:
y(t)Lx(t)
(3.10)
但在实际中,由于误差的存在,上述过程应表示为如下形式:
y(t)eo(t)Lx(t)ei(t)
(3.11)
其中,ei(t)表示输入端的噪声,eo(t)表示噪声通过线性系统后的情况。如果输入噪声能量大于输出噪声能量,则说明该线性系统能够抑制噪声,具有较好的稳定性;反之,则称系统稳定性较差。
数学上,稳定系统和不稳定系统也称为“良态”问题和“病态”问题。 定义3.7:若原始数据有很小的变化δx,对应的输出变化δy也很小,则称该数学问题是良态问题;反之,若δy很大,则称为病态问题。 例3.16:线性方程组:
2x13x28
2x13.00001x28.00002
的解为x1=,1x2=2
若方程系数有一个小的扰动,
2x13x28
2x12.99999x28.00003
解此方程得
x1=8.5,x2=-3
可以看出,即使方程系数变化很小,方程的解也可能产生很大的变化。从系统角度分析,系统输入很小误差,对结果产生很大影响。
需要注意的是,上述两个方程的解都是其对应问题的真实解。也就是说,问题病态与否,完全取决于该数学问题本身的属性,在采用数值方法求解之前就存在,与数值方法无关。
下面讨论线性方程组微小扰动对方程组解的影响。 情况一:b存在扰动:
问题:给定方程组Axb,解为x*;另给定包含误差方程组Axbe,其解为x',分析x*和x'的差异。
解:xA1bx'A1bA1e
xx'x
A1ex
A1eb/A
A1A
e
b
公式中,
e
反映了常数向量与施加于其上的加性噪声的相对关系,A1A则b
反映了线性系统对该相对误差的抑制或增强能力。 情况二:A存在扰动:
问题:给定方程组Axb,其解为x*,另给定包含误差方程组(AE)xb(AE可逆),其解为x',分析x*和x'的差异。 解:令x'x*x,
将方程(AE)x'b与Ax*b做差,并展开,得到:
(AE)xExx(AE)1Ex(IA1E1)A1Ex 根据相容性,可得:
||x||
(IA1E)1A1E||x||
||A1||||E||
1||A1E||E||A||
E1
1||A||||A||
||A||||A1||||A||
当误差较小时,||A1||||A||
E
0,上述公式可近似为: ||A||
E||x||
||A1||||A||
||x||||A||
E同样地,公式中反映了系数矩阵与施加于其上的加性噪声的相对关系,
||A||
A1A则反映了线性系统对该相对误差的抑制或增强能力。
上述推导过程中用到了如下性质: 性质 5:(IB)1证明:
I(IB)(IB)1(IB)1BIB)1(IB)1I
1
1
1B
B(IB)1
(IB)11B(IB)1(IB)
(1B)111B
(IB)1
情况三:A,b都存在扰动:
问题:给定方程组Axb,其解为x*,另给定包含误差方程组(AE)xbe,其
解为x',分析x*和x'的差异 解:不加推导地给出其近似公式:
||x||
||x||
eEeE)||A1||||A||()
A||b||||A||||b||||A||
1||A1||||A||
||A||
(
||A1||||A||
通过上述三种情况的分析可以发现,线性系统的稳定性由||A1||||A||决定,称为(该(相容)矩阵范数对应的)条件数。 常用的条件数包括: (1)(2
)
cond(A,)AA1
cond(A, 2)A2A
1
2
条件数具有如下性质:
对任何非奇异矩阵A,cond(A)1
对非奇异矩阵A和常数c0,有cond(cA)cond(A) 对正交矩阵A,2-范数对应的条件数cond(A,2)1 当条件数很大时,方程组Axb是病态问题; 当条件数较小时,方程组Axb是良态问题。
需要注意的是,条件数是系数矩阵的特征,与算法的选择无关。另一方面,条件数与所选择的范数有关,不同范数计算的条件数不同。 典型的病态矩阵—Hilbert矩阵,其定义如下:
11Hn2
1n
1
2131n1
1n
1...
n1
11
...n22n1
...1314
利用matlab函数“hilb”,产生3阶、5阶、7阶… … Hilbert 矩阵,用matlab函数“cond”计算相应的条件数,编程分析矩阵误差对结果的影响。
3.5 迭代法
3.5.1 迭代法的基本概念
本章第二节介绍的解线性方程组的直接法利用有限次操作得到问题解。迭代法是从初始猜测值开始,通过依次逼近获得问题解的方法。迭代法是非线性方程组求解和最优化问题的基本方法,数值计算中的Jacobi法、Gauss-Seidel法,最优化中的Newton迭代法、最快下降法(Steepest Descent)、数字信号处理:维纳滤波、卡拉曼滤波等都用到了迭代法。 迭代序列的基本过程为: 给定函数f():
n
n
和迭代初值,迭代过程描述为:x0a
x(1)f(x(0))x(2)f(x(1))......
x(k1)f(x(k))......
迭代过程具有如下性质:
迭代法原理:如果迭代序列收敛,则其极限点为方程f(x)x的解。 证明:由于{x(k)}收敛,则x(k1)f(x(k))有:
0, K,使得∀k>K时,满足:x(k)x*
2
x(k1)x(k)(x(k1)x*)(x(k)x*) x
(k1)
x+x
*(k)
x
*
即:x(k1)-x(k)收敛,
根据迭代公式,x(k1)f(x(k)),则:
f(x(k))-x(k)
(k)(k)
即limf(x)-x0 k
因此,极限点满足f(x)x*0。
其中,x*称为函数f(x)的不动点(fixed point)。
在上述推导过程中,对函数f(x)并没有特别的要求。因此,迭代法可用于各类方程(组)(线性方程(组)、非线性方程(组))求解。这也是迭代法较之直接法最大的优势之一。
迭代法可以利用简单的规则,通过多次迭代得到时分复杂的结果,分形结构是对迭代法最直观的展现。
分形是一种粗糙的或分裂的几何形状,其每一部分都是整体的缩小比例的复制(自相似性)。下图展示了分形结构的生成过程,从简单的图形出发,通过简单的迭代规则,最终获得十分复杂的图形。
图3-5 分形结构迭代过程
利用分形可以绘制十分精美的图案,下面给出几张利用迭代方法得到的分形图。
图3-6 分形图像
另外,分形技术还应用于三维虚拟现实以及超宽带天线设计等领域。NOVA纪录片“寻找隐藏的维度”(hunting the hidden dimension)对分形技术的应用进行了详细地介绍,读者有兴趣可以了解更多的相关内容。
3.5.2 迭代法的构造
根据迭代法的原理,如果迭代序列x(k1)f(x(k))收敛,则其极限点为方程
f(x)x的解。利用迭代法求解首先要将原方程等价转化为f(x)x的形式。对
于同一方程,上述转化过程并不唯一。 例3.17:用迭代法求解方程2x1
解:首先构造与之等价的方程,很容易看出,下列方程都与原方程等价。
2x1x1x
1x
x
3
x2(11.5x)
110x
x
12
利用上述方程可分别构造对应的迭代公式:
x(k1)1x(k)
2(x)(1x)/3x(k1)(1x(k))/3
3(x)2(11.5x)x(k1)2(11.5x(k))4(x)(110x)/12x(k1)(110x(k))/12
各迭代序列的收敛过程如图3-7所示。
8
1(x)1x
f(x) = 2 * (1 - 1.5 x)
f(x) = (1 + 10x) / 12
图3-7 不同迭代过程
各种迭代方法前四步的迭代结果如下表所示:
x0 x1 x2 x3 x4
1-x 0 1 0 1 0
(1+x)/3 0 1/3 0.4444 0.4815 0.4938
2(1-1.5x) 0
2 -4 14 -40
(1+10x)/12 0 0.0833 0.1528 0.2106 0.2589
… … … … … … … … … …
可以看出,对于不同的迭代公式,有的序列收敛,有的序列发散,有的序列振荡。对于收敛序列,不同序列的逼近真值0.5的速度也不同,对于序列2(x),经过4次迭代,其值已经很接近0.5,而对于迭代序列4(x),经过4次迭代,其值距离0.5仍相差很大。
从上例可以看出,迭代法的关键是如何构造迭代函数f(x)x,使得迭代序列收敛,且快速收敛。
3.5.3 Jacobi迭代
对于线性方程组Axb,其迭代公式的一般形式可表示为:
xMxg
其中,M称为迭代矩阵。
同样,对于同一线性方程组Axb,其对应的构造方法多种多样。 例3.18:
AxbAxb0 (AI)xbx
其对应的迭代公式为:x(k1)(AI)x(k)b 另外,
A1AxA1bx0xAb
1
其中,迭代矩阵M等于0。当然,本迭代矩阵并没有实际应用价值,但其说明一点:虽然迭代方法种类繁多,但构造迭代矩阵的开销不能大于矩阵求逆的开销。否则该方法在运算量和计算精度方面较矩阵求逆的方法没有优势。 下面介绍一种常用的求解线性方程组的迭代公式—Jacobi迭代公式。
设Axb,并将A矩阵划分为三部分:ALDUD(LU)(注意:此处的L,U和D矩阵与矩阵分解中介绍的L、U和D矩阵有着完全不同的含义!)
0a021
其中,La31a32
an1an2
an,n1
a11
a22
,D
0
a13a230
...a1n
...a2n
......
...0
, ann
0a120U
于是原方程组Axb变成:[D(LU)]xb,即:
Dx(LU)xb
由于det(D)0,所以有xD1(LU)xD1b. 令MJD1(LU),gJD1b,得到如下迭代公式:
x(k1)MJx(k)gJ11MJD(LU),gJDb
上式称为Jacobi迭代式(矩阵形式)。
(3.12)
由于矩阵D为对角矩阵,D1可以很方便的计算得到。另外,在实际编程中,上述过程一般采用分量的形式实现,即:
x1(k1)(b10a12x2(k)a13x3(k),,a1nxn(k))11,
x2(k1)(b2a21x1(k)0a23x3(k),,a2nxn(k))22,
(k1)(k)(k)(k)(k)x(baxaxax,,ax0)nn.nnn11n22n33n,n1n1
因此,在Jacobi算法求解过程中,并不需要复杂的矩阵求逆操作。
(3.13)
例3.19 用Jacobi迭代法解线性方程组Axb,其中
54TA,b(1,1). 45
解: (1)建立迭代公式
x(k1)MJx(k)gJ
00.81T
其中,MJD1(LU) gDb(0.2,0.2),J
0.80
各次迭代的结果如下所示。
迭代次数 0 1 2 3 4 5 6 7 „„
x的第一分量 0.0 0.6 0.36 0.744 0.5904 0.8362 0.7379 0.8951 „„
x的第二分量 0.5 0.2 0.68 0.488 0.7952 0.6723 0.8689 0.7903 „„
原方程的精确解x*(1,1)T,由表可以看出,迭代向量x(k)逐步逼近x*,当迭代次数足够多时,迭代序列收敛到方程的精确解。
3.5.4 高斯-赛德尔迭代法
一般认为新近似解要比老近似解更接近真实解,将已计算出的x(k1)分量替换Jacobi迭代公式中x(k)相应分量即可得到Gauss-Seidel迭代。 对Jacobi迭代公式的分量形式作如下变形:
xi
(k1)
(biaijxj(k))ii
j1ji
n
(biaijxj(k)
j1
n1
ji1
ax
ij
n
(k)j
)ii,i1,2,
,n.
其中,若将x(k1)的分量逐个替换x(k)的相应分量,使得
xi
(k1)
(biaijxj
j1
i1
k(
1)
ji1
ax
ij
n
kj))ii,i1,2,
,n. (3.14)
称上式为Gauss-Seidel迭代式的分量形式。 下面利用例子说明Gauss-Seidel迭代法的过程。
例3.20 用Gauss-Seidel迭代法解线性方程组Axb,其中
54T
A,b(1,1).
45
解: (1)建立Jacobi迭代公式
x(k1)MJx(k)gJ
00.8T
g(0.2,0.2)其中,MJ ,J
0.80
选择初值为:x(0)[0,0.5]。 一次迭代:
首先,利用Jacobi迭代公式计算x(1)的第一个分量:
0
x1(1)00.80.20.6
0.5
此时,如果继续利用Jacobi迭代公式计算x(1)的第二个分量:
0(1)
x20.800.20.2
0.5
则得到了Jacobi迭代的一次迭代结果,如表?所示。
但是与Jacobi迭代不同的是,在计算x(1)的第二个分量时,Gauss-Seidel利用刚刚计算得到的x1(1)替换x1(0),从而得到如下公式:
0.6(1)
x20.800.20.68
0.5
按照上述步骤得到的x(1)0.60.68即为Gauss-Seidel迭代法的结果。 二次迭代:
首先,利用Jacobi迭代公式和第一步得到的Gauss-Seidel迭代结果,
x(1)0.60.68计算x(2)的第一个分量:
06
x1(2)00.80.20.744
0.68
然后,利用本步得到的x1(2)的第一个分量0.744,计算x(2)的第二个分量:
0.744(2)
x20.800.20.7952
0.68
重复上述过程,可以得到Gauss-Seidel的各步骤结果:
迭代次数 0 1 2 3 4 5 6 7 „„
x的第一分量 0.0 0.6 0.744 0.8362 0.8951 0.9328 0.9570 0.9725 „„
x的第二分量 0.5 0.68 0.795 0.8689 0.9161 0.9463 0.9656 0.9780 „„
对比Jacobi迭代的迭代过程可以发现,Gauss-Seidel经过7次迭代,与真实值[1, 1]的误差要小于Jacobi迭代误差。
另外,从上例的计算过程可以看出,Gauss-Seidel迭代中没计算一个分量,都需要本次迭代中前面分量的信息。因此,Gauss-Seidel迭代较Jacobi迭代更不适于并行实现,但另一方面,由于新值的使用,Gauss-Seidel迭代可以用一组内存空间存储迭代向量,相对Jacobi迭代更节省内存空间。
Gauss-Seidel的迭代过程也可以表示为矩阵形式。 给定Jacobi迭代公式:
x(k1) bx(k)bx(k)bx(k)g1221331nn11
(k1)(k1)(k)(k)x2b21x1 b23x3b2nxng2
(k1)(k1)(k1)(k1)
bn1x1bn2x2bn,n1xn1 gxnn
0
令:Lb21
bn1
0000
b
n2
00
0 U0
00
00
b
12
000
b
b
2,n 0
1n
则Gauss-Seidel的迭代过程可表示为:
x(k1)Lx(k1)Ux(k)g(IL)x(k1)Ux(k)g
x(k1)(IL)1Ux(k)(IL)1g
也可用矩阵A的元素来表示上述过程。 令:
00
a021
La31a32
an1an2
00a12
0 U
0
0
a13a23
a1n
a2n
an1n0
ann1
则有:LD1L和UD1U
由x(k1)(IL)1Ux(k)(IL)1g可得:
x(k1)(DL)1Ux(k)(DL)1b
(3.15)
因此,Gauss-Seidel的迭代矩阵为M(DL)1U,g(DL)1b。
3.5.5 松弛迭代法
在Gauss-Seidel迭代法的基础上,还进一步演化出了松弛迭代法(SOR迭代),又称为加速的Gauss-Seidel公式。该方法包括两部分构成: 首先:迭代步,用Gauss-Seidel公式给出预报值xi
(k1)
。
xi
(k1)
(biaijxj
j1
i1
(k1)
ji1
ax
ij
n
(k)j
)ii,i1,2,
(k1)
,n.
其次,加速步,即将当前迭代值xi(k)与预报值xi即:
xi(k1)xi
(k1)
线性组合得到校正值xi(k1),
(1)xi(k),i1,2,
,n.
其中,(0,2)称为松弛因子,当(0,1)时,称为松弛公式;当(1,2)称为超松弛公式;特别地,当1时,则为Gauss-Seidel公式。 合并以上两式得到SOR迭代的迭代公式为:
xi(k1)xi(k)(xi
xi(k)[xixi
(k)
(k1)
xi(k))aii(k)
xi]aij
i1
(k1)
j1
(k1)
aijxj(k)),i1,2,
j1n
aii
(biaijxj
,n.
或者写成:
xi(k1)xi(k)xi(k),xi
(k)
aii
(biaijxj
j1
i1
(k1)
aijxj(k)),i1,2,
j1
n
,n
利用Gauss-Seidel公式可以得到松弛迭代法的矩阵表示:
x(k1)D1Lx(k1)D1Ux(k)D1b
x = x(k1)x(k)=D1Lx(k1)(D1U)x(k)D1bx(k)x(k1)x(k)x
x(k)(D1Lx(k1)D1Ux(k)D1bx(k)) (1)x(k)D1Lx(k1)D1Ux(k)D1b
-1
(I-D-1L)由于I-D-1L1,故和(D-L)-1存在,有:
x(k1)(DL)1[(1)DU]x(k)(DL)1b 因此,SOR法的迭代公式为:
M(DL)1[(1)DU], g (DL)1b
(3-16)
下列代码为Jacobi、Gauss-Seidel和SOR法迭代过程的演示程序,其结果如图3-8
所示,红色为Jacobi法结果,绿色为Gauss-Seidel法结果,蓝色为SOR法结果,松弛系数0.4。读者可自己运行本代码,观察不同迭代法的收敛过程以及SOR法松弛因子对迭代过程的影响。
图3-8 三种迭代方法的迭代过程
3.5.6 迭代法的收敛性
前面的讨论可以看出,对于同一方程,可以构造不同的迭代序列,但不同的迭代序列是否收敛,以及收敛速度各不相同。本节将讨论不同迭代方法的收敛性能。 我们首先观测迭代序列与线性方程组解的真实值之间的差异,可以很容易得到如下性质:
性质:如果线性方程组的精确解x*存在,即x*满足:x*Mx*g,则迭代序列与的x*差满足如下关系:x(k)x*Mk(x(0)x*)。 该性质可以很容易的通过迭代推导获得:
x(k)x*Mx(k1)gMx*g
M
(x(k1)x*)M2(x(k2)x*)
Mk(x(0)x*)
假设迭代矩阵的特征值分解存在,即M可以写成MQ1Q的形式,则迭代序列可进一步写为:
x(k)x*Mk(x(0)x*)
=Q1QQ1Q...Q1Q(x(0)x*) =Q1kQ(x(0)x*)
其中,为对角矩阵,
1
2
1k
k,...
n
2k
...
nk
k0,明显地,当M矩阵的所有特征值的绝对值都小于1时,随着迭代的继续,
相应地,迭代结果与精确解的误差x(k)x*也趋近于0。
另一方面,当存在M矩阵的某个特征值的绝对值大于1时,k,迭代结果与精确解的误差x(k)x*也趋近于。
由于并非所有矩阵都能够进行特征值分解,因此,上述分析过程不能涵盖所有情况。幸运的是,通过进一步严密分析,上述结论对任意线性方程组迭代都成立,因而得到如下定理。
定理3.6 设x(k1)Mx(k)g,其中,M为迭代矩阵,对任意x(0)Rn,迭代公式收敛的充分必要条件是(M)1。
根据定理3.6可以看出:迭代法收敛与否只决定于迭代矩阵的谱半径,与初始向量和右端项无关。该定理的证明是如下引理的直接结果。 引理:设A为n阶方阵,则limAk0的充分必要条件为(A)1。
k
证明:
必要性:若limAk0
k
由矩阵收敛定义可知:limAk0 ,又因ρ(A)≤A,故有:
k
0(Ak)Ak
由夹逼定理,则有:lim[(Ak)]=0 。
k
又由于 (Ak)[(A)]k 则:lim[(A)k]=0(A)1
k
充分性:若(A)1,取存在矩阵范数
,使得:
1-(A)
0, 2
A(A)
k
k
1(A)
1 2
则有:limA0
由算子范数的相容性,可得:
0AkA
k-1
A...A
k
k
由夹逼定理可得:limAk0。
上述证明过程中用到了如下定理。
定理3.7:设A为任意n阶方阵,则对于任意正数,存在矩阵范数
,使得:
A(A)
该定理的证明超出本书范围,有兴趣的读者可寻找相关书籍,自行证明。 另外,利用矩阵范数的性质:(A)A,可得到关于线性迭代收敛性的如下推论。
推论1:对任意初始向量x(0)和右端项g,若M1,由迭代公式
Mx(k1)Mxk()g产生的向量序列x(k)收敛。
需要注意的是,推论1为充分条件,而非必要条件,即:
如果M1条件满足,则可以判定该迭代过程收敛;但反之,如果M1,则不能认为迭代过程发散,需要利用谱半径,做进一步的判断。 例3.21:对于迭代矩阵
0.8968 0.0337 -0.0678
M 0.0337 0.5434 0.7166
-0.0678 0.7166 -0.5402
其谱半径为(M)0.9,其范数为:M断该迭代矩阵是否收敛。
例22:讨论下述方程组Jacobi和Gauss-Seidel方法的收敛性。
1.3247。利用无穷范数,无法判
x12x2-2x31
x1x2x32 2x2xx3123
解:首先写出Jacobi和Gauss-Seidel方法的迭代矩阵:
MJa
0-22I-D-1A-10-1
-2-20
2-2
其对应的特征方程为:I-B1130
2
2
因此有1230,故(B)01。 因此Jacobi迭代法收敛。
对于Gauss-Seidel方法,其迭代矩阵为:
0-22
M(D-L)-1U02-3
002
2
其对应的特征方程为: I-M0-2
-2
3(-2)20
-2
因此有10,232,故(M)21。 因此Gauss-Seidel迭代法发散。
从本例中可以看出,虽然在引入Gauss-Seidel迭代法时我们提到:一般认为新近似解要比老近似解更接近真实解,但在实际中,也存在Jacobi迭代法收敛而Gauss-Seidel迭代法不收敛的情况。
因此,迭代法的收敛性必须通过分析其谱半径获得。 对于松弛法,其收敛性满足如下必要条件: 推论2:松弛法收敛的必要条件是02。 证明:设松弛法迭代矩阵M有特征值1,2,...,n 因为:det(M)12...n[(M)]n
有收敛性定理,松弛法收敛必有:det(M)1 又因为:
det(M)(D-L)-1(1-)DU(D-L)-1
1a11a22...ann
(1-)DU(1-)na11a22...ann
因此有:det(M)(1-)n1,即:02。
该推论给出了松弛法中松弛因子的选择范围:02。
除了收敛定理外,对于特殊的方程,还存在一些其他的收敛性判别方法: 定义3.8 (强超条件与弱超条件)若A(aij)nn满足aiiaij,
j1
jin
i1,2,,n.,
则称A具有严格对角优势,或者说A满足主对角元强超条件(简称强超条件); 若A满足aiiaij,
j1
jin
i1,2,,且至少有一个i值满足严格不等式,,n(.4-22)
则称A具有对角优势,或者说A满足主对角元弱超条件(简称弱超条件)。
832
例3.23 A4111满足式(4-21),称A满足强超条件;
6312421
B242满足式(4-22),称B满足弱超条件;
1230021
1210
也满足弱超条件。 同理,C
01210012
定义3.9(不可约矩阵)若A(aij)nn,通过行对换和相应的列对换变成
A11
0A12
,其中A11,A22为方阵,则称A为可约矩阵;反之为不可约矩阵。例如: A22
a13
a23为可约矩阵, a33a13a23a33a43
a14a24
为不可约矩阵, a34a44
a11a120a
220a32a11a12aa21220a32
0a42
又如,在例3.23中,A为强超不可约矩阵;B为弱超不可约矩阵。 定理3.8 (非奇异条件)设矩阵A(aij)nn. 若A满足主对角元强超条件,则det(A)0;
若A满足主对角元弱超条件且不可约,则det(A)0. 定理3.9 (收敛的充分条件)设Axb.
(1)若A满足主对角元强超条件,则对x(0)Rn,Jacobi迭代法及Gauss-Seidel迭代法均收敛;
(2)若A满足主对角元弱超条件且不可约,则对x(0)Rn,Jacobi迭代法及Gauss-Seidel迭代法均收敛;
例3.24 以下列矩阵作为线性方程组的系数矩阵,试判别它们对Jacobi迭代法及Gauss-Seidel迭代法的收敛性。
521
(1)A1132。
112
解 该矩阵弱超不可约,故据定理3.8,以它为系数矩阵构成的线性方程组,对
x(0),Jacobi迭代法及Gauss-Seidel迭代法均收敛;
122
(2)A2111。
221
解 该矩阵不满足强超或弱超条件,不能应用定理3.9,可考虑用定理3.6。对Jacobi迭代法,设A2的迭代矩阵为BJ,且
122
MJ101,
220
MJ
1,M11
22
det(IMJ)11248。
22
除了判断迭代方程收敛性以外,算子范数也可用于判断迭代方程的收敛速度,其相关定理如下:
定理3.10:设有迭代公式x(k1)Mx(k)g,若M1,x(k)收敛于x*,则有误差