數(shù)值分析方法面向“四新”人才培養(yǎng)普通高等教育系列教材第二章

數(shù)值代數(shù)基礎(chǔ)2.1線性方程組的直接解法2.2向量與矩陣的范數(shù)2.3線性方程組的迭代解法2.4矩陣特征值計(jì)算2.5

Python程序在數(shù)值代數(shù)中的應(yīng)用

目錄/Contents

2.1線性方程組的直接解法

記2.1.1Gauss消去法消元過程

Gauss消去法的消元過程由n-1步完成第一步：則消元過程

Gauss消去法的消元過程由n-1步完成第二步：則消元過程

Gauss消去法的消元過程由n-1步完成重復(fù)上述步驟直至n-1步回代求解Gauss消去法總的計(jì)算量為例用Gauss消去法求解線性方程組解

Gauss消去法的消元過程是將方程組的增廣矩陣進(jìn)行下列一系列的初等行變換由此得到與方程組同解的上三角形方程組回代求解，得適用條件定理：若A的各階順序主子式不等于0，則高斯消去法能順序進(jìn)行消元，得到唯一解。存在問題如果某個(gè),但很小，會(huì)引入較大的誤差。1.在高斯消去法消去過程中可能出現(xiàn)的情況，這時(shí)高斯消去法將無法進(jìn)行2.1.2Gauss列主元素消去法基本思想在每輪消元之前，選列主元素（絕對值最大的元素）具體步驟如此至多經(jīng)過n-1步,就得到與之同解的上三角形方程組的增廣矩陣,再用回代過程即可得方程組的解.Gauss全主元素消去法例用Gauss列主元消元法和全主元消去法解線性方程組解

（1）列主元素消去法第一步，選擇-20作為該列的主元素，（2）全主元素消去法在所有系數(shù)中選擇絕對值最大的40作為主元素，交換第一、二行和交換第一、二列使該主元素位于對角線的第一個(gè)位置上，得2.1.3矩陣的三角分解法若矩陣A有分解：A=LU，其中L為下三角陣，U為上三角陣，則稱該分解為A的LU分解，又稱為Doolittle分解若矩陣A有分解A=LU，則解線性方程組Ax=b就等價(jià)于求解定理利用矩陣的乘法運(yùn)算介紹Doolittle分解法求解線性方程組的具體計(jì)算過程。

設(shè)A=LU由矩陣乘法規(guī)則得由此可得U的第1行元素和L的第1列元素:

①

計(jì)算U的第k行元素從而有②計(jì)算L的第k列元素，由于于是由得自此就可以利用上述公式逐步求出U與L的各元素，因此Gauss消元法解線性方程組就等價(jià)于解兩個(gè)三角方程組

例

用Doolittle分解法求解方程組解：（1）用分解公式（2.1.8）~（2.1.11）計(jì)算A=LU，令（2）解方程Ly=b（3）解方程Ux=y2.1.4對稱矩陣的楚列斯基分解（平方根法）

，i=1,2,...,n;j=i+1,...,n對稱矩陣的楚列斯基分解（平方根法）

其中D為對角陣，U0為上三角陣，于是有A=L1U=L1DU0，又A為對稱陣，A=AT=(

L1DU0)T=U0T

DT

L1T=

U0T

D

L1T，有分解唯一性得U0T

=L1，即A=

L1DL1T.因?yàn)锳是正定陣，A的各階順序主子式為正，故于是對角陣D可以分解為利用待定系數(shù)法來計(jì)算Cholesky分解中下三角陣L的元素lij.由根據(jù)矩陣乘法運(yùn)算，比較上式兩邊元素：由，因?yàn)锳是正定陣，設(shè)A為n

階對稱正定陣，由定理可知，A可以唯一的分解為A=LDLT

其中矩陣L為單位下三角陣，D為非奇異陣.設(shè)A=LDLT

為了減少乘除運(yùn)算次數(shù)，刪除重復(fù)計(jì)算，令

，則

例

