计算方法与实习-线性方程组数值解法

第三章线性方程组数值解法第一节引言工程中几乎有一半的问题涉及到线性方程组的求解。设线性方程组a11x1+a12x2+……+a1nxn=b1

a21x1+a22x2+……+a2nxn=b2

AX=B……an1x1+an2x2+……+annxn=bn

a11a12……a1n

X1

b1A=a21a22……a2nX=X2B=b2……...…an1an2……annX3bn

A称为方程组的系数矩阵，当是n阶非奇异矩阵时，既A0，此时方程组有唯一解。

X阵是解向量，B阵是常向量。在线性代数中学过用克莱姆法则求解，它是一种直接法（属于解析法），但随着n计算工作量

直接法：经有限次计算得准确解（在无舍入误差下），实际上舍入误差客观存在，得到的依然还是近似解。方法有主元|消去法、三角分解法。计算量小、精度较|高，但程序设计复杂。解线性方程组的数值解法|迭代法：将问题构成一个无穷序列，逼近|准确解。方法有雅可比法、高斯-赛得尔法|算法简单、易于编程，但存在收敛性及收敛|速度的|问题，只对具有某些性质的系数矩阵的方程组才适用。第二节消去法一、高斯消去法1.基本思想属于直接法，一般由“消元”和“回代”两过程组成，基本思想是将系数矩阵转化成“三角系数矩阵”。以三阶方程组说明如下2x1+4x2-2x3=6x1-x2+5x3=04x1+x2-2x3=224-2设A（1）=1-15=（aij(1)）=(aij)=A(i=1,2,3,j=1,2,3)41-2

6B(1)=0=(bi)(i=1,2,3)2令li1=ai1(1)/ai1(2)(i=1,2,3)，称li1为消元因子，用li1消去除第一个方程外其余各方程的x1，即将第一个方程乘以-li1后加到第i个方程上去（i=2，3），则有24-2x160-36x2=-30-72x3-10此时记为A（2）X=B(2)

其中aij(2)=aij(1)-li1a1j(1)(i=j=2,3)

bi(2)=bi(1)-li1b1(1)令li2=ai2(2)/a22(2)(i=3)，将第二个方程乘以-li2后加到第i个方程上去（i=3）则有24-2x160-36x2=-300-12x33记为A（3）X=B(3)

其中aij(3)=aij(2)-li2a2j(2)(i=j=3)

bi(3)=bi(2)-li2b2(2)(i=3)现在A（3）已成上三角阵，上述过程称之为消元过程，下面进行回代过程，所以从第三个方程可解出x3后，回代入第二个方程得x2，再将x2、x3回代到第一个方程解出x1。在上述过程中要求aii0，若aii=0，则可通过行、列变换A来满足aii0，aii称为主元素。aii越小，作除数引起的误差就越大，所以尽可能选用maxaii，这样就有了全主元、列主元消去法。现讨论一般n阶方程组AX=B记A（1）X=B(1)，其中A(1)=(aij(1))=（aij）=AB(1)=(bi(1))=（bi）T=B消元过程：设a110，将第一个方程乘以-li1li1=ai1(1)/ai1(1)(i=2,3,……,n)

加到第i个方程上去，就消去了第2~第n个方程中的x1，得到等价方程组。a11（1）a12（1）……a1n（1）

x1b1(1)0a22（2）……a2n（2）

x2b2(2)……...=…0an2（n）……ann（n）

xnbn

(n)

A（2）X=B（2）aij(2)=aij(1)-li1a1j(1)(i=j=2，3，4……n)bi(2)=bi(1)-li1b1(1)(i=2，3，……，n)实际上，令10……0

-l211……0L1=…01…-ln10…1

则10……0

a11（1）a12（1）……a1n（1）

-l211……0a21（1）a22（1）……a2n（1）

L1A(1)=…01………-ln10…1

an1（1）an2（1）……ann（1）a11（1）a12（1）……a1n（1）

-l21a11（1）+a21（1）-l21a12（1）+

a22（1）……-l21a1n（1）

+a2n（1）=……-ln1

