数值计算（迭代法解方程组）
# -*- coding: utf-8 -*-
import numpy as np # a package to calculate
from import identity # to get a identity matrix
import time # calculate time of diffient method
(2) # set seed to make x0 unchanged
def Create_Tridiagonal_Matrices(a, b, c, n): # 创建一种特殊的三对角矩阵，主对角线元素b，比主对角线低一行的对角线上a，比主对角线高一行的对角线上c
A = ((n, n))
if (b <= c or b <= a or b < a + c or n < 2): # 判断是否符合三对角矩阵的基本定义,以及n的个数，1行的三对角毫无意义
print "this is not a Create_Tridiagonal_Matrices,please check a,b,c "
A[0][0] = b;
A[0][1] = c
A[n - 1][n - 1] = b;
A[n - 1][n - 2] = a
for j in range(1, n - 1):
A[j][j - 1] = a
A[j][j] = b
A[j][j + 1] = c
return A
def Jacobi(A, b):
n = [0]#初始化DLU为0
D = (, dtype="float")
L = (, dtype="float")
U = (, dtype="float")
for i in range(n):#根据规则A=D-L-U计算
for j in range(i + 1, n):
U[i][j] = -A[i][j]
L[j][i] = -A[j][i]
for i in range(n):
D[i][i] = A[i][i]
X_k1 = ((n, 1))#初始化X_k1
#X_k2 = (n, 1) # 初始化X_k2，无所谓
D_inv = (D)
B = (D_inv, L + U)#计算相应的B和f
f = (D_inv, b)
#X_k2 = (B, X_k1) + f
while (((((A,X_k1)-b, 2)))/(((b, 2))) > 1e-6): #迭代法不断的趋近
#X_k1 = X_k2
X_k1 = (B, X_k1) + f
return X_k1
def Gauss_Seidel(A, b): #思路类似，在迭代的过程中立即用到了上一个解
# D=((A),0)
n = [0]
D = (, dtype="float")
L = (, dtype="float")
U = (, dtype="float")
for i in range(n):#根据规则A=D-L-U计算
for j in range(i + 1, n):
U[i][j] = -A[i][j]
L[j][i] = -A[j][i]
for i in range(n):
D[i][i] = A[i][i]
# 初始化X_k1，k2
X_k1 = ((n, 1))
#X_k2 = (n, 1)
D_Minus_L_inv = (D - L)
B_G = (D_Minus_L_inv, U) #算出相应B
f_G = (D_Minus_L_inv, b)#算出相应f
#X_k2 = (B_G, X_k1) + f_G
while (((((A, X_k1) - b, 2))) / (((b, 2))) > 1e-6):
#X_k1 = X_k2
X_k1 = (B_G, X_k1) + f_G
return X_k1
def S_O_R(A, b, w): #思路类似，增加了w来加速迭代过程
# D=((A),0)
n = [0]
D = (, dtype="float")
L = (, dtype="float")
U = (, dtype="float")
for i in range(n):
for j in range(i + 1, n):
U[i][j] = -A[i][j]
L[j][i] = -A[j][i]
for i in range(n):
D[i][i] = A[i][i]
X_k1 = ((n, 1))
D_Minus_wL_inv = (D - w * L)
B_w = (D_Minus_wL_inv, (1 - w) * D + w * U)#算出相应B
f_w = w * (D_Minus_wL_inv, b)#算出相应f
while (((((A, X_k1) - b, 2))) / (((b, 2))) > 1e-6):
X_k1 = (B_w, X_k1) + f_w
return X_k1
def timecmp(A,b,fun1,fun2,fun3):#通过运行100次的平均时间较各个函数的效率
time1=((10,1))
time2 = ((10, 1))
time3 = ((10, 1))
for i in range(10):
t1=()
Jacobi(A, b)
time1[i][0]=()-t1
t2 = ()
Gauss_Seidel(A, b)
time2[i][0] = () - t2
t3 = ()
S_O_R(A, b, 1.5)
time3[i][0] = () - t3
return (time1),(time2),(time3)
n=100
A = Create_Tridiagonal_Matrices(-1, 2, -1, n)
print "\n 生成的系数行列式是", A
b = ((n, 1))
print "\n 生成的b是", b
print "\n 调用函数算出来解为", (A, b)
t1 = ()
X = Jacobi(A, b)
print "计算花费时间",()-t1
print 'Jacobi计算出来的解:',X
t1 = ()
X = Gauss_Seidel(A, b)
print "\n计算花费时间",()-t1
print 'Gauss_Seidel计算出来的解:',X
print X
t1 = ()
X = S_O_R(A, b, 0.5)
print "\n计算花费时间",()-t1
print 'S_O_R计算出来的解:',X
print "over"
print "Jacobi,Gauss_Seidel,S_O_R 在10次运行的平均时间",timecmp(A,b,Jacobi,Gauss_Seidel,S_O_R)