数值分析 第三章 基于MATLAB的科学计算—线性方程组1

及其基于ＭＡＴＬＡＢ的程序实现与分析 三、 解线性方程组（线性矩阵方程）

解线性方程组是科学计算中最常见的问题。所说的“最常见”有两方面的含义：

１） 问题的本身是求解线性方程组；

２） 许多问题的求解需要或归结为线性方程组的求解。 关于线性方程组

AxBxA1B

其求解方法有两类：

（１）

１） 直接法：高斯消去法（Gaussian Elimination）； ２） 间接法：各种迭代法（Iteration）。

１、高斯消去法

１） 引例

考虑如下（梯形）线性方程组：

x12x2x32121x12x12x2x323.5034x1x114x1 3x4x12332232x312x30.500x31高斯消去法的求解思路：把一般的线性方程组（１）化成（上或下）梯形的形式。 ２）高斯消去法——示例

考虑如下线性方程组：

x1x2x312x12.001x22.01x35x19x22x33060111x11122.0012.01x5

2921x330601１） 第一个方程的两端乘

2加到第二个方程的两端，第一个方程的两端乘 1－１加到第三个方程的两端，得

11x11100.0010.01x3

21010x330600２） 第二个方程的两端乘10加到第三个方程的两端，得 0.00111x11100.0010.01x3

201010x360600３） 从上述方程组的第三个方程依此求解，得

x1x2x312401x2100030.01x33000 x6003３）高斯消去法的不足及其改进——高斯（全、列）主元素消去法

在上例中，由于建模、计算等原因，系数2.001而产生0.0005的误差，实际求解的方程组为

1122.000591x11x12565.22.01x5x23014.9 2x450.702x33060131注：数值稳定的算法

高斯列主元素消去法就是在消元的每一步选取（列）主元素—一列中绝对值最大的元取做主元素，高斯列主元素消去法是数值稳定的方法。

列主元素消去法的基本思想：在每轮消元之前，选列主元素（绝对值最大的元素）,使乘数lik1.

列主元素消去法的步骤：设已经完成第1步到第k1步的按列选主元、交换两行、消元计算，得到矩阵