a11（1）+an1（1）-ln1a12（1）+an2（1）……-ln1a1n（1）+

ann（1）

a11（1）a12（1）……a1n（1）

0a22（2）……a2n（2）

=……=A（2）0an2（2）……ann（2）

1

b1（1）

b1（1）

-l211b2（1）

-l21b1（1）+b2（1）

L1B(1)=…………=……=B(2)-ln1

1bn（1）

-ln1

bn（1）+bn（12）

即用初等矩阵L1左乘A（1）和B（1）后得到A（2）和B（2）。接着要消去第3~第n个方程的x2，即第2个方程乘以-li2li2=ai2(2)/a22(2)(i=3,4,……,n)，且a22(2)

0加到第i个方程上去（i=3,4,……,n），得a11（1）a12（1）……a1n（1）

x1b1(1)0a22（2）……a2n（2）

x2b2(2)|a33（3）…a3n（3）x3=b2(3)|…….…..….…..||…||…|00……ann（3）

xnbn

(3)

A（3）X（3）其中aij(3)=aij(2)-li2a2j(2)(i,j=3,4,…,n)

bi(3)=bi(2)-li2b2(2)(i=3,4,…,n)实际上

1

01令L2=0-l321……0-ln2

……

1则L2A（2）=A（3），L2B（2）=

B（3）设第k-1步已作过k-1次消元，则等价方程组

A（k）X=B（k）

a11（1）a12（1）……a1n（1）

b1(1)0a22（2）……a2n（2）

b2(1)A（k）=……b（3）=…0akk（k）…akn（k）bk(k)

………00ank（k）…ann（k）bn

(k)作第k步消元，若akk(k)

0，则将第k个方程乘以-liklik=aik(k)/akk(k)(i=k+1,k+2,……n)加到第i个方程上(i=k+1,k+2,……n)，得A（k+1）x=B（k+1）其中aij(k+1)=aij(k)-likakj(k)(i，j=k+1,k+2,…,n)

bi(k+1)=bi(k)-likbk(k)(i=k+1,k+2,…,n)令1

01Lk=0-lk+1k1……0-lnk

……

1则Lk

A（k）=A（k+1），Lk

B（k）=

B（k+1）按上述步骤直到第n-1步做完，到得

A（n）X=B（n）即a11（1）a12（1）……a1n（1）

x1b1(1)0a22（2）……a2n（2）

x2b2(1)……...=…00……ann（n）xnbn

(n)现在的A（n）是一个上三角系数矩阵，消元过程结束。逆序回代过程xn=bn(n)/ann（n）

xn-1=(bn-1(n-1)-an-1n（n-1）xn)/an-1n-1（n-1）……x1=(b1(1)-a1n（1）xn-a1n-1（1）xn-1-…-a12（1）x2)/a11（1）

n∴xk=(bk(k)-∑akj（k）xj)/akk（k）j=k+1由于A（n）=Ln-1Ln-2……L2L1A（1）

B（n）=Ln-1Ln-2……L2L1B（1）

令A（n）=U，B（n）=Y，则u11u12……u1ny10u22……u2ny2

U=……Y=…0……unnyn所以Ln-1Ln-2……L2L1A=U。因为det(

Lk)=1，故Lk-1存在1

01Lk-1

=0

lk+1k1……0

lnk

……

1有A=L1-1L2-1……Ln-2-1Ln-1-1U令L=L1-1L2-1……Ln-2-1Ln-1-11

l211=

l31l321……ln1

ln2……

lnn-11所以A=LU同理L1

L2……Ln-2

Ln-1B=YB=L1-1L2-1……Ln-2-1Ln-1-1Y=LY结论：经消元过程后将A分解成单位下三角阵L和上三角阵U的乘积，这种分解成简单矩阵之积的方法称为矩阵分解法。因为

AX=B

LY=B所以LUX=

B

UX=YLY=B求解AX=B求解

UX=Y这两个三角方程组问题，难度降低。先从第一个方程解出Y，再在第二个方程求出X即从

b1

1y1

b2

l211y2

…

=……...

bnln1

ln2……

lnn-11yn

k-1解得y1=b1,，