用Cholesky分解法求解方程組解：利用公式對系數(shù)矩陣進(jìn)行Cholesky分解平方根法與改進(jìn)的平方根法不僅計(jì)算量是Gauss消元法的一半，其數(shù)值穩(wěn)定性良好，是求解中小型稠密對稱正定線性方程組的好方法.2.1.5解三對角線性方程組的追趕法該類方程組中除了對角線和兩條相鄰的對角線外，其他元素都為零，稱該類方程組為三對角方程組解三對角線方程組的追趕法，其計(jì)算公式如下：第一步

對三對角陣A做Doolittle分解，則有第二步解方程組Ly=f（“追”的過程），得（2.1.20）（2.1.21）第三步解方程組Ux=y（“趕”的過程）得追趕法的基本思想與Gauss消元法及三角分解法相同，只是由于系數(shù)矩陣的特殊性，使得求解的計(jì)算公式簡化，計(jì)算量減少，僅為5n-4次乘除法.此法是一種具有最優(yōu)運(yùn)算量且具有良好的數(shù)值穩(wěn)定性.（2.1.22）例

用追趕法解下面三對角方程組解：由公式(2.1.20)得于是系數(shù)矩陣分解為由公式（2.1.21）解方程組Ly=f

得由公式（2.1.22）解方程組Ux=y

得3、遞推化：把復(fù)雜的計(jì)算歸結(jié)為簡單過程的多次重復(fù)計(jì)算，易于用循環(huán)結(jié)構(gòu)來實(shí)現(xiàn)（如迭代法）謝謝！數(shù)值分析方法面向“四新”人才培養(yǎng)普通高等教育系列教材主編

李冬果李林高磊首都醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院智能醫(yī)學(xué)工程學(xué)學(xué)系第二章

數(shù)值代數(shù)基礎(chǔ)2.1線性方程組的直接解法2.2向量與矩陣的范數(shù)2.3線性方程組的迭代解法2.4矩陣特征值計(jì)算2.5

Python程序在數(shù)值代數(shù)中的應(yīng)用

目錄/Contents

2.2向量和矩陣的范數(shù)2.2.1向量范數(shù)向量范數(shù)性質(zhì)2.2.2矩陣范數(shù)

證明：要證明（2.2.2）是一種矩陣范數(shù)，需要證明其滿足定義2.2中的非負(fù)性、齊次性、三角不等式、和相容條件。

2.2.3病態(tài)方程組與矩陣的條件數(shù)

這類方程組的解對右端常數(shù)項(xiàng)的擾動(dòng)很敏感，我們把這類方程組稱為病態(tài)方程組

考察常數(shù)項(xiàng)b的微小誤差對解的影響

設(shè)A是精確的，

b是有誤差（或擾動(dòng))的，誤差或擾動(dòng)記為δb.

，則

上述結(jié)論說明矩陣A的條件數(shù)越大，解的相對誤差

也越大，此時(shí)方程組的病態(tài)程度也越高n階Hilbert矩陣3、遞推化：把復(fù)雜的計(jì)算歸結(jié)為簡單過程的多次重復(fù)計(jì)算，易于用循環(huán)結(jié)構(gòu)來實(shí)現(xiàn)（如迭代法）謝謝！數(shù)值分析方法面向“四新”人才培養(yǎng)普通高等教育系列教材主編

李冬果李林高磊首都醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院智能醫(yī)學(xué)工程學(xué)學(xué)系第二章

數(shù)值代數(shù)基礎(chǔ)2.1線性方程組的直接解法2.2向量與矩陣的范數(shù)2.3線性方程組的迭代解法2.4矩陣特征值計(jì)算2.5

Python程序在數(shù)值代數(shù)中的應(yīng)用

目錄/Contents

2.3線性方程組的迭代解法2.3.1迭代法的基本思想（2.3.1）（2.3.2）2.3.2迭代法的收斂條件事后估計(jì)事前估計(jì)對線性方程組構(gòu)造迭代格式

2.3.3Jacobi迭代法據(jù)此建立迭代格式為:即為解方程組的Jacobi迭代公式Jacobi迭代法的基本思想對線性方程組

若（2.3.5）Jacobi迭代法的的矩陣表示

A=L+D+U

據(jù)此得到迭代公式為式(2.3.6)稱為Jacobi迭代公式，B稱為Jacobi迭代矩陣（2.3.6）

例

