数值分析实验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

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)