求解非奇异线性代数方程组 Ax=b
1.Jacobi迭代法
令 A = D - L - U
其中D中的对角线元素与A相等,其余为零。L中元素的行数大于列数时和A对应元素相同,其余为零。U中元素的行数小于列数时和A对应元素相同,其余为零。则方程组可以写成 x = B*x + g,其中B = D(逆)*(L+U),g = D(逆)*b。若在右边给定初始向量,则可以计算一个新的向量,以此类推,即为Jacobi迭代法。以下为其代码:
import numpy as np print("解非奇异线性方程组‘Ax=b’") print("请输入非奇异矩阵的行标(n):") n = int(input()) print("请输入非奇异矩阵(A)(一行一行的输入):") A = np.array([([float(0) for i in range(n)]) for j in range(n)]) for i in range(n): for j in range(n): A[i][j] = float(input()) print("请输入向量b:") b = np.array([([0]) for i in range(n)]) for i in range(n): b[i][0] = float(input()) D = np.array([([float(0) for i in range(n)]) for j in range(n)]) for i in range(n): D[i][i] = A[i][i] L = np.array([([float(0) for i in range(n)]) for j in range(n)]) for i in range(n): for j in range(n): if i > j: L[i][j] = -A[i][j] U = np.array([([float(0) for i in range(n)]) for j in range(n)]) for i in range(n): for j in range(n): if i < j: U[i][j] = -A[i][j] Dn = np.linalg.inv(D) # Dn是D的逆 M = Dn@(L+U) t, T = np.linalg.eig(M) # t是特征值,T是特征向量 for i in range(n): if abs(t[i]) > 1: print("迭代法不收敛") g = Dn@b print("请输入向量x(任取解作为初始值):") x = np.array([([float(0)]) for i in range(n)]) k = int(0) while k < 1: x = M@x + g r = A@x - b num = float(0) for i in range(n): num = num + abs(r[i]) if num < 0.0001: # 判断迭代完成的条件 k = 1 print(x)
2.Gauss-Seidel迭代法
注意到Jacobi迭代法中各个分量的计算顺序是没有关系。现在,在计算x(k)用x(k-1)的各个分量计算,但当计算x(k)的第二分量时,由于x(k)的第一个分量已经算出,则用它来代替x(k-1)的第一个分量,其余分量不变。类似的,来计算其它分量,则这种迭代格式为Gauss-Seidel迭代法,简称G-S迭代法。与Jacobi迭代法相比,G-S迭代法编写程序时储存量会更少。同Jacobi迭代法的计量相同,G-S迭代法可以写成 x(k) = (D-L)(逆)*U*x(k-1) + (D-L)(逆)*b。以下为其代码:
import numpy as np
print("解非奇异线性方程组‘Ax=b’")
print("请输入非奇异矩阵的行标(n):")
n = int(input())
print("请输入非奇异矩阵(A)(一行一行的输入):")
A = np.array([([float(0) for i in range(n)]) for j in range(n)])
for i in range(n):
for j in range(n):
A[i][j] = float(input())
print("请输入向量b:")
b = np.array([([0]) for i in range(n)])
for i in range(n):
b[i][0] = float(input())
D = np.array([([float(0) for i in range(n)]) for j in range(n)])
for i in range(n):
D[i][i] = A[i][i]
L = np.array([([float(0) for i in range(n)]) for j in range(n)])
for i in range(n):
for j in range(n):
if i > j:
L[i][j] = -A[i][j]
U = np.array([([float(0) for i in range(n)]) for j in range(n)])
for i in range(n):
for j in range(n):
if i < j:
U[i][j] = -A[i][j]
Dn = np.linalg.inv(D - L) # Dn是D-L的逆
M = Dn@U
g = Dn@b
x = np.array([([float(0)]) for i in range(n)])
k = int(0)
while k < 1:
x = M@x + g
r = A@x - b
num = float(0)
for i in range(n):
num = num + abs(r[i])
if num < 0.0001: # 判断迭代完成的条件(在一定的误差内)
k = 1
print(x)