用Jacobi迭代法求解方程組解：將原方程組化為等價(jià)方程組由迭代公式得Jacobi迭代格式為其對應(yīng)的Jacobi迭代矩陣為每次迭代充分利用當(dāng)前最新的迭代值，即在求第k+1次迭代時(shí)的第i個(gè)分量xi(k+1)，用新分量x1(k+1),x2(k+1),...,xi-1(k+1)代替舊分量x1(k),x2(k),...,xi-1(k)，在整個(gè)迭代過程中只要使用一組單元存放迭代向量，這樣建立起來的迭代格式稱為Gauss-Seidel迭代.其迭代法格式如下：2.3.4Gauss-Seidel迭代法（2.3.9）

Gauss-Seidel迭代法的矩陣表示Gauss-Seidel迭代法的收斂性定理2.13：解線性方程組的Gauss-Seidel迭代法收斂的充分必要條件為其迭代矩陣的譜半徑小于1例

用Gauss-Seidel迭代法求解方程組解

由迭代公式（2.3.9）得Gauss-Seidel迭代格式為：其對應(yīng)的Gauss-Seidel迭代矩陣為例

給定線性方程組討論采用Jacobi迭代法和Gauss-Seidel迭代法求解時(shí)的收斂性，該方程組精確解為解

將原方程組化為等價(jià)方程組（1）Jacobi迭代法的迭代格式為（2）Gauss-Seidel迭代法的迭代格式為定理2.14：線性方程組Ax=b的系數(shù)矩陣A為嚴(yán)格對角占優(yōu)矩陣，則解此方程組的Jacobi法和Gauss-Seidel迭代法都是收斂的.證明

由于系數(shù)矩陣A為嚴(yán)格對角占優(yōu)矩陣，所以A為非奇異矩陣，Jacobi迭代公式的迭代矩陣為利用對角占優(yōu)知由定理2.11知Jacobi迭代法收斂.2.3.5超松弛迭代法（SOR方法）超松弛迭代法(SuccessiveOverRelaxationMethod)是Gauss-seidel迭代法的改進(jìn)，是解大型稀疏矩陣線性方程組的有效方法。

（2.3.15）（2.3.16）將（2.3.15）代入（2.3.16）化簡得超松弛迭代公式

SOR迭代法的矩陣表示例

用SOR迭代法解線性方程組

解：SOR法的迭代格式為3、遞推化：把復(fù)雜的計(jì)算歸結(jié)為簡單過程的多次重復(fù)計(jì)算，易于用循環(huán)結(jié)構(gòu)來實(shí)現(xiàn)（如迭代法）謝謝！數(shù)值分析方法面向“四新”人才培養(yǎng)普通高等教育系列教材主編

李冬果李林高磊首都醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院智能醫(yī)學(xué)工程學(xué)學(xué)系第二章

數(shù)值代數(shù)基礎(chǔ)2.1線性方程組的直接解法2.2向量與矩陣的范數(shù)2.3線性方程組的迭代解法2.4矩陣特征值計(jì)算2.5

Python程序在數(shù)值代數(shù)中的應(yīng)用

目錄/Contents

2.4矩陣特征值計(jì)算2.4.1冪法與反冪法冪法是求矩陣按模最大特征值（主特征值）及相應(yīng)特征向量的簡單而有效的方法，特別適用于大型矩陣或稀疏矩陣，也是計(jì)算譜半徑的有效方法。

則

例

用冪法求解矩陣的按模最大特征值及其相應(yīng)的特征向量，精確值

加速收斂：引進(jìn)矩陣B

反冪法是冪法的變形，是用來計(jì)算A的按模最小的特征值和相應(yīng)的特征向量。則

在反冪法中也可用原點(diǎn)位移來加速迭代或求其他特征值例

求解矩陣A的給定近似特征值p=-13的特征值及其相應(yīng)的特征向量，其中2.4.2基于Householder變換的QR分解

例

試給出一個(gè)Householder矩陣QR分解是工程中應(yīng)用最為廣泛的一類矩陣分解，它是將一個(gè)矩陣分解為正交矩陣Q和上三角矩陣R的乘積第一步：對n階實(shí)方陣A的第一列進(jìn)行Householder變換，取矩陣的QR分解

第二步：

