數(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]是否為絕對(duì)值最大的元素，如果不是，則交換兩行，

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每行獲取主對(duì)角線元素，并通過(guò)行初等變換將其下的各行首個(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]是否為絕對(duì)值最大的元素，如果不是，則交換兩行，

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。例如針對(duì)例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中的方程組得到了求解。對(duì)三對(duì)角矩陣可以直接應(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例

用追趕法解下面三對(duì)角方程組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）后程序沒(méi)有得到結(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_method(A,np.array([[1.],[1.],[1.]]),1e-5)第1次循環(huán)時(shí),m_k=8.00000000第2次循環(huán)時(shí),m_k=9.25000000第3次循環(huán)時(shí),m_k=9.54054054第4次循環(huán)時(shí),m_k=9.59490085第5次循環(huán)時(shí),m_k=9.60407440第6次循環(huán)時(shí),m_k=9.60542900第7次循環(huán)時(shí),m_k=9.60557200第8次循環(huán)時(shí),m_k=9.60556710冪法與反冪法的實(shí)現(xiàn)、基于QR分解的特征值求法defhouseholder(w):

n=len(w)

w=w.copy().reshape((n,1))

m=la.norm(w,2)#盡量使新w[0]的絕對(duì)值變大，于是按原w[0]的符號(hào)計(jì)算新w[0]

ifw[0]<0:

w[0]=w[0]-m

else:

w[0]=