单步线性定常迭代法

求解非奇异线性代数方程组                      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)