QR算法

如此繼續(xù)下去，可得

3、遞推化：把復(fù)雜的計(jì)算歸結(jié)為簡單過程的多次重復(fù)計(jì)算，易于用循環(huán)結(jié)構(gòu)來實(shí)現(xiàn)（如迭代法）謝謝！數(shù)值分析方法面向“四新”人才培養(yǎng)普通高等教育系列教材主編

李冬果李林高磊首都醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院智能醫(yī)學(xué)工程學(xué)學(xué)系第二章

數(shù)值代數(shù)基礎(chǔ)2.1線性方程組的直接解法2.2向量與矩陣的范數(shù)2.3線性方程組的迭代解法2.4矩陣特征值計(jì)算2.5

Python程序在數(shù)值代數(shù)中的應(yīng)用

目錄/Contents

2.5Python程序在數(shù)值代數(shù)中的應(yīng)用線性方程組的直接解法的實(shí)現(xiàn)Gauss消去法、矩陣的三角分解法的實(shí)現(xiàn)（Doolittle/Cholesky分解）線性方程組的迭代解法的實(shí)現(xiàn)Jacobi迭代法、Gauss-Seidel迭代法、超松弛迭代法(SOR)法的實(shí)現(xiàn)矩陣特征值的Python計(jì)算冪法與反冪法的實(shí)現(xiàn)、基于QR分解的特征值求法SciPy工具包linalg模塊（la.inv，la.solve，la.det，la.norm，la.eig，la.lu,la.svd,la.qr）Gauss消去法defgauss_elimination(A,b):

A_aug=np.c_[A,b].astype(np.float64)

rows,cols=A_aug.shape

foriinrange(rows-1):

#判斷A_aug[i,i]是否為絕對值最大的元素，如果不是，則交換兩行，

max_index=np.argmax(np.abs(A_aug[i:,i]))

ifmax_index!=0:

max_index=max_index+i

A_aug[[i,max_index],:]=A_aug[[max_index,i],:]

m0=A_aug[i,i]

forjinrange(i+1,rows):

m1=-A_aug[j,i]/m0

forkinrange(i,cols):

A_aug[j,k]=m1*A_aug[i,k]+A_aug[j,k]

returnA_aug函數(shù)gauss_elimination_0會(huì)依次在A_aug每行獲取主對角線元素，并通過行初等變換將其下的各行首個(gè)非0系數(shù)修改為0，從而得到最后的上三角矩陣。例如A=np.array([[2.,-1.,3.],[4.,2.,5.],[1.,2.,2.]])b=np.array([[1.],[1.],[1.]])B=gauss_elimination_0(A,b)結(jié)果為：array([[2.

,-1.

,

3.

,

1.

],

[0.

,

4.

,-1.

,-1.

],

[0.

,

0.

,

1.125,

1.125]])>>>Gauss消去法back_substitution(B)array([[-1.],[0.],[1.]])defback_substitution(A_aug):

rows,cols=A_aug.shape

A=A_aug[:,:-1].copy()

b=A_aug[:,-1].copy().reshape((rows,1))

res=np.zeros((rows,1),np.float64)

foriinrange(rows-1,-1,-1):

s=0

forjinrange(i+1,rows):

s+=A[i,j]*res[j]

res[i]=(b[i]-s)/A[i,i]

returnres在Gauss消去法的基礎(chǔ)上增加一個(gè)搜索列主值的步驟就得到了列主元素法defgauss_elimination(A,b):

A_aug=np.c_[A,b].astype(np.float64)

rows,cols=A_aug.shape

foriinrange(rows-1):

#判斷A_aug[i,i]是否為絕對值最大的元素，如果不是，則交換兩行，

max_index=np.argmax(np.abs(A_aug[i:,i]))

ifmax_index!=0:

max_index=max_index+i

A_aug[[i,max_index],:]=A_aug[[max_index,i],:]

m0=A_aug[i,i]

forjinrange(i+1,rows):

m1=-A_aug[j,i]/m0

forkinrange(i,cols):

A_aug[j,k]=m1*A_aug[i,k]+A_aug[j,k]