A(k),b(k)(1)a11(1)))a12a1(1a1(1kn(2)(2)(2)a2a2a22kn(k)(k)akkakn(k)(k)ankannb1(1)(2)b2.bk(k)(k)bn

第k步计算如下：对于k1,2,,n1，

(k)aik（1） 选列主元素，即确定i0使ai(kk)max； kin0（2） 如果ai(kk)0，则方程组解不唯一，或者A接近奇异矩

0阵，停止运算；

（3） 如果i0k，则交换[A,b]第i0行与第k行元素； （4） 消元计算：

(k1)(k)(k)(k)(k)aijaijlikakj,likaik/akk, (k1)(k)(k)bilikbk,i,jk1,k2,,n.bi（5） 回代计算：

xnbn/ann,n x(bax)/a,in1,n2,,2,1.iijjiiiji1完全主元素消去法即是每次选主元时，依次按行、列选取绝对值最大的元素作为主元素，然后交换两行、两列，再进行消元计算.

完全主元素消去法的步骤：设已经完成第1步到第k1步的选主元、交换行和列、消元计算，得到矩阵

(1)a11(1)))a12a1(1a1(1b1(1)kn(2)(2)(2)(2)a2a2a22b2kn. (k)(k)akkaknbk(k)(k)(k)(k)ankannbnA(k),b(k)第k步计算选主元素的范围为ki,jn，即确定i0,j0使

(k). ai(0k,)j0maxaijki,jn第k步计算如下：对于k1,2,,n1，

(k)（1） 选主元素，即确定i0,j0使ai(k,)jmaxaij；

00ki,jn（2） 如果ai(k,)j0，则方程组解不唯一，或者A接近奇异矩

00阵，停止运算；

（3） 如果i0k，则交换[A(k),b(k)]第i0行与第k行元素；如果

j0k，则交换[A(k),b(k)]第j0列与第k列元素；

（4） 消元计算：

(k1)(k)(k)(k)(k)aijaijlikakj,likaik/akk, (k1)(k)(k)bilikbk,i,jk1,k2,,n;bi（5） 回代求解.

【注】 完全主元消去法是解低阶稠密矩阵方程组的有效方法，但完全主元消去法解方程组，在选主元素时要化费较多的计算机时间，行主元消去法与列主元消去法运算量大体相同，实际计算时，用列主元消去法即可满足一定的精度要求.

对同一数值问题，用不同的计算方法，所得结果的精度大不一样.对于一个算法来说，如果计算过程中舍入误差能得到控制，对计算结果影响较小，则称此算法是数值稳定的；否则，如果计算过程中舍入误差增长迅速，计算结果受舍入误差影响较大，则称此算法为数值不稳定的.因此，我们解数值问题时，应选择和使用数值稳定的算法，否则如果使用数值不稳定的算法，就可能导致计算失败. 例 用高斯列主元消去法解方程组

x12x23x31,2x14x25x30, 3x5x6x0.2311231356015/320r1r3解 [A,b]2450245002/310

3560123101/311101/201001013/200103. 001/210012所以，方程组的解为x(1,3,2).

T

４）高斯列主元素消去法的MATLAB实现：xA\\b，意为

xA1b．

例 LinearEquiation02.m

open LinearEquiation02 LinearEquiation02

一个典型的例子: Hilbert矩阵:Hij1

1注: 非奇异矩阵的条件数:

５）ＬＵ分解 （ＬＵ Factorization）（高斯消去法、Doolittle分解）

高斯消去法的消元过程，从代数运算的角度看就是用一个下三角矩阵L左

乘方程组的系数矩阵Ａ，且乘积的结果为上三角矩阵U，即

L1AULybAxLUxb， (２) UxyALU可通过直接用A元素计算矩阵A的三角分解矩阵L和U.这种直接计算A的三角分解的方法有实用上的好处.下面利用矩阵乘法规则来确定三角矩阵L和U.

1u11u12u1nl211u22u2nAl31l321unnll1n1n2. 第一步：利用A的第一行、第一列元素确定U的第一行、L的第一列元素.由矩阵乘法，

a1j(1,0,0,,0)(u1j,u2j,uij,0,,0)Tu1jai1(li1,li2,lii1,1,,0)(u11,0,,0)Tli1u11(j1,2,,n)， (i2,3,,n)，

得到

u1ja1j(j1,2,,n),li1ai1/u11(i2,3,,n). （3.7）

设已经计算出U的第1至r -1行元素，L的第1至r -1列元素，现在要计算U的第r行元素及L的第r列元素.

第r步：利用A的第r行、第r列剩下的元素确定U的第r行、L的第r列元素(r2,3,,n).由矩阵乘法，有

arj(lr1,lr2,,lrr1,1,0,,0)(u1j,u2j,,ujj,0,,0)T

lrkukjurjk1r1(jr,r1,,n)，

得U的第r行元素为

urjarjlrkukjk1r1(jr,r1,,n,r2,3,,n). （3.8）

由air(li1,li2,,lii1,1,0,,0)(u1r,u2r,,urr,0,,0)T

likukrlirurrk1r1(ir1,r2,,n)，

得

lir(airlikukr)/urrk1r1(r2,3,,n1,ir1,,n). （3.9）

例5 用LU分解法求解方程组

223x23477x1. 2245x37解 对系数矩阵A进行LU分解，

u112,u122,u133,l212,l311.

由u2ja2jl21u1j，有u223,u231.

l32(a32l31u12)/u222，u33(a33l31u13l32u23)6.

因此

1ALU2112122331. 6解方程组Lyb，得y13,y212y15,y37y12y26. 解方程组Uxy，得x31,x2(5x3)/32,x1(32x23x3)/22. 6） LU 分解的MATLAB实现：[L,U]lu(A)或[L,U,P]lu(A) 例

A=rand(5); [L,U，P]=lu(A)

A=rand(5); [L,U,P]=lu(A) L=P\\L

当A是主对角占优的三对角矩阵时，基于Doolittle分解可得到解这类方程组的追赶法。

２、Cholesky分解 (Cholesky Factorization)

对称正定矩阵A的Cholesky分解和以A为系数矩阵地的线性方程组的改进的平方根法：

设n阶方程组Axb，A是对称正定矩阵(Positive Definite Matrix)，则A有三角分解

1lu11u12u1n121u22u2nALUl31l321unnll1n1n2. 再将U分解为

u11u12u1nu22u2nUunnu11u22u1nu121uu1111u12nDU0, u22unn1则ALDU0.

（1） 对称正定矩阵有唯一的分解ALDLT

TDLT，且对称阵ATA，则有 这是由于ALDU0，ATU0TLDU0U0DLT

再利用三角分解的唯一性，得U0TL.因此，对称正定矩阵有唯一的分解ALDLT.

（2） D是正定对角阵（即uii0,i1,2,,n）

由于A对称正定的充要条件是BABT对称正定，其中B是n阶可逆方阵.取BLT，就推知D是正定对角阵.

因此D的对角元素uii0(i1,2,,n)，记DDD，其中

u111D2u22, unn1212则ALDDL(LD)(LD).

T12121212T（3） 乔莱斯基（Cholesky）分解

将LD记为L，则ALLT称为Cholesky分解.利用Cholesky直接分解公式，推导出的解方程组方法，称为Cholesky方法或平方根法.

（4） 解方程组的平方根法（Cholesky方法） 由Cholesky分解，有

l11l11l21ln1llll212222n2LLT. （3.10） Allllnnnnn1n212利用矩阵乘法，逐步确定L的第i行元素lij(j1,2,,i).

由aijlikljklijljj（当jk时，ljk0），有分解公式：对于

k1j1i1,2,,n

lij(aijlikljk)/ljj(j1,2,,i1),lii(aiilik)1/2 （3.11）

2k1k1j1i1将对称正定矩阵A作Cholesky分解后，则解方程组Axb就转化为解两个三角方程组Lyb,Lxy.

T例7 用Cholesky方法解方程组

114414.252.75x6. 12.753.57.25解 对系数矩阵作Cholesky分解得到

114220.50.514.252.75LLT0.52. 21.5112.753.50.51.51解Lyb，得y(2,3.5,1)T. 解LTxy，得x(1,1,1)T.

cholesky分解的MATLAB的实现：L=chol(A)。 3、追赶法

在许多实际问题中，如，常微分方程两点边值问题、三次样条插值方法等，往往遇到线性方程组Axf的求解，其中

b1c1abc222,f(f1,,fn)T. （3.13） Aabcn1n1n1anbn称具有公式（3.13）形式的系数矩阵A为三对角阵,称相应的线性方程组为三对角方程组（Tridiagonal Linear Systems).具有这种形式的方程组在实际问题中是经常遇到的，而且往往是对角占优（Diagonally Dominant）的.A满足条件：

b1c10，

biaici0,aici0(i2,3,,n1)，

bnan0.

这类方程组Axf的解存在唯一（A非奇异），可以直接利用高斯消去法或直接分解法，而其解答可以用极其简单的递推公式表示出来，即下面介绍的追赶法.追赶法通常是数值稳定的.

对A作LU分解（Doolitle分解），可以发现L、U具有非常简单的形式.

1m2ALU1m31mn1l1u1lu22. un1ln由矩阵乘积，得

l1ml21cn1bnu1m2u1l2u2mnln1un1mnun1ln. b1c1ab22c2Aan1bn1an比较等式两端，得到

i1,2,,n1,uici,lb11 （3.14） ma/l,i2,3,,n,ii1ilibimiui1, i2,3,n.因为上述分解ALU，则方程组Axf的求解转化为解两个简单的三角方程组Lyf和Uxy，从而得到求解方程组的算法公式.

先解Lyf，即

y1f1,yifimiyi1(i2,3,,n). （3.15）

再解Uxy，即

xnyn/ln,xi(yiuixi1)/li(in1,n2,,1). （3.16）

这种把三对角方程组的解用递推公式（3.14）、（3.15）、（3.16）表示出来的方法形象化地叫做追赶法，其中（3.14）、（3.15）是关于下标i由小到大的递推公式称为追的过程，而（16）却是下标i由大到小的递推公式称为赶的过程，一追一赶构成了求解

Axf的追赶法.

例9 用追赶法解三对角方程组

411141x3. 142解 系数矩阵分解得到

1L10.250.266741,U3.75. 113.7333解Lyf，得y(1,3.25,2.8667)T. 解Uxy，得x(0.5179,1.0714,0.7679)T. 调用函数LU_Factorization.m解例9.输入 A=[4 -1 0;-1 4 -1;0 -1

4];b=[1;3;2];[x,L,U,index]=LU_Factorization(A,b) 得到方程组的解及相应的LU分解矩阵：

x = 0.5179 L= 1.0000 0 0 U= 4.0000 -1.0000 0 1.0714 -0.2500 1.0000 0 0 3.7500 -1.0000 0.7679 0 -0.2667 1.0000 0 0 3.7333

为了对线性方程组的直接法作出误差分析，为了讨论方程组迭代法的收敛性，需要对向量和矩阵的大小进行度量，进而引入了

范数━

用于度量“量”的大小的概念

引言

实数的绝对值：a是数轴上的点a到原点0的距离；

复数的模：abia2b2是平面上的点a,b到原点0,0的距离； 还有其他刻画复数大小的方法（准则）：如 １）ab； ２）maxa,b

向量的内积、范数及n维空间距离的度量

令P是一数域，Pn是P上的向量空间，如果函数

x，y:PnPnP有如下性质：

１、共轭对称性：x，yPn，y，xx，y；

２、非负性：xPn，x，x0，x，x0x0； ３、线性性：x，y，zPn，a，bP，

axby，zax，zby，z；

则称x，y是Pn上的一个向量内积（inner product），向量空间

Pn上的向量内积通常用符号x，y表示，定义了内积的向量空间Pn称为内积空间（inner product space）。记做Pn，·，·表示。 例１．１ QQTPnn，Q0，x，yPn，容易验证函数

x，yxTQy （１．１）

定义了Pn上的一个内积。

令P是一数域，Pn是P上的向量空间，如果函数x:PnP有如下性质：

１、非负性：xPn，x0，x0x0； ２、齐次性：xPn，aP，axax； ３、三角不等式：x，yPn，xyxy；

则称x是Pn上的一个向量范数（norm），向量空间Pn上的范数通常用符号x表示。定义了范数的向量空间Pn称为赋范空间（normed space）。记做Pn，·表示。

例１．２ QQTPnn，Q0，xPn，容易验证函数

xxTQx （１．２）

定义了Pn上的一个范数，这样定义的范数称为由内积（１．１）诱导的范数。

例１．３ RnCn上常用的向量范数：

xx1，x2，，xnRnCn，

T１、１—范数：x1xk；

k1n２、２—范数：x2xTx； ３、—范数：xmaxx1，x2，，xn；

1kn令P是一数域，Pn是P上的向量空间，如果实值函数

dx，y:PnPnR有如下性质：

１、对称性：x，yPn，dy，xdx，y；

２、非负性：x，yPn，dx，y0，dx，y0xy

３、三角不等式：x，y，zPn，dx，zdx，ydy，z； 则称dx，y是Pn上的一个距离（函数）（distance function）或度量（metric），定义了度量的向量空间Pn称为度量空间（metric space），记做Pn，d表示。

例１．4 RnCn上常用的（由范数诱导的）度量：x，yRnCn， １、１—范数诱导的度量：dx，yxy1； ２、２—范数诱导的度量：dx，yxy2； ３、—范数诱导的度量：dx，yxy； 矩阵的范数

矩阵APnm是线性映射（当nm时为线性变换）：PmPn的一种表现形式。因此，除了可以把矩阵看做向量而定义其范数外，更为基本、更为重要的是表征其线性映射的算子范数（operator norm），以APnn的情况为例：

Asupx0Axx （１．３）

其中（１．３）右端的范数·是赋范空间Pn，·中向量的范数，由矩阵算子范数的定义（１．３）容易证明（对映像大小的估计）不等式：

AxAx， xPn， （１．４）

称满足不等式（１．４）的矩阵范数是与对应的向量范数相容的。 例１．５ 常用的矩阵范数： １、１—范数（列范数）：

na，j1，2，，n A1max； ij1jni1２、２—范数（谱范数）： A2maxATA； ３、—范数（行范数）： Anmaxaij，i1，2，，n； 1ini1上述三种范数是如下定义的矩阵p—范数的特例：

４、由向量的p—范数：xpxk，1p，定义：

k1np1pApsupAxp （１．５）

x0５、Ｆ—范数（Frobenius）: AF2aij； i1j1nn12例 设A12，求范数A1,A2,A,AF. 34解 A1max{13,24}6.

Amax{12,34}7. AF14916305.4772.

由于ATAAT1014，所以1420Amax15221,1522115221，因此A2152215.4650.

向量和矩阵范数的MATLAB实现

在MATLAB中，norm( )函数求向量与矩阵的范数，其命令格式为norm(X,p).

当X为向量或矩阵时，norm(X,p)表示向量或矩阵X的p范数，例如，norm(X,1)表示X的1范数，norm(X,2)表示X的2范数，

norm(X,inf)表示X的范数.对于矩阵X，norm(X,'fro')表示X的F范数.

矩阵的条件数与直接法的误差分析

解线性方程组的直接法产生误差的主要原因：1）不同的算法及舍入误差的影响；2）方程组本身固有的问题（病态或良态）本节我们将分析方程组的状态并估计算法的误差，即原始数据扰动对解的影响.

考虑n阶线性方程组Axb，其中A为非奇异矩阵.

由于A（或b）的数值是测量得到的，或者是计算的结果，在第一种情况下A（或b）常带有某些观测误差，在后一种情况A（或

b）又包含有舍入误差.因此我们处理的实际矩阵是AA（或bb），下面我们来研究数据A（或b）的微小误差对解的影响.首

先考虑一个例子.

6x182例 设方程组Axb，即，它的精确解为(1,1)T. 26.00001x28.00001现在考虑系数矩阵和右端项的微小变化对方程组解的影响，即考察方程组

6x18225.99999x8.00002，

2其解变为(10,2)T.这种现象的出现是完全由方程组的性态决定的. 定义：如果矩阵A或常数项b的微小变化，引起方程组Axb解的巨大变化，则称此方程组为病态方程组，矩阵A称为病态矩阵（相对

于方程组而言），否则称方程组为良态方程组，矩阵A称为良态矩阵.需要一种能刻划矩阵和方程组“病态”程度的标准. 线性方程组的误差分析 设线性方程组为

Axb，

其中ARnn，x,bRn，且A非奇异.x*为精确解，x为解的误差，记

xxx.

设A为A的误差，b为b的误差.下面分别讨论x与A，b的关系.

b有误差而A无误差情形

将带有误差的右端项和带误差的解向量代入方程组，则

A(xx)bb.

由于Axb，而得到xA1b，从而xA1b. 另一方面，由(3.24)式取范数，有bAx，即

A1(b0).可得 xb定理 设A是非奇异矩阵，Axb0，且A(xx)bb，则有误差估计式

xxcond(A)bb,

其中cond(A)A1A称为方阵A的条件数.

【注】 （1） 解的相对误差是右端项b的相对误差的cond(A)倍；

（2） 如果条件数越大，则解的相对误差就可能越大；

（3） 条件数成了刻划矩阵的病态程度和方程组解对A或b扰动的敏感程度.

【定义7】 称条件数很大的矩阵为“病态”矩阵；称病态矩阵对应的方程组为病态方程组.反之，则称矩阵为良态矩阵，对应的方程组为良态方程组. A及b都有误差的情形

定理7 设方程组Axb0，A为非奇异矩阵，A及b都有误差，若A的误差A非常小，使A1A1，则有误差估计式

xx*bA. (3.26) AAb1cond(A)Acond(A)例 已知方程组Axb中

121A314131415113124471,b102， 560137660若b(0.001,0.001,0.001)T时，估计解的相对误差.

解 由于b估计

xx13102 108.3333,b120.001,由式（3.25）有误差

20150.0010.0186，

1310212xx比右端项的相对误差

bb9.2308106扩大了2015倍.

例 下列希尔伯特（Hilbert）矩阵是一族著名的病态矩阵.

1/21/n11/21/31/(n1). Hn1/n1/(n1)1/(2n1)用MATLAB函数计算条件数Hn.输入 for n=3:8 cond(hilb(n)) end

得到3至8阶希尔伯特矩阵的条件数分别为：

524.0568 1.5514e+004 4.7661e+005 1.4951e+007 4.7537e+008 1.5258e+010

由此可见，随着n的增加，Hn的病态可能越严重.Hn常出现在数据拟合和函数逼近中.

对于病态方程组，通常的方法无法得到它的准确解，需要采用一些特殊的处理方法. 病态方程组的处理

对于病态方程组，可采用高精度的算术运算，如双精度或扩充精度，以改善或减轻方程组的病态程度.如果用无限精度运算即不存在舍入误差的话，即使条件数很大，也没有病态可言.我们也可对病态方程组作预处理，使改善方程组系数矩阵的条件数.

例 设方程组Axb，即

1300064/10220003/106191000x2000， 3/1062/106试验证其为病态方程组，且对其作预处理PAxPb，使

cond(PA)cond(A).

解 （1）用MATLAB函数计算系数矩阵的条件数，得

cond(A)1.92011010，显然方程组为病态方程组.

（2） 令Pdiag(1,103,106)，使PAxPb，其中

1291,Pb2， PA3214323则有cond(PA)87.354cond(A).显然，经过预处理后的方程组

PAxPb是良态的.

奇异值分解法解病态方程组.

奇异值分解（Singular-Value Decomposition）简称SVD，对于n阶矩阵A，必存在正交矩阵U，V和对角阵S，使得A有奇异值分解

AUSVT. （3.27）

奇异值分解非常有用，对于矩阵n阶矩阵A，存在n阶正交矩阵

U,V，n阶对角阵S，满足AUSVT。U和V中分别是A的奇异向

量，而S是A的奇异值。AA的正交单位特征向量组成U，特征值组成STS，ATA的正交单位特征向量组成V，特征值（与AAT相同）组成STS。因此，奇异值分解和特征值问题紧密联系（注：奇异值分解对ARmn也适用）。

在MATLAB中，函数svd( )表示矩阵的奇异值分解，其命令格式为

T[U,S,V]=svd(A)

其中，U，V为正交矩阵，S为对角阵.例如，求4阶Hilbert矩阵的奇异值分解，输入

[U,S,V]=svd(hilb(4)) 得到

U = -0.7926 0.5821 -0.1792 -0.0292 -0.4519 -0.3705 0.7419 0.3287 -0.3224 -0.5096 -0.1002 -0.7914 -0.2522 -0.5140 -0.6383 0.5146

S = 1.5002 0 0 0 0 0.1691 0 0 0 0 0.0067 0 0 0 0 0.0001

V = -0.7926 0.5821 -0.1792 -0.0292 -0.4519 -0.3705 0.7419 0.3287 -0.3224 -0.5096 -0.1002 -0.7914 -0.2522 -0.5140 -0.6383 0.5146 矩阵S对角线上的元素称为奇异值，由大到小排列.