yk=bk-lkjyj（k=2,…,n）j=1这实际上就是消元过程。再求UX=Yu11u12……u1nx1y10u22……u2nx2y2………=…0……unnxnyn得xn=yn/

unnn

xk=（

yk-

ukj

xj）/ukk（k=n-1,n-2,…,2,1）j=k+1这也是回代过程。结论：（1）高斯消去法分消元、回代两过程。（2）从矩阵分解角度看，消去是解一个下三角方程组，回代是

解一个上三角方程组。（3）消去法顺利进行必须满足akk

(k)

0,（k=1,2,…,n），若出现akk

(k)

=0，则可交换行列后再进行消元。2.运算量（运算次数）在第k步消元中，计算lik=aik(k)/akk(k)(i=k+1,k+2,……n)，需作n-(k+1)+1=n-k次除法计算aij(k+1)=aij(k)-likakj(k)(i,j=k+1,…,n)需作（n-k）2次乘法和（n-k）2次减法计算bi(k+1)=bi(k)-likbk(k)(i=k+1,…,n)需作（n-k）次乘法和（n-k）次减法。所以消元过程总运算量为：

n-1乘除法次数：M1=[（n-k）+(n-k)2+(n-k)]k=1

n-1=[2（n-k）+(n-k)2]=(n-1)n(2n-1)/6+2(n-1)n/2k=1=n3/3+n2/2-5/6*n=O(n3)（以n3为数量级的欧函数O(x)）n-1加减法次数S1=[（n-k)2+(n-k)]=(n-1)n(2n-1)/6+k=1

(n-1)n/2=n3/3-n/3=O(n3)而回代过程的运算量

nn-1M2=i=n(n+1)/2=O(n2),

S2=i=n(n-1)/2=O(n2)i=1i=1总量M=M1+

M2=n3/3+n2-n/3=O(n3)

S

=S1+

S2=n3/3+n2/2

-5n/6

=O(n3)乘除法耗时大大多于加减法耗时，故高斯消元法的计算量为O(n3)。二、列主元高斯消去法及程序设计（1）方法以akk

(k)0，（k=1,2,…,n）为条件，若akk

(k)

0会使误差增大，结果失真。可以先选列主元，然后消元，现用增广矩阵a11（1）a12（1）……a1n（1）

b1(1)0a22（2）……a2n（2）

b2(2)………00……akn（k）

bk(k)00……ak+1n（k）

bk+1(k)………00……ann（k）

bn

(k)选出第k列的主元素，即maxaik

(k)=amk

(k)kin若amk

(k)<则停止，否则继续a11（1）a12（1）…… a1n（1）

b1(1)0a22（2）…… a2n（2）

b2(2)……… 00…amk（k）…akn（k）

bm(k)00…ak+1k（k）…ak+1n（k）bk+1(k)………00……ank（k）…ann（k）

bn

(k)(k=1,2,3,…,n)从而使消元因子lik

=

aik(k)/akk(k)

1，稳定性提高。例0.00287.1387.15

4.453-7.2637.27因为4.453>0.002

故选a21为主元素，r1r24.453-7.2637.274.453-7.2637.27| || |0.00287.1387.15087.1387.13解得x1=1.000，

x2=10有不必选主元情况，即若系数矩阵是对称严格的对角占优阵，从而能够保证对角线的元素满足下列条件

n aii

>

aij

j=i j≠i高斯列主元消去法N-S图消元次数K=1k<=n?选列主元amk|amk|

<epsNrmrk进行第k步消元输出方程组的解X建立增广矩阵A(n,n+1)m=kYNYk=k+1回代求解停止三、高斯-约当消去法（Gauss-Jordan）这是一种无回代的方法，使计算工作量减小，在第k次消元时不仅把akk(k)

化成1，不仅将aik（i=k+1，……，n）化成0，而且也将aik（i=1，2，……，k-1）也化为0，最后消元过程结束时的系数矩阵为11

……1得消元公式akj

(k)=akj(k-1)/akk(k-1)（j=1,2,…,n+1）aij

(k)=aij