returnA_aug回代部分不需要修改。練習(xí)：用該函數(shù)解線性方程組A=np.array([[10.,-19.,-2.],

[-20.,40.,1],

[1.,4.,5.]],dtype=np.float64)b=np.array([3.,4.,5.]

,dtype=np.float64).reshape((3,1))B=gauss_elimination(A,b)結(jié)果為：>>>print(B)[[-20.

40.

1.

4.

]

[

0.

6.

5.05

5.2

]

[

0.

0.

-2.34166667

4.13333333]]回代后有：>>>back_substitution(B)array([[4.41637011],

[2.35231317],

[-1.76512456]])Doolittle分解defdoolittle(A):

n=A.shape[0]

L=np.zeros((n,n),np.float64)

U=np.zeros((n,n),np.float64)

foriinrange(n):

forkinrange(i,n):

s1=0

#求

L(i,j)*U(j,k)的和

forjinrange(i):

s1+=L[i,j]*U[j,k]

U[i,k]=A[i,k]-s1

forkinrange(i,n):

ifi==k:

L[i,i]=1

else:

s2=0#求

L(k,j)*U(j,i)的和

forjinrange(i):

s2+=L[k,j]*U[j,i]

L[k,i]=(A[k,i]-s2)/U[i,i]

returnL,U函數(shù)doolittle可以把系數(shù)矩陣A分解為單位下三角陣L和非奇異上三角陣U。例如針對例2.1.4可以使用以下doolittle函數(shù)進(jìn)行分解：>>>A=np.array([[2,5,-6],[4,13,-19],[-6,-3,-6]],dtype=np.float64)>>>L,U=doolittle(A)>>>print(L)[[1.

0.

0.]

[2.

1.

0.]

[-3.

4.

1.]]>>>print(U)[[2.

5.-6.]

[0.

3.-7.]

[0.

0.

4.]]defback_sub_for_dool(L,U,b):

rows,cols=L.shape

y=np.zeros((rows,1),np.float64)

foriinrange(rows):

s=0

forjinrange(i):

s+=L[i,j]*y[j]

y[i]=(b[i]-s)

x=np.zeros((rows,1),np.float64)

foriinrange(rows-1,-1,-1):

s=0

forjinrange(i+1,rows):

s+=U[i,j]*x[j]

x[i]=(y[i]-s)/U[i,i]

returnx,yDoolittle分解：回代函數(shù)則繼續(xù)上例：>>>b=np.array([[10],[19],[-30]],dtype=np.float64)>>>x,y=back_sub_for_dool(L,U,b)>>>xarray([[3.],

[2.],

[1.]])>>>yarray([[10.],

[-1.],

[4.]])于是例2.1.4中的方程組得到了求解。對三對角矩陣可以直接應(yīng)用Doolittle分解，也可以編程實(shí)現(xiàn)公式2.1.20-2.1.22：deftridiagmatrixalg(A,f):

rows,cols=A.shape

u=np.zeros((rows,),dtype=np.float64)

l=np.zeros((rows-1,),dtype=np.float64)

a=[A[i+1,i]foriinrange(rows-1)]

b=[A[i,i]foriinrange(rows)]

c=[A[i,i+1]foriinrange(rows-1)]

u[0]=b[0]

foriinrange(1,rows):

l[i-1]=a[i-1]/u[i-1]

u[i]=b[i]-l[i-1]*c[i-1]

y=np.zeros((rows,1),dtype=np.float64)

y[0]=f[0]

foriinrange(1,rows):

y[i]=f[i]-l[i-1]*y[i-1]

x=np.zeros((rows,1),dtype=np.float64)

x[-1]=y[-1]/u[-1]

foriinrange(-2,-rows-1,-1):

x[i]=(y[i]-c[i+1]*x[i+1])/u[i]

returnx,y,l,u例

用追趕法解下面三對角方程組A=np.array([[2,1,0,0],

[1,3,1,0],

[0,1,1,1],

[0,0,2,1]],dtype=np.float64)b=np.array([1,2,2,0],dtype=np.float64).reshape((4,1))x,y,l,u=tridiagmatrixalg(A,b)print(x)Jacobi迭代法defsolver_jacobi(A,b,ep=1e-8,M=100):

k=0

n=len(b)

x=np.zeros((n,1),dtype=np.float64())

