解线性方程组的直接法

主讲教师：高小辉E-mail:fzlcstar@126.com第四章线性方程组AX=B的数值解法(TheSolutionofLinearSystemsAX=B)1求解4.1引言许多实际问题可归结为线性（代数）方程组机械设备、土建结构的受力分析；经济计划输电网络、管道系统的参数计算；企业管理大型的方程组需要有效的数值解法。数值解法的稳定性和收敛性问题需要注意。2小行星轨道计算问题3天文学家要确定一小行星的轨道,在轨道平面建立以太阳为原点的直角坐标系.在坐标轴上取天文单位(地球到太阳的平均距离),对小行星作5次观察,测得坐标数据

(x1,y1),(x2,y2),(x3,y3),(x4,y4),(x5,y5)将数据代入椭圆的一般方程:a1x2+a2xy+a3y2+a4x+a5y+1=0得a1xk2+a2xkyk+a3yk2+a4xk+a5yk+1=0(k=1,2,3,4,5)4即求解方程组可确定椭圆方程(小行星轨道方程)5对一般线性方程组:AX=b,其中当系数矩阵A的行列式|A|≠0时,则方程组有唯一解.6求解线性方程组:AX=b的一般过程:输入：A，b解方程组算法输出：

X直接法：经过有限步算术运算求得精确解迭代法：从初始解出发，逐步求出近似解来逼近在求解小型(未知数较少)方程组时,直接法很有效.在求解大型方程组时,迭代法是最有效的方法.74.2高斯消去法(GaussianElimination

)1.上三角线性方程组(Upper-triangularLinearSystem)

定义4.1NN矩阵A=[aij]中的元素满足对所有i>j,有aij=0,则称NN矩阵A=[aij]为上三角矩阵。如果A中的元素满足对所有i<j,有aij=0,则称NN矩阵A=[aij]为下三角矩阵。4.2.1顺序高斯消去法89AX=B上三角线性方程组表示为：a11x1+a12x2+a13x3+…+a1n-1xn-1+a1nxn=b1a22x2+a23x3+…+a2n-1xn-1+a2nxn=b2a33x3+…+a3n-1xn-1+a3nxn=b3………….an-1n-1xn-1+an-1nxn=bn-1

annxn

=bn2.回代(BackSubstitution)设AX=B是上三角线性方程组，如果：akk0,k=1,2..n,则方程组存在唯一解。10证明：a11x1+a12x2+a13x3+…+a1n-1xn-1+a1nxn=b1a22x2+a23x3+…+a2n-1xn-1+a2nxn=b2a33x3+…+a3n-1xn-1+a3nxn=b3………….an-1n-1xn-1+an-1nxn=bn-1

annxn

=bn唯一用归纳法可证明xn-1,xn-2….x1是唯一的11例4.2：证明下列线性方程组无解

4x1–x2+2x3+3x4=207x3-4x4=-76x3+5x4=43x4=6例4.1：利用回代法求解线性方程组

4x1–x2+2x3+3x4=20–2x2+7x3-4x4=-76x3+5x4=43x4=6

x4=6/3=2x3=(4-5x4)/6=-1x2=(-7-7x3+4x4)/-2=-4x1=(20+x2-2x3-3x4)/4=3x4=2x3=-1x3=1/7a22=00x2

+12例4.3：证明下列线性方程组有无穷解

4x1–x2+2x3+3x4=20

0x2+7x3-0x4

=-76x3+5x4=43x4=6x4=2x3=-1x3=-1x2=4x1-16无穷解a22=013复习如A是NN矩阵a11a12a13…a1na21a22a23…a2n

……an1an2an3…annA=则A的行列式为：划掉A的第i行和第j列后的行列式第i行展开第j列展开1438-45-17-69A=j=2,det(A)=-3+5+6=77-4-179287928-4-1i=1,det(A)=2-3+8

=2(45-6)-3(-36+7)+8(24-35)=775-1-69-4-179-457-6153.

如果NN矩阵A=[aij]是上三角矩阵或下三角矩阵，则：

