定义矩阵乘法
def mult(h, x):
result = []
for x in h:
summ = 0
for index, y in enumerate(x):
summ += y * x[index]
result.append(summ)
return result
创建希尔伯特矩阵
def create_hobert(n):
h=[]
for x in range(1, n + 1):
row = []
for y in range(1, n + 1):
row.append(1 / (x + y - 1))
h.append(row)
return h
雅克比迭代法
def jacobi(A,B,sigma,N):
n = len(A)
x0 = []
x = []
for i in range(0,n):
x0.append(0)
x.append(0)
for k in range(1,N+1):
R = 0
for i in range(0,n):
sum_ax = 0
for j in range(0,n):
sum_ax = sum_ax + A[i][j] * x0[j]
x[i] = x0[i] + ((B[i] - sum_ax)/A[i][i])
if abs(x[i] - x0[i]) > R:
R = abs(x[i] - x0[i])
if R <= sigma:
print("精确度等于",sigma,"时,jacobi迭代法需要迭代",k,"次收敛")
return (x,k)
for i in range(0,n):
x0[i] = x[i]
return (x,k)
Hauss-Seidel迭代法
def gauss_seidel(A,B,sigma,N):
n = len(A)
x0 = []
x = []
for i in range(0,n):
x0.append(0)
x.append(0)
for k in range(1,N+1):
R = 0
for i in range(0,n):
sum_ax = 0
for j in range(0,n):
if j >= i:
sum_ax = sum_ax + A[i][j] * x0[j]
else:
sum_ax = sum_ax + A[i][j] * x[j]
x[i] = x0[i] + ((B[i] - sum_ax)/A[i][i])
if abs(x[i] - x0[i]) > R:
R = abs(x[i] - x0[i])
if R <= sigma:
print("精确度等于",sigma,"时,gauss_seidel迭代法需要迭代",k,"次收敛")
return (x,k)
for i in range(0,n):
x0[i] = x[i]
return (x,k)
SOR迭代法
def sor(A,B,sigma,N,omega):
n = len(A)
x0 = []
x = []
for i in range(0,n):
x0.append(0)
x.append(0)
for k in range(1,N+1):
R = 0
for i in range(0,n):
sum_ax = 0
for j in range(0,n):
if j >= i:
sum_ax = sum_ax + A[i][j] * x0[j]
else:
sum_ax = sum_ax + A[i][j] * x[j]
x[i] = x0[i] + omega * ((B[i] - sum_ax)/A[i][i])
if abs(x[i] - x0[i]) > R:
R = abs(x[i] - x0[i])
if R <= sigma:
print("精确度等于",sigma,"时,松弛因子为",omega,"时,sor迭代法需要迭代",k,"次收敛")
print(x)
return (x,k)
for i in range(0,n):
x0[i] = x[i]
return (x,k)
测试希尔伯特矩阵H
if __name__ == "__main__":
n = 10
H = create_hobert(n)
x = [1 for x in range(n)]
b = mult(H,x)
print('雅克比迭代法:',jacobi(H, b, 0.1, 200))
print('SOR迭代法:',sor(H, b, 0.001, 20000, 1.5))
# print('Guass-Seidel迭代法:',gauss_seidel(H,b,0.00001,20))
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。