数值分析实验2_求解线性方程组直接法

一 实验目的
1．掌握求解线性方程组的高斯消元法及列主元素法；
2. 掌握求解线性方程组的克劳特法；
3. 掌握求解线性方程组的平方根法。

二 实验内容
1．用高斯消元法求解方程组（精度要求为610-=ε）：
1231231
233272212240x x x x x x x x x -+=⎧⎪-+-=-⎨⎪-+=⎩ 2．用克劳特法求解上述方程组（精度要求为610-=ε）。

3. 用平方根法求解上述方程组（精度要求为610-=ε）。

4. 用列主元素法求解方程组（精度要求为610-=ε）：
1231231
233432222325x x x x x x x x x -+=⎧⎪-+-=⎨⎪--=-⎩ 三 实验步骤（算法）与结果
1.
程序代码（Python3.6）：
import numpy as np
def Gauss(A,b):
n=len(b)
for i in range(n-1):
if A[i,i]!=0:
for j in range(i+1,n):
m=-A[j,i]/A[i,i]
A[j,i:n]=A[j,i:n]+m*A[i,i:n]
b[j]=b[j]+m*b[i]
for k in range(n-1,-1,-1):
b[k]=(b[k]-sum(A[k,(k+1):n]*b[(k+1):n]))/A[k,k]
print(b)
运行函数：
>>> A=np.array([[3,-1,2],[-1,2,-2],[2,-2,4]],dtype=np.float) >>> b=np.array([7,-1,0],dtype=np.float)
>>> x=Gauss(A,b)
输出：
结果：解得原方程的解为x1=3.5,x2=-1,x3=-2.25
2
程序代码（Python3.6）：
import numpy as np
A=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(A[0])
i=0
while i<n:
j=i
while j<n:
U[i,j]=A[i,j]-sum(L[i,0:i]*U[0:i,j])
j+=1
k=i+1
while k<n:
L[k,i]=(A[k,i]-sum(L[k,0:i]*U[0:i,i]))/U[i,i]
k+=1
i+=1
print('L=',L)
print('U=',U)
def solvey(L,b):
n=len(y)
y[0]=b[0]
for i in range(1,n):
y[i]=b[i]-sum(L[i,0:i]*y[0:i])
print('y=',y)
def solvex(U,y):
n=len(x)
x[n-1]=y[n-1]/U[n-1,n-1]
for i in range(n-2,-1,-1):
x[i]=(y[i]-sum(U[i,(i+1):n]*x[(i+1):n]))/U[i,i]
print('x=',x)
运行函数：
>>> LU(A)
>>> solvey(L,b)
>>> solvex(U,y)
输出：
结果：同样解得原方程的解为x1=3.5,x2=-1,x3=-2.25
3程序代码（Python3.6）：
import numpy as np
A=np.array([[3,-1,2],[-1,2,-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,0],dtype=float)
def Cholesky(A):
n=len(A[0])
for k in range(n):
L[k,k]=pow(A[k,k]-sum(L[k,0:k]*L[k,0:k]),0.5)
for i in range(k+1,n):
L[i,k]=(A[i,k]-sum(L[i,0:i]*L[k,0:i]))/L[k,k]
print('L=',L)
def solvey(L,b):
n=len(y)
for i in range(n):
y[i]=(b[i]-sum(L[i,0:i]*y[0:i]))/L[i,i]
print('y=',y)
def solvex(L,y):
n=len(x)
for i in range(n-1,-1,-1):
x[i]=(y[i]-sum(L[(i+1):n,i]*x[(i+1):n]))/L[i,i]
print('x=',x)
运行函数：
>>> Cholesky(A)
>>> solvey(L,b)
>>> solvex(L,y)
输出：
结果：同样解得原方程的解为x1=3.5,x2=-1,x3=-2.25
4
程序代码（Python3.6）：
import numpy as np
A=np.array([[3,-1,4],[-1,2,-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=0
for i in range(k,n):
if abs(A[i,k])>A_max:
A_max=abs(A[i,k])
r=i
if A_max<1e-6:
print('系数矩阵奇异，无法求解方程！')
break
if r>k:
for j in range(k,n):
s=A[k,j]
A[k,j]=A[r,j]
A[r,j]=s
t=b[k]
b[k]=b[r]
b[r]=t
for j in range(k+1,n):
A[k,j]=A[k,j]/A[k,k]
b[k]=b[k]/A[k,k]
for i in range(n):
if i!=k:
for j in range(k+1,n):
A[i,j]=A[i,j]-A[i,k]*A[k,j]
b[i]=b[i]-A[i,k]*b[k]
print('x=',b)
运行函数：
>>> Main_Gauss(A,b)
输出：
结果：解得原方程的解为x1=1,x2=2,x3=0.5
四实验收获与教师评语
实验收获：掌握了求解线性方程组的高斯消元法、列主元素法、克劳
特法和平方根法等算法流程，以及能够运用Python、MA TLAB等语言实现并解出方程组的根。