一、高斯消元法
import numpy as np
def gauss(A, b):
'''
高斯消元法求解Ax=b，返回x向量
'''
m,n=A.shape[0],A.shape[1]
B=np.concatenate((A,b.reshape(m,1)),axis=1) # 构造增广矩阵
for k in range(m):
if B[k,k]==0:
for i in range(k+1,m):
if B[i,k]!=0:
B[[k,i],:]=B[[i,k],:]
for j in range(k+1,n):
B[j,:]=B[j,:]-(B[j,k]/B[k,k])*B[k,:]
x = np.empty(n)
for k in range(n-1,-1,-1):
x[k]=(B[k,n]-B[k,k+1:n].dot(x[k+1:]))/B[k,k]
return x
二、雅可比迭代法
关于几种迭代算法的迭代公式
def gauss_seidel(A, b):
'''
gauss_seidel迭代法求解Ax=b，返回x向量
'''
n=A.shape[0]
x=np.zeros(n)
res=np.linalg.norm(A.dot(x)-b)
k=1
while res>1e-5 and k<1000:
for i in range(n):
x[i]=x[i]+(b[i]-A[i,:].dot(x))/A[i,i]
k=k+1
res=np.linalg.norm(A.dot(x)-b)
return x
三、SOR迭代法
def sor(A, b, w):
'''
SOR迭代法求解Ax=b，返回x向量
'''
n=A.shape[0]
x=np.zeros(n)
res=np.linalg.norm(A.dot(x)-b)
k=1
while res>1e-5 and k<1000:
for i in range(n):
x[i]=x[i]+w*(b[i]-A[i,:].dot(x))/A[i,i]
k=k+1
res=np.linalg.norm(A.dot(x)-b)
return x
四、实例以及执行结果
A=np.array([[10,-1,-2],[-1,10,-2],[-1,-1,5]])
b=np.array([7.2,8.3,4.2])
print(gauss(A, b))
print(jacobi(A, b))
print(gauss_seidel(A, b))
print(sor(A, b, 0.5))