(k)-aik(k-1)akj(k)（ik,i=1,2,…,n，j=1,2,…,n+1）第三节矩阵的直接分解及其在解方程组中的应用一、矩阵分解的紧凑格式当A0且aij

(k)0时，作LU分解（经消元后分解），另可由aij及递推关系直接得L及U阵的各元素。

A=LUa11a12……a1n1u11u12…u1j…u1na21a22……a2nl211u22..u2j..u2n………………ai1ai2……ain=li1li2…1ii+1………………………an1

an2……annln1

ln2……1……unn比较两边第一行元素，有

a11=

u11，a12=u12……

a1n=u1n 即

u1j=

a1j(j=1,2,….,.n)比较两边第一列元素，有

a21=

l21u11，a31=l31u11，…...an1=

ln1u11 即ai1=li1u11(i=2,3,…,n)所以

li1=

ai1/u11(i=2,3,…,n)由此得

a1j（j=1,2,…n)-->

u1j（j=1,2,…n)li1=

ai1/u11

ai1(i=2,3,…,n)再比较第2行剩下的元素(j=2,….,.n)

a2j=

l21u1j+

u2j(j=2,….,.n) 即u2j=

a2j-

l21u1j(j=2,….,.n)而第二列剩下的元素(i=3,4,…,n) ai2=li1u12+

li2u22(i=3,4,…,n) 即li2=（ai2-li1u12)/

u22(i=3,4,…,n)所以通过 a2j(j=2,3,….,.n) l21--->u2j

u1j(j=2,3,….,.n)

ai2--->li2

li1

u12设求出U的第一行到第K-1行元素u11u12…u1k…u1j…u1n与L1

00u22...

u2k...u2j…

u2nl2110……0……l32……0...uk-1k…uk-1j…uk-1n……100……lk1

lk2…lkk-11…0……00…………0……ln1

ln2.....

lnk-1现求

ukj(j=k,k+1,…,n)及lik

(i=k+1,…,n)nk-1因为akj

=

lkq

.uqj=

lkq.uqj

+ukj

(j=k,k+1,…,n)q=1q=1

k-1所以ukj=

akj-

lkq

.uqj

(j=k,k+1,…,n)q=1

nk-1因为aik=

liq

.uqk=

liq

.uqk

+

lik

ukk(i=k+1,…,n)

q=1q=1

k-1所以lik=(aik-

liq

.uqk

)/ukk(i=k+1,…,n)q=1这样就可以进行直接分解。又因为LY=B y1=b1 k-1

yk=bk

-

lkq

.yqk

q=1这与计算ukj具有相同的格式。实际上在具体计算时可以采用紧凑格式，见P55。二、追赶法1.方法

这种方法适用于A为三对角矩阵情况。即b1c1

x1d1

a2b2c2

x2

d2

a3

b3

c3……x3d3……

...=…an-1

bn-1

cn-1xn-1dn-1anbnxnd6

b1x1+

c1x2=d1

所以 x1=

d1/

b1-c1/

b1*x2=

p1-q1x2这里令 p1=

d1/

b1

q1=

c1/

b1

因为a2x1+

b2x2

+

c2x3=

d2，将x1代入得x2=

(d2-a2p1)

/

(

b2-a2q1)-c2x3/

(

b2-a2q1)

=

p2-q2x3其中

p2=(d2-a2p1)

/

(

b2-a2q1) q2=c2/

(

b2-a2q1)同理x3=

(d3-a3p2)

/

(

b3-a3q2)-c3x4/

(

b3-a3q2)

=

p3-q3x4其中

p3=(d3-a3p2)

/

(

b3-a3q2) q3=c3/

(

b3-a3q2)第k步时xk=

(dk

-akpk-1)

/

(

bk-akqk-1)-ckxk+1/

(

bk-akqk-1)

=pk-qkxk+1其中

pk=(dk-akpk-1)

/

(

bk-akqk-1) qk=ck/

(

bk-akqk-1)而第n步时xn=

pn=

(dn

-anpn-1)

/

(

bn-anqn-1)消元为i从小到大为，这是追，回代按xn--->

xn-1---……- -->

