理论推导请参考《数值分析》李庆扬
完整代码放置在GitHub
Gauss消去法解方程组
import numpy as np
A = np.array([[0.4096,0.1234,0.3678,0.2943,0.4043],[0.2246,0.3872,0.4015,0.1129,0.1550],
[0.3645,0.1920,0.3781,0.0643,0.4240],[0.1784,0.4002,0.2786,0.3927,-0.2557]])
print(A)
rows,columns = A.shape
ans = np.zeros(shape=(rows,1))
for c1 in range(rows-1):
for c2 in range(rows-c1-1):
if A[c1+c2+1 , c1] != 0:
a = A[c1+c2+1 , c1] / A[c1 , c1]
A[c1+c2+1:c1+c2+2 , :] = A[c1+c2+1:c1+c2+2 , :] - A[c1:c1+1 , :] * a
print(A)
#--------------------代回求解-------------------------
for i in range(rows):
ans[rows-i-1] = A[rows-i-1 , columns-1] / A[rows-i-1 , columns-2-i]
A[: , columns-1:columns] = A[: , columns-1:columns] - ans[rows-i-1] * A[: , rows-i-1:rows-i]
print(ans)输出：
[[ 0.4096 0.1234 0.3678 0.2943 0.4043]
[ 0.2246 0.3872 0.4015 0.1129 0.155 ]
[ 0.3645 0.192 0.3781 0.0643 0.424 ]
[ 0.1784 0.4002 0.2786 0.3927 -0.2557]]
[[ 4.09600000e-01 1.23400000e-01 3.67800000e-01 2.94300000e-01
4.04300000e-01]
[-2.77555756e-17 3.19534863e-01 1.99820605e-01 -4.84764160e-02
-6.66937988e-02]
[ 7.13898500e-18 0.00000000e+00 -5.98156607e-04 -1.85126813e-01
8.13706787e-02]
[-1.14249713e-15 0.00000000e+00 1.38777878e-17 3.07244855e+01
-1.37247573e+01]]
[[-0.18191778]
[-1.66303081]
[ 2.21722832]
[-0.44670422]]列主元素消去法求解方程组
import numpy as np
A = np.array([[0.4096,0.1234,0.3678,0.2943,0.4043],[0.2246,0.3872,0.4015,0.1129,0.1550],
[0.3645,0.1920,0.3781,0.0643,0.4240],[0.1784,0.4002,0.2786,0.3927,-0.2557]])
np.set_printoptions(precision=9)
print(A)
rows,columns = A.shape
ans = np.zeros(shape=(rows,1))
for c1 in range(rows-1):
k = c1
for j in range(c1+1 , rows):
if A[k,c1] < A[j,c1]:
k = j
if k != c1:
A[[c1 , k]] = A[[k , c1]]
for c2 in range(rows-c1-1):
if A[c1+c2+1 , c1] != 0:
a = A[c1+c2+1 , c1] / A[c1 , c1]
A[c1+c2+1] = A[c1+c2+1] - A[c1] * a
for i in range(rows):
ans[rows-i-1] = A[rows-i-1 , columns-1] / A[rows-i-1 , columns-2-i]
A[: , columns-1:columns] = A[: , columns-1:columns] - ans[rows-i-1] * A[: , rows-i-1:rows-i]
print(ans)输出：
[[ 0.4096 0.1234 0.3678 0.2943 0.4043]
[ 0.2246 0.3872 0.4015 0.1129 0.155 ]
[ 0.3645 0.192 0.3781 0.0643 0.424 ]
[ 0.1784 0.4002 0.2786 0.3927 -0.2557]]
[[-0.181917777]
[-1.663030813]
[ 2.217228318]
[-0.446704219]]追赶法求解方程组
#----------矩阵分解----------
L = np.zeros(shape=(n,n))
U = np.eye(n)
b = A[off0_i , off0_j].reshape((n,1))
c = A[off1_i , off1_j].reshape((n-1,1))
gama = A[offd1_i , offd1_j].reshape((n-1,1))
alpha = np.zeros(shape=(n , 1))
beta = np.zeros(shape=(n-1 , 1))
L[offd1_i , offd1_j] = A[offd1_i , offd1_j]
alpha[0 , 0] = b[0 , 0]
beta[0 , 0] = c[0 , 0] / alpha[0 , 0]
#注意gama的数组大小不同
for i in range(1 , n):
alpha[i , 0] = b[i , 0] - gama[i-1 , 0] * c[i-1 , 0] / alpha[i-1 , 0]
if i != n-1:
beta[i , 0] = c[i , 0] / alpha[i , 0]
L[off0_i , off0_j] = alpha.reshape((n,))
U[off1_i , off1_j] = beta.reshape((n-1,))
#可以验证分解是否正确
print(np.dot(L,U))
#---------回代求解-----------
y = np.zeros(shape=(n,1))
x = np.zeros(shape=(n,1))
y[0,0] = B[0,0] / alpha[0,0]
for i in range(1,n):
y[i,0] = (B[i,0] - gama[i-1,0] * y[i-1,0]) / alpha[i,0]
x[n-1 , 0] = y[n-1 , 0]
for i in range(n-2 , -1 , -1):
x[i,0] = y[i,0] - beta[i,0] * x[i+1,0]
print(x.reshape((10,-1)))Jacobi迭代法求解方程组
import numpy as np
AB = np.array([[8,-3,2,20],[4,11,-1,33],[6,3,12,36]])
np.set_printoptions(precision=12)
rows,columns = AB.shape
A = AB[: , 0:rows]
B = AB[: , rows:columns]
ans1 = np.zeros(shape=(rows,1))
ans2 = np.zeros(shape=(rows,1))
a = np.diagonal(A) #取出对角线的元素
a = a.reshape((3, 1)) #将一维数组转换成二维数组
A = -1*A
i, j = np.diag_indices_from(A) #取出对角线元素的坐标索引值
A[i,j] = 0
A = A/a
B = B/a
err = 1
while err >= 1.e-10 :
ans2 = np.dot(A,ans1) + B
err = np.max(np.abs(ans1 - ans2))
ans1 = ans2.copy()
print(ans2)输出：
[[3.000000000013]
[2.000000000013]
[0.999999999992]]Gauss_seidel迭代法求解方程组
def gauss_seidel_iteration_matrix(A):
n = len(A)
D = np.diag(np.diag(A))
L = np.tril(A, k=-1)
U = np.triu(A, k=1)
D_inv = np.linalg.inv(D)
M = np.dot(D_inv, L + U)
return M
def spectral_radius(matrix):
eigenvalues, _ = np.linalg.eig(matrix)
return max(abs(eigenvalues))
def main(A):
# 计算 Gauss-Seidel 迭代矩阵
iteration_matrix = gauss_seidel_iteration_matrix(A)
# 计算谱半径
radius = spectral_radius(iteration_matrix)
print("Gauss-Seidel 迭代矩阵:")
print(iteration_matrix)
print("谱半径:", radius)
#----------------计算谱半径---------------------
if __name__ == '__main__':
main(A)
#---------------------------------------------
ans1 = np.zeros(shape=(n,1))
ans2 = np.zeros(shape=(n,1))
a = np.diagonal(A).reshape((n, 1)) #取出对角线的元素
A = -1*A
i, j = np.diag_indices_from(A) #取出对角线元素的坐标索引值
A[i,j] = 0
A = A/a
B = B/a
err = 3
k = 0
while err >= 1.e-5:
ans_tem = ans1.copy()
for i in range(n):
ans2[i,0] = np.dot(A[i:i+1] , ans1) + B[i,0]
ans1[i,0] = ans2[i,0]
err = np.max(np.abs(ans_tem - ans2))
k = k + 1
#print(k) #打印已迭代次数
x = ans2.reshape((10,-1))
print(x)