定理：A是NN方阵，线性方程组AX=B有唯一解当且仅当det(A)016高斯消元法：思路首先将A化为上三角阵，再回代求解。=174初等变换(ElementaryTransformation)下列三种变换可使一个线性方程组变换成另一个等价的线性方程组交换变换：对调方程组的两行比例变换：用非零常数乘方程组的某一行替换变换：将方程组的某一行乘一个常数再加到另一行

18例4.4：求抛物线y=A+Bx+Cx2,它经过三点(1,1),(2,-1),(3,1)。

(2)-(1)A+B+C=1(1)B+3C=-2(2)2B+8C=0(3)(3)-2(2)A+B+C=1(1)B+3C=-2(2)2C=4(3)C=2,B=-8,A=7抛物线方程y=7-8x+2x2A+B+C=1(1)A+2B+4C=-1(2)A+3B+9C=1(3)(3)-(1)(3-1)(3-2)19上述求解方程组的方法就是高斯(Gauss)消元法。从式(4-1)到(4-2)的过程称为消元过程而由(4-2)求出C、B、A的过程称为回代过程。因此用高斯消去法求解性方程组要经过消元和回代两个过程。20可将线性方程组AX=B的系数存在N(N+1)增广矩阵中a11a12a13…a1nb0a21a22a23…a2nb1

……an1an2an3…ann

bn[A|B]=215初等行变换(ElementaryRowOperation)对增广矩阵进行如下变换可得到一个等价的线性方程组交换行变换：对调矩阵两行比例行变换：用非零常数乘矩阵某一行所有的元素替换行变换：将矩阵的某一行的所有元素乘一个常数再加到另一行对应的元素上去。

224.2.2主元素（Pivot）高斯消去法

系数矩阵A中的元素arr用来消去akr(k=r+1,r+2,…,N),arr为第r个主元，第r行为主元行。例4.6

用增广矩阵表示下列线性方程组，并求等价上三角线性方程组和方程组的解。x1+2x2+x3+4x4=132x1+0x2+4x3+3x4=284x1+2x2+2x3+x4=20-3x1+x2+3x3+2x4=6

23x1+2x2+x3+4x4=132x1+0x2+4x3+3x4=284x1+2x2+2x3+x4=20-3x1+x2+3x3+2x4=6

1214

132043

28

422120-3132

60761445-主元行＊2-主元行＊4-主元行＊-31214

130–42-5

2

0–6–2–15-3224-主元行＊-1.75-主元行＊1.507614451214

130–42-5

2

0–6–2–15-32009.55.2548.51214

130–42-5

2

00–5–7.5-350000–9-181214

130–42-5

2

00–5–7.5-35-主元行＊1.9回代：x4=2,x3=4,x2=-1,x1=3消元过程25有回代的高斯消去法(GaussianEliminationwithBackSubstitution)如果A是NN非奇异矩阵（存在A-1），则存在线性方程组UX=Y与线性方程组AX=B等价，这里U是上三角矩阵，并且akk0。当构造出U和Y后，可用回代法求解UX=Y，并得到方程组的解X。26消元记27Step1：设，计算因子将增广矩阵第i行mi1

第1行，得到其中28一般地，假定已完成了(k-1)步消元，即已将[A(1)|B(1)]转化为以下形式：29Stepk：设，计算因子共进行？步n

1且计算30回代31消元法步骤(以4阶为例)：增广矩阵

32计算3个消元因子(乘子向量)(1)设a11≠0,[m21

m31

m41]T=[a21

a31

a41]T/a11

用-m21乘矩阵第一行后加到矩阵第二行;用-m31乘矩阵第一行后加到矩阵第三行;用-m41乘矩阵第一行后加到矩阵第四行;消元33用-m31乘矩阵第一行后加到矩阵第三行;消元34用-m41乘矩阵第一行后加到矩阵第四行;消元35第二轮消元后(2)设a22(1)≠0,计算第二个消元因子[m32

m42]T=[a32(1)

a42(1)]T