x1就是赶。2.程序设计三、消去法的应用1.解线性方程组系AX=BL(L=1,2,…,m)m个方程此时，增广矩阵为[A|B1B2…...

Bm

]消去过程相同，但回代时分别用BL（n)代，所得解向量Xi即为

常向量为BL时的解。2.求逆阵当A=(aij)n*n为满秩矩阵时，令X=A-1n*n，则因为

AA-1=AX=E，即相当于求解方程组系x111

x12

0

x1n0

Ax12=0,Ax22=0,…,Ax2n=0…...…………

x1n0

xn20

xnn

1

3.求行列式值将行列式转化成对角阵，则|A|为对角线元素之积。第四节迭代法及其收敛性利用逼近的基本思想，将

AX=B变形为

X=CX+F(C称为迭代矩阵)构造公式 X（k+1)=CX(k)+F当给定一组初值X（0）后代入公式得X（1），依次求得序列{X（k）}，当C不同时，方法有多种。收敛问题是limX(k)

=X*这是极限向量。

k即limxi(k)=xi*(i=1,2,…,n;k=1,2,…)k问题：迭代格式的建立、收敛条件、收敛速度和误差估计。P63-64

的定理3.2和定理3.3说明了收敛条件一、雅可比（Jacobi）迭代法线性方程组AX=B若|A|不等于0，且aii不等于

0（i=1,2,…,n)将方程组变形为x1=(b1-a12x2-a13x3…-amxn)/a11x2=(b2-a21x1-a23x3…-a2nxn)/a22…...xn=(bn-an1x1-an2x2…-ann-1xn-1)/ann

n即xi=-

aijxj/aii+bi/aii(i=1,2,…,n)

j=1ji给出X(0)=（x1（0）,x2（0）,x3（0）,…,

xi（0）,…,xn（0）)T,利用上式反复迭代，可得{X（k）}，若{X（k）}-->X*,则{xi（k）}-->xi*，xi*即为解，k为迭代次数。实际上若xi（k+1）-xi（k）<(i=1,2,……n),则停止迭代。xi（k+1）即为方程组的解。P67例7Jacobi法简单，只作矩阵的一次乘法，但需要保存X（k）来求X（k+1）。定理：

n若μ=max

aij/aii<1，则Jacobi法收敛。1in

j=1ji雅可比迭代法N-S图

k=1k<=kmax?e=0,i=1|sum-y[i]|>eNsum=sum-a[i][j]*x[j]读入n,kmax,eps，建立a(n,n),b(n),y[n]i≠jYNYi<=n?sum=b[i],j=1j<=n?j=j+1sum=sum/a[i][i]

e=|sum-y[i]|x[i]=y[i],输出k,x1,x2,…,xnN输出"无解"停止Yy[i]=sum,i=i+1e<eps?k=k+1结束二、高斯－赛德尔法及程序设计1.方法是改进的Jacobi法。当计算xi（k+1）时，x1（k+1）…...xi-1（k+1）已求出，所以可将x1（k+1）…...xi-1（k+1）及

xi+1（k）…...xn（k）代入迭代公式求出xi（k+1），从而加快迭代收敛速度，迭代次数也下降了。用一组存储单元存放X，即

i-1n

xi（k+1）=(bi－

aij

xi（k+1)－aij

xi

（k))/aiij=1j=1+1

（i=1，2，…，n）实际上若xi（k+1）-xi（k）<(i=1,2,……n),则停止迭代。xi（k+1）即为方程组的解。P69例8。高斯赛德尔法收敛条件有：

n aiiaij或A为对称正定阵。

j=1ji2.高斯－赛德尔迭代法N-S图

k=1k<=kmax?e=0,i=1|y-x[i]|>eNsum=sum-a[i][j]*x[j]读入n,kmax,eps，建立a(n,n),b(n)i≠jYNYi<=n?sum=b[i],j=1j<=n?j=j+1y=sum/a[i][i]

e=|y-x[i]|输出k,x1,x2,…,xnN输出"无解"停止Yx[i]=y,i=i+1e<eps?k=k+1结束3.化AX=B为易收敛形式 AX=B-->X=CX+D若