y=x.copy()

whilek<M:

foriinrange(n):

s=0

forjinrange(i):

s+=A[i][j]*x[j]

forjinrange(i+1,n):

s+=A[i][j]*x[j]

y[i]=(b[i]-s)/A[i][i]

ifmax(np.abs(x-y))<ep:

returny,k+1

k+=1

x=y.copy()

print(f"stopedwheniterationisgreaterthan{M:d}")例

用Jacobi迭代法求解方程組>>>A=np.array([[1,2,-2],...

[1,1,1],...

[2,2,1]],dtype=np.float64)>>>b=np.array([[6],[6],[11]],dtype=np.float64)>>>x,k=solver_jacobi(A,b)>>>print(x)[[2.]

[3.]

[1.]]>>>print(k)4Gauss-Seidel迭代法defsolver_gs(A,b,ep=1e-8,M=100):

k=0

n=len(b)

x=np.zeros((n,1),dtype=np.float64())

y=x.copy()

whilek<M:

foriinrange(n):

s=0

forjinrange(i):

s+=A[i][j]*y[j]

forjinrange(i+1,n):

s+=A[i][j]*x[j]

y[i]=(b[i]-s)/A[i][i]

ifmax(np.abs(x-y))<ep:

returny,k+1

k+=1

x=y.copy()

print(f"stopedwheniterationisgreaterthan{M:d}")例

用Jacobi迭代法求解方程組>>>A=np.array([[1,2,-2],...

[1,1,1],...

[2,2,1]],dtype=np.float64)>>>b=np.array([[6],[6],[11]],dtype=np.float64)>>>solver_gs(A,b)stopedwheniterationisgreaterthan100在循環(huán)達(dá)到了默認(rèn)的循環(huán)上限（100）后程序沒有得到結(jié)果，由2.3.4節(jié)的討論知，此時(shí)Gauss-Seidel迭代法發(fā)散。例

給定線性方程組討論采用Jacobi迭代法和Gauss-Seidel迭代法求解時(shí)的收>>>A=np.array([[8,-1,1],...

[2,10,-1],...

[1,1,-5]],dtype=np.float64)>>>b=np.array([[8],[11],[-3]],dtype=np.float64)斂性，該方程組精確解為>>>solver_jacobi(A,b,ep=1e-4)(array([[1.0000052],

[1.00000423],

[0.99999586]]),8)>>>solver_gs(A,b,ep=1e-4)(array([[1.00000141],

[0.99999922],

[1.00000013]]),5)超松弛迭代法(SOR)法defsolver_sor(A,b,w=1.5,ep=1e-8,M=100):

k=0

n=len(b)

x=np.zeros((n,1),dtype=np.float64())

y=x.copy()

whilek<M:

foriinrange(n):

s=0

forjinrange(i):

s+=A[i][j]*y[j]

forjinrange(i+1,n):

s+=A[i][j]*x[j]

xk=(b[i]-s)/A[i][i]

y[i]=(1-w)*x[i]+w*xk

ifmax(np.abs(x-y))<ep:

returny,k+1

k+=1

x=y.copy()

print(f"stopedwheniterationisgreaterthan{M:d}")例

用SOR迭代法解線性方程組>>>A=np.array([[4,-2,-4],...

[-2,17,10],...

[-4,10,9]],dtype=np.float64)>>>b=np.array([[10],[3],[-7]],dtype=np.float64)>>>solver_sor(A,b,w=1.46,ep=1e-6)(array([[1.9999978],

[1.00000144],

[-1.0000026]]),21)冪法與反冪法的實(shí)現(xiàn)defpow_method(A,u0,ep):

k=1

m0=np.max(abs(u0))

whileTrue:

v=u0/m0

u0=A@v

mk=np.max(np.abs(u0))

print(f"第{k:2d}次循環(huán)時(shí),m_k={mk:10.8f}")

ifnp.abs(mk-m0)<ep:

return(mk,v)

m0=mk

k+=1例

用冪法求解矩陣的按模最大特征值及其相應(yīng)的特征向量，精確值>>>A=np.array([[7,3,-2],...

[3,4,-1],...

[-2,-1,3]],dtype=np.float64)>>>l,v=pow_met