/a22(1)

用-m32乘矩阵第二行后加到矩阵第三行;用-m42乘矩阵第二行后加到矩阵第四行;36第三轮消元后用-m43乘矩阵第三行后加到矩阵第四行;(3)设a33(2)≠0,计算第三个消元因子:m43=a43(2)/a33(2)对应三角形方程组37消元过程中,消元因子可排列为一矩阵用回代算法即可得出方程组的解.38由初等变换的知识可得.L×U=A在消元过程中,我们得到两个矩阵此即为矩阵A的LU分解.39高斯消元法的主要缺陷1.计算消元因子mik=aik/akk当分母akk为零时,计算无法进行;2.分母akk的绝对值非常小时,也将会对计算结果误差产生很大影响解决办法40选主元以避免app(p)=0

如果app(p)=0,寻找第p行下满足akp(p)

0的第一行，设为第k行，然后交换行k和行p，使新app(p)

0。如果app(p)≠0则不交换。此方法称为平凡选主元方法。41将列主元所在行与第k行交换后,再按上面的高斯消元法进行下去,此法即称为列主元消元法。选|aik(k)|(i=k,…,n)最大的一个(列主元)选主元以减少误差在每一步消元过程中取系数子矩阵的第一列元素中绝对值最大者作主元。42例4.7:

值x1=x2=1.000是如下方程组的解：

1.133x1+5.281x2=6.41424.14x1-1.210x2=22.93

用4位有效数字及选主元的高斯法求方程的解。行2-行1＊24.14/1.133:1.133x1+5.281x2=6.414-113.7x2=-113.8

x2=1.001,x1=0.995643例4.8：24.14x1-1.210x2=22.931.133x1+5.281x2=6.414

行2-行1＊1.133/24.14:24.14x1-1.210x2=22.935.338x2=5.338x2=1.000,x1=1.00044例用列主元法解第一列中绝对值最大是10，取10为主元45第二列的后两个数中选出主元2.546X3=6.2/6.2=1X2=(2.5-5X3)/2.5=-1X1=(7+7x2-0x3)/10=0X1=0X2=-1X3=1N阶方程组第k轮消元时，选第k列的后(n-k+1)个元素中绝对值最大者作主元。474.3矩阵三角分解法（TrianguarFactorization）定义

如果非奇异矩阵A可表示为下三角矩阵L和上三角矩阵U的乘积A=LU，则A存在一个三角分解。4.3.1高斯消去法和矩阵三角分解法a11x1+a12x2+a13x3+…+a1n-1xn-1+a1nxn=b1a21x1+a22x2+a23x3+…+a2n-1xn-1+a2nxn=b2a31x1+a32x2+a33x3+…+a3n-1xn-1+a3nxn=b3……an1x1+an2x2+an3x3+…+ann-1xn-1+annxn

=bn48a11a12a13…a1na21a22a23…a2na31a32a33…a3n……an1an2an3…ann=100…0

m2110…0

m31m321…0

……mn1mn2mn3…1u11u12u13…u1n0u22u23…u2n00

u33…u3n……00

0

…unn49a11a12a13…a1na21a22a23…a2na31a32a33…a3n……an1an2an3…annx1x2x3xn=b1b2b3bn100…0

m2110…0

m31m321…0

……mn1mn2mn3…1u11u12u13…u1n0u22u23…u2n00

u33…u3n……00

0

…unnx1x2x3xn=b1b2b3bny1y2y3yn50100…0

m2110…0

m31m321…0

……mn1mn2mn3…1y1y2y3yn=b1b2b3bnu11u12u13…u1n0u22u23…u2n00

u33…u3n……00

0

…unnx1x2x3xn=y1y2y3yn

UX=YLY=bYX51由初等变换的知识可得.L×U=A在高斯消元过程中,我们得到两个矩阵此即为矩阵A的LU分解.52例4.3.2

求解：x1+2x2+4x3+x4=212x1+8x2+6x3+4x4=523x1+10x2+8x3+8x4=794x1+12x2+10x3+6x4=821241

