线性代数方程组的数值解法
![线性代数方程组的数值解法]()
n阶线性代数方程组的一般形式
⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧a11x1+a12x2+⋯+a1nxn=b1a21x1+a22x2+⋯+a2nxn=b2⋮an1x1+an2x2+⋯+annxn=bn
或写成矩阵——向量形式:Ax=b,其中A为系数矩阵,x为解向量,b为右端常向量,分别为
A=⎝⎜⎜⎛a11⋮an1⋯⋱⋯a1n⋮ann⎠⎟⎟⎞,x=⎝⎜⎜⎜⎜⎛x1x2⋮xn⎠⎟⎟⎟⎟⎞,b=⎝⎜⎜⎜⎜⎛b1b2⋮bn⎠⎟⎟⎟⎟⎞
直接消去法
顺序高斯消去法
设
akk(k)=0
对k=1,2,…,n-1计算
⎩⎪⎪⎪⎨⎪⎪⎪⎧mik=akk(k)aik(k)aij(k+1)=aij(k)−mikakj(k)bi(k+1)=bi(k)−mikbk(k)i,j=k+1,k+2,…,n
⎩⎪⎪⎨⎪⎪⎧xn=ann(n)bn(n)xi=(bi(i)−∑j=i+1naii(i)aij(i)xj)/aii(i)i=n−1,…,2,1
-
按消元规则运算后,对角线以下元素为0,故运算中对于对角线以下元素不进行计算,以减少计算量
-
对角线下元素对回代求解无影响,故将乘数放在该处,即aik/akk→aik,i=k+1,k+2,…,n以节省存储单元
-
对角线以上元素和常数项采用“原地”工作方式,即经变换后的元素仍放在原来的位置上
{aij−aikakj→aijbi−aikbk→bii,j=k+1,k+2,…,n
以节省存储单元
-
回代后的值仍被放在常数项存储单元
annbn→bn
(bi−j=i+1∑naijxj)/aii→bi,i=n−1,n−2,…,1
以节省存储单元,这时b1,b2,…,bn单元中存放的就是输出值x1,x2,…,xn
将第k步消xk时用的方程及其系数分别叫作第k步的主方程和主行,而xk的系数akk(k)叫作第k步的主元素,要求各步主元素不为0
![高斯消元法算法框图]()
定理
-
方程组系数矩阵的顺序主子式全不为零,则高斯消去法能实现方程组的求解
顺序主子式:设方程组系数矩阵A=(aij)n,其顺序主子式
Di=∣∣∣∣∣∣∣∣a11⋮ai1⋯⋱⋯a1i⋮aii∣∣∣∣∣∣∣∣,i=1,2,…,n
非奇异矩阵:|A|≠0
-
设方程组Ax=b,如果系数矩阵A为严格对角占优矩阵,则用高斯消去法求解时,主元素akk(k)全不为零
严格对角占优矩阵:矩阵A=(aij)n每一行对角元素的绝对值都大于同行其他元素绝对值之和
∣aii∣>j=1j=i∑n∣aij∣,i=1,2,…,n
-
高斯消元法求解n阶线性方程组共需乘除法次数近似为1/3n3
列主元高斯消去法
顺序高斯消元法有如下两弊端:
- 顺序消元过程中可能出现主元素akk(k)=0
- akk(k)≠0,但其绝对值很小,用它作为除数也会导致其他元素的数量级急剧增大和使舍入误差扩大,影响计算结果精度
为避免在消元过程确定乘数时所用除数是零或绝对值小的数,即零主元或小主元,在每一次消元之前,要增加一个选主元的过程,将绝对值大的元素交换到主对角线的位置上。根据交换的方法又分成全选主元和列选主元两种选主元方法
全选主元是当变换到第k步时,从右下角n-k+1阶子阵中选取绝对值最大的元素,然后通过行交换与列交换将其交换到akk的位置上,并且保留下交换的信息,以供后面调整解向量中分量的次序时使用
列选主元是当消元到第k步时,从k列的akk以下(包括akk)的各元素中选出绝对值最大的,然后通过行交换将其交换到akk的位置上。交换系数矩阵中的两行(包括常数项),只相当于两个方程的位置交换了,因此,列选主元不影响求解的结果
列选主元:设主元在第l(k≤l≤n)个方程,即
∣alk∣=k⩽i⩽nmax∣aik∣
若l≠k,将l和k方程互易位置,使新的akk成为主元,然后继续进行
列选主元比全选主元的运算量小,但一般可以满足精度要求,所以列选主元更常被采用
严格对角占优矩阵可保证akk(k)就是主元
![列选主元框图]()
定理
- 列主元消去法在解方程组时,还可求出系数行列式。设系数矩阵
A=⎝⎜⎜⎛a11⋮an1⋯⋱⋯a1n⋮ann⎠⎟⎟⎞
用列主元消去法将其转换为上三角形矩阵,对角线上元素为a11(1),a22(2),…,ann(n),于是行列式det(A)=(−1)ma11(1)a22(2)⋯ann(n)
其中m为所进行的行交换次数
高斯——若尔当消去法
前述的高斯消去法有消元和回代两过程,对消元过程稍加改变使方程组化为对角形方程组Dx=b形式,其中矩阵D为对角形矩阵
D=⎝⎜⎜⎜⎜⎛a11(1)a22(2)⋱ann(n)⎠⎟⎟⎟⎟⎞
此时求解无需回代
每次消元时利用主元将其所在列的其余元素全部消为0,即在第k次消元时,把k列中(k,k)位置以下和以上的元素消为0,但主行不进行消元
高斯——若尔当消去法执行一次归一化过程需要进行n-k+1次除法,而执行一次消去过程需要进行(n-1)(n-k+1)次乘法,因此高斯——若尔当消去法进行的乘除法总计算工作量为
k=1∑n(n−k+1)n=2n2(n+1)≈21n3
高斯——若尔当消去法比高斯消去法多了约1/6的计算工作量,但高斯——若尔当消去法无需回代,算法结构稍许简单
用归一化的高斯——若尔当消去法求矩阵的逆比较方便
设A非奇异,则A-1存在,记
A−1=X=⎝⎜⎜⎛x11⋮xn1⋯⋱⋯x1n⋮xnn⎠⎟⎟⎞
则求A的逆矩阵A-1等价于求n阶矩阵X,使AX=I,其中I为n×n的单位矩阵。将X和I分块
X=(x1,x2,⋯,xn),I=(e1,e2,⋯,en)
其中
xi=(x1i,x2i,⋯,xni)T
e1=(1,0,⋯,0)T,⋯,en=(0,⋯,0,1)T
则求解AX=I等价于求解系数矩阵相同的n个方程组Axj=ej,j=1,2,…,n,即
Ax1=⎝⎜⎜⎜⎜⎛10⋮0⎠⎟⎟⎟⎟⎞,Ax2=⎝⎜⎜⎜⎜⎛01⋮0⎠⎟⎟⎟⎟⎞,…,Axn=⎝⎜⎜⎜⎜⎛0⋮01⎠⎟⎟⎟⎟⎞
可用上述高斯——若尔当消去法(或再加上选主元技术)求解这n个方程组
三角分解法
LU
矩阵A各阶主子式不为零,则可唯一地分解成一个单位下三角阵L和一个非奇异的上三角阵U的乘积
各阶主子式不为零也可说成顺序主子阵非奇异。上述分解即是杜里特尔分解,对于克洛特分解所需条件也是一样的
顺序主子阵非奇异,因为非奇异矩阵就不一定存在LU分解,设
A=(0110)
则|A|=-1≠0,A非奇异。若A有LU分解,即存在数a,b,c,d,使
(0110)=(1a1)(bdc)
比较等式两边第一列,有b=0,ab=1
上两式不能同时成立,即A不存在LU分解
下三角形矩阵的乘积仍是下三角形矩阵
杜里特尔分解
⎝⎜⎜⎜⎜⎛a11a21⋮an1a12a22⋮an2⋯⋯⋱⋯a1na2n⋮ann⎠⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎛1l21⋮ln11⋮ln2⋱⋯1⎠⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎛u11u12u22⋯⋯⋱u1nu2n⋮unn⎠⎟⎟⎟⎟⎞
由矩阵乘法规则
a1i=u1i,i=1,2,…,n
ai1=li1u11,i=2,3,…,n
得U的第一行元素和L的第一列元素
u1i=a1i,i=1,2,…,n
li1=u11ai1,i=2,3,…,n
当得出U的前r-1行元素和L的前r-1列元素,,对于i=r,r+1,…,n,有
ari=k=1∑nlrkuki=k=1∑r−1lrkuki+uri
air=k=1∑nlikukr=k=1∑r−1likukr+lirurr
得计算U的第r行和L的第r列元素的公式
uri=ari−k=1∑r−1lrkuki,i=r,r+1,…,n
lir=(air−k=1∑r−1likukr)/urr,i=r+1,r+2,…,n
A的行列式|A|=|L||U|=u11u22…unn
解线性方程组
求解线性方程组Ax=b时,当对A进行LU分解时,等价于求解LUx=b,可归结为利用递推计算相继求解两个三角形(系数矩阵为三角矩阵)方程组,用顺代,由Ly=b,求出y,再利用回代,由Ux=y求出x
-
Ly=b,即计算
{y1=b1yi=bi−∑k=1i−1likyk,i=2,3,…,n
-
Ux=y,即计算
{xn=yn/unnxi=(yi−∑k=i+1nuikxk)/uii,i=n−1,…,2,1
用LU直接分解的方法求解线性方程组的计算量为1/3n3,和高斯消去法需要的计算量的数量级基本相同
求逆矩阵
- 对矩阵A和单位矩阵I组成的增广矩阵AI进行杜里特尔分解AI→LUY,其中A,I,L,U,Y∈Rn×n,L为单位下三角阵,U为上三角阵
- 对j=1,2,…,n求解方程组系Ux=yi,其中yi是矩阵Y的第j列。解记为x=aj,j=1,2,…,n
- A-1=(a1,a2,…,an)
克洛特分解
⎝⎜⎜⎛a11⋮an1⋯⋱⋯a1n⋮ann⎠⎟⎟⎞=⎝⎜⎜⎜⎜⎛l11l21⋮ln1l22⋮ln2⋱⋯lnn⎠⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎛1u121u13u23⋱⋯⋯u1nu2n⋮1⎠⎟⎟⎟⎟⎞
可得
li1=ai1,i=1,2,⋯,n
u1i=l11a1i,i=2,3,⋯,n
且当L的前r-1列元素和U的前r-1行元素已经算出以后,对于i=r,r+1,…,n有
lir=air−k=1∑r−1likukr,i=r,r+1,⋯,n
uri=(ari−k=1∑r−1lrkuki)/lrr,i=r+1,r+2,⋯,n
-
计算li1,i=1,2,…,n和uli,i=2,3,…,n
-
对于r=2,3,…,n计算lir,i=r,r+1,…,n和uri,i=r+1,r+2,…,n
-
求解Ly=b,即计算
{y1=l11b1yi=(bi−∑k=1i−1likyk)/lii,i=2,3,⋯,n
-
求解Ux=y,即计算
{xn=ynxi=yi−∑k=i+1nuikxk,i=n−1,⋯,2,1
平方根法
对称正定矩阵:矩阵A∈Rn×n满足AT=A,且对任意非零向量x∈Rn,有(Ax,x)=xTAx>0
- 对称正定矩阵A的对角元为正,即aii>0,i=1,2,…,n
- 实对称矩阵A正定的充要条件是A的所有特征值为正
- 对称正定矩阵非奇异,其逆亦为对称正定矩阵
- 实对称矩阵A正定的充要条件是A的所有顺序主子式为正
- 正定矩阵的顺序主子阵是正定的
乔累斯基分解
- 对称正定矩阵A存在唯一的单位下三角矩阵L和对角矩阵D,使A=LDLT
乔累斯基分解:对称正定矩阵A存在唯一的对角元素均为正数的下三角矩阵L,使A=LLT
乔累斯基分解所需要的乘除次数约为1/6n3数量级,差不多比LU分解节省一半的工作量,但要进行n次开方运算
将A=LLT写成
⎝⎜⎜⎜⎜⎜⎜⎛a11⋮an1⋱⋯aii⋯⋱a1n⋮ann⎠⎟⎟⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎛l11l21⋮ln1l22ln2⋱⋯lnn⎠⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎛l11l21l22⋯⋯⋱ln1ln2⋮ln⎠⎟⎟⎟⎟⎞
对于i=1,2,…,n,有
⎩⎪⎨⎪⎧lii=(ai−∑k=1i−1lik2)21lji=(aj−∑k=1i−1ljklik)/li,j=i+1,i+2,⋯,n
∣A∣=∣∣∣LLT∣∣∣=∣L∣∣∣∣LT∣∣∣=i=1∏nli2
将A=LLT代入Ax=b,有LLTx=b,令LTx=y,有
{Ly=bLTx=y
相应的求解公式
{yi=(bi−∑k=1i−1likyk)/lii,i=1,2,…,nxi=(yi−∑k=i+1nlkixk)/lii,i=n,n−1,…,1
由矩阵分解公式可知
aii=k=1∑ilik2
∣lik∣≤aii,k≤i
平方根法所求得的中间量lkk是可控的,即舍入误差的增长也可控
平方根法解正定方程组的缺点是需要进行开方运算
改进平方根法
将对称正定矩阵A=(aij)n×n进行LDLT分解,可避免开方运算,其中D=diag(di),且di>0,L为单位下三角矩阵,则
⎝⎜⎜⎜⎜⎜⎜⎛a11⋮an1⋱⋯aij⋯⋱a1n⋮ann⎠⎟⎟⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎛1l21⋮ln11⋮ln2⋱⋯1⎠⎟⎟⎟⎟⎞⎝⎜⎜⎜⎛d1d2⋱dn⎠⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎜⎛1l211l31l32⋱⋯⋯⋱1ln1ln2⋮lnj1⎠⎟⎟⎟⎟⎟⎟⎞
对=1,2,…,n有
{lij=(aij−∑k=1j−1likdkljk)/dj,di=aii−∑k=1i−1dklik2j=1,2,…,i−1
按此方式进行LDLT分解避免了开方运算,但在计算每个元时多了相乘的因子,为避免重复计算,进行如下变换:A=LDLT=TLT,其中T=LD,辅助变量tij=lijdj
代入得对于i=1,2,…,n计算
⎩⎪⎪⎨⎪⎪⎧tij=aij−∑k=1j−1tikljk,j=1,2,⋯,i−1lij=djtij,j=1,2,⋯,i−1di=aii−∑k=1i−1tiklik
Ax=b等价于求解L(DLTx)=b,可分解成由Ly=b求y,再由DLTx=y求x
{yi=bi−∑k=1i−1likyk,i=1,2,⋯,nxi=diyi−∑k=i+1nlkixk,i=n,n−1,⋯,1
范数
向量范数
∥x∥1=i=1∑n∣xi∣=∣x1∣+∣x2∣+⋯+∣xn∣
∥x∥2=(i=1∑nxi2)21=x12+x22+⋯+xn2
∥x∥∞=1≤i≤nmax∣xi∣=imax{∣x1∣,∣x2∣,⋯,∣xn∣}
矩阵范数
∥A∥∞=1≤i≤nmaxj=1∑n∣aij∣
∥A∥1=1≤j≤nmaxi=1∑n∣aij∣
∥A∥2=λmax(ATA)
其中λmax(ATA)表示矩阵ATA的最大特征值
假设我们有一个2×2的矩阵A:
A=(3113)
想要找到这个矩阵的特征值。首先,我们需要解特征方程,即det(A−λI)=0,其中I是单位矩阵,λ是特征值
矩阵A的特征方程是
det(3−λ113−λ)=0
(3−λ)(3−λ)−1⋅1=0
λ2−6λ+8=0
λ1=4,λ2=2
则矩阵A的最大特征值是λ1=4
∥A∥F=i=1∑mj=1∑naij2
迭代法
已知线性方程组Ax=b,其中A∈Rn×n,b∈Rn,A非奇异,设找到一个等价方程组x=Bx+f,从而建立迭代格式
x(k+1)=Bx(k)+f,k=0,1,⋯
当给定初始向量x(0)后,可按上式迭代,得到迭代向量序列{x(k)},当k→+∞时,x(k)→x*,则x*是线性方程组Ax=b的解
对线性方程组Ax=b用等价方程x=Bx+f建立迭代格式x(k+1)=Bx(k)+f,k=0,1,…逐步求解的方法叫迭代法。若limk→+∞x(k)=x∗,则称迭代法收敛,x*即方程组的解,否则称此迭代法发散
雅可比迭代
j=1∑naijxj=bi,i=1,2,⋯,n
设aii≠0,从上式中分离出变量xi,将它改写成
xi=aii1⎝⎛bi−j=1,j=i∑naijxj⎠⎞,i=1,2,⋯,n
据此建立得到雅可比迭代公式
xi(k+1)=aii1⎝⎛bi−j=1,j=i∑naijxj(k)⎠⎞,i=1,2,⋯,n
xi(k+1)=aii1(bi−j=1∑i−1aijxj(k)−j=i+1∑naijxj(k)),i=1,2,⋯,n
![雅可比迭代算法框图]()
迭代的收敛性
矩阵A∈Rn×n的所有特征值λi(i=1,2,…,n)的模的最大值称为矩阵A的谱半径ρ(A)
ρ(A)=1≤i≤nmax∣λi∣