數(shù)值分析實驗2_求解線性方程組直接法

1、一實驗?zāi)康? .掌握求解線性方程組的高斯消元法及列主元素法；2 .掌握求解線性方程組的克勞特法；3 .掌握求解線性方程組的平方根法。二實驗內(nèi)容1 .用高斯消元法求解方程組(精度要求為8 = 10)：3x1 - x2 2x3 7 7-x1 2x2 - 2x3 - 1 2x1 -2x2 4x3 = 02 .用克勞特法求解上述方程組(精度要求為w = 10)3 .用平方根法求解上述方程組(精度要求為 君=10)4 .用列主元素法求解方程組(精度要求為名= 10屆)：3x1 -x2 4x3 = 3-x1 2x2 -2x3 =2 2 x1 - 3 x2 - 2x3 = -5三 實驗步驟(算法)與結(jié)果1.

2、程序代碼(Python3.6 ):import numpy as npdef Gauss(A,b):n=len(b)for i in range(n-1):if Ai,i!=0:for j in range(i+1,n):m=-Aj,i/Ai,iAj,i:n=Aj,i:n+m*Ai,i:n bj=bj+m*bi for k in range(n-1,-1,-1):bk=(bk-sum(Ak,(k+1):n*b(k+1):n)/Ak,kprint(b)運行函數(shù):>>> A=np.array(3,-1,2,-1,2,-2,2,-2,4,dtype=np.float)>>

3、;> b=np.array(7,-1,0,dtype=np.float)>>> x=Gauss(A,b)輸出：»> _= RESTART: £；*數(shù)值分析9=t仃=> >> A=np. array (3, - 1T 2, -1,2,-2, 2r-2, 4 J, d-type=np. f loat.)> » b-np. array (7,-1,0, dtype-np. float)> >> x=Gauss b!3.5 -L -2,25結(jié)果：解得原方程的解為x1=3.5,x2=-1,x3=-2.2

4、52程序代碼(Python3.6 ):import numpy as npA=np.array(3,-1,2,-1,2,-2,2,-2,4,dtype=float) L=np.array(1,0,0,0,1,0,0,0,1,dtype=float) U=np.array(0,0,0,0,0,0,0,0,0,dtype=float) b=np.array(7,-1,0,dtype=float) y=np.array(0,0,0,dtype=float) x=np.array(0,0,0,dtype=float)def LU(A):n=len(A0) i=0 while i<n: j=i w

5、hile j<n:UiJ=AiJ-sum(Li,0:i*U0:i,j) j+=1 k=i+1 while k<n:Lk,i=(Ak,i-sum(Lk,0:i*U0:i,i)/Ui,i k+=1 i+=1print('L=',L)print('U=',U)def solvey(L,b):n=len(y)y0=bfor i in range(1,n):yi=bi-sum(Li,0:i*y0:i) print('y=',y)def solvex(U,y):n=len(x)xn-1=yn-1/Un-1,n-1for i in range(n-

6、2,-1,-1):xi=(yi-sum(Ui,(i+1):n*x(i+1):n)/Ui,i print('x=',x)運行函數(shù)：>>> LU(A)>>> solvey(L,b)>>> solvex(U,y)輸出：»>»> LUCA)=RESTART: E:，數(shù)值分析Vtest.pvL= LD.d1-0,33333333L0,0,66666667r(L 8L1U= 3.-1.2,o.o.1.66666667 0.-1. 33333333：L61>> solvey (Li b)y= 7

7、,L 33333333 -3. 6>» salves (UT y)x= 3.5 -L -24 25»>結(jié)果：同樣解得原方程的解為 x1=3.5,x2=-1,x3=-2.253 程序代碼(Python3.6 ):import numpy as npA=np.array(3,-1,2,-12-2,2,-2,4,dtype=float) L=np.array(0,0,0,0,0,0,0,0,0,dtype=float) b=np.array(7,-1,0,dtype=float) y=np.array(0,0,0,dtype=float) x=np.array(0,0

8、,0,dtype=float)def Cholesky(A):n=len(A0)for k in range(n):Lk,k=pow(Ak,k-sum(Lk,0:k*Lk,0:k),0.5)for i in range(k+1,n):Li,k=(Ai,k-sum(Li,0:i*Lk,0:i)/Lk,kprint('L=',L)def solvey(L,b):n=len(y)for i in range(n):yi=(bi-sum(Li,0:i*y0:i)/Li,iprint('y=',y)def solvex(L,y):n=len(x)for i in rang

9、e(n-1,-1,-1):xi=(yi-sum(L(i+1):n,i*x(i+1):n)/Li,iprint('x=',x)運行函數(shù)：>>> Cholesky(A)>>> solvey(L,b)>>> solvex(L,y)輸出：»>= = RESTART: E，數(shù)值分析emt*py => >> Cholesky(A)L= L 73205081 0.0.-0.57735027 1.29099445 0,L 15470054 -L 03279556 1.26491106> »

10、solvey (L, b)y= 4,04145198 1.03279556 -2.34604989> » solvex (Li y)x= 3.5 -L -2. 25»>結(jié)果：同樣解得原方程的解為x1=3.5,x2=-1,x3=-2.254程序代碼(Python3.6 ):import numpy as npA=np.array(3,-1,4,-12-2,2,-3,-2,dtype=float) b=np.array(3,2,-5,dtype=float)def Main_Gauss(A,b):n=len(b)for k in range(n):A_max=0fo

11、r i in range(k,n):if abs(Ai,k)>A_max: A_max=abs(Ai,k) r=iif A_max<1e-6:print('系數(shù)矩陣奇異，無法求解方程！)breakif r>k:for j in range(k,n):s=Ak,jAk,j=Ar,jAr,j=st=bkbk=brbr=tfor j in range(k+1,n):Ak,j=Ak,j/Ak,kbk=bk/Ak,kfor i in range(n):if i!=k:for j in range(k+1,n):Ai,j=Ai,j-Ai,k*Ak,j bi=bi-Ai,k*bk print('x=',b)運行函數(shù)：>>> Main_Gauss(A,b)輸出：»>=-=:=-=-= RESTART: E: 薊 值分析t