2864

310884121061000

2100

311041211241

04-22

00-23000-6A===LU531000

2100

311041211241

04-22

00-23000-6x1x2x3x4=215279821000

2100

31104121y1y2y3y4=21527982y1=21y2=10y3=6y4=-2454x4=4x3=3x2=2x1=11241

04-22

00-23000-6x1x2x3x4=21106-24554.3.2直接三角分解法例4..3.1使用高斯消去构造下列矩阵的三角分解：43-1

-2-45126A=

00-0.5100.250143-1

0-2.54.501.256.2500010001=43-1

-2-45126行2-行1*(-2/4)行3-行1*(1/4)

00-0.5100.25–0.5143-1

0-2.54.5008.5行3-行2*(1.25/-2.5)43-1-0.5-5–0.58.556定理设无行交换变换的高斯消去法可求解一般线性方程组AX=B，则矩阵A可分解为一个下三角矩阵L和一个上三角矩阵U的乘积：A=LU，而且L对角线元素为1，U的对角线元素非零。得到L和U后，可通过如下步骤得到X：

1、利用前向替换法对方程组LY=B求解Y2、利用回代法对方程组UX=Y求解X574.4置换矩阵

可能一个非奇异矩阵不能直接分解为A=LU例4.4.1证明下列矩阵不能直接分解为A=LU：126

48-1-235A=

58126

48-1-235A=

=100

m2110m31m321u11

u12

u13

0u22

u2300u33

1=u114=m21u11=m21-2=m31u11=m31

2=u128=m21u12+u22=4*2+03=m31u12+

m32u22

=(-2)*2+m32*0

=-459定义一个NN置换矩阵P是一个在每一行和每一列只有一个元素为1，而其它元素为0的矩阵。P的列是单位矩阵行的置换，可表示为：P=[E’k1E’k2……

E’kn

]’Pij=j=ki0其它情况如：0100

100000010010

P==[E’2E’1E’4

E’3]’60定理设P=[E’k1E’k2……

E’kn

]’是一个置换矩阵。PA是一个新的矩阵，它的行是将A中的行按行k1A,行k2A，…行knA调整顺序后形成的。例4.4.2设A为44矩阵，P为上式中的置换矩阵，则PA矩阵中的行是将A中的行调整顺序后形成的，顺序为第1，2，3，4行对应于A中的第2，1，4，3行。

610100

100000010010

a11a12a13a14a21a22a23a24a31a32a33a34a41a42a43a44=a21a22a23a24a11a12a13a14a41a42a43a44a31a32a33a3462定理如果A是非奇异矩阵，则存在一个置换矩阵P，使得PA存在三角分解：PA=LU例4.4.3设A为非奇异矩阵126

48-1-235A=126

48-1-235PA=100

001010=126

-23548-1126

071700-25可分解U634.5扩展高斯消去过程定理（非直接分解：PA=LU）设A是一NN矩阵。假设高斯消去法可求解经过行变换的一般线性方程组AX=B。则存在一个置换矩阵P，使得PA可分解为一个下三角矩阵L和一个上三角矩阵U：PA=LU。而且L的主对角线元素为1，U的主对角线元素非零。可用如下4步求出X：1、构造矩阵L，U，P2、计算列向量PB3、用前向替换法对方程组LY=PB求解Y4、用回代法对方程组UX=Y求解X如何求P64a11a12a13a14a21a22a23a24a31a32a33a34a41a42a43a440010

010010000001P=a31a32a33a34a21a22a23a24a11a12a13a14a41a42a43a44通过选列主元0010

000110000100P=PAX=PBAX=B65x1+2x2+4x3+x4=212x1+8x2+6x3+4x4=523x1+10x2+8x3+8x4

=794x1+12x2+10x3+6x4=821241

2864

31088412106A=1、构造矩阵L，U，P4121062864

3108812410001

010000101000P=4121061/2211

3/410.53.51/4–11.5–0.5664121061/2211

3/410.