数值分析实验2_求解线性方程组直接法
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
一 实验目的
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 j=i while j U[i,j]=A[i,j]-sum(L[i,0:i]*U[0:i,j]) j+=1 k=i+1 while k 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)