计算方法0

高斯列主元消去法

import numpy as np

size = int(input())
a = []
b = []
for i in range(size):
    a.append(list())
    for j in range(size):
        a[i].append(int(input()))
    b.append(int(input()))
for k in range(size):
    c = list(zip(a, b))
    c.sort(key=lambda x: x[0], reverse=True)
    a, b = zip(*c)
    a = list(a)
    b = list(b)
    for i in range(k + 1, size):
        tmp = a[i][k]
        for j in range(size):
            a[i][j] = a[i][j] - tmp/a[k][k]*a[k][j]
        b[i] = b[i] - tmp/a[k][k]*b[k]

a = np.array(a)
b = np.array(b)
for k in range(size - 1, -1, -1):
    tmp = b[k]
    b[k] = 0
    b[k] = (tmp - sum(a[k]*b))/a[k][k]
print(b)

# 示例
# 1 2 3 14
# 2 5 2 18
# 3 1 5 20
# 1 2 3

雅可比迭代法

import numpy as np

size = int(input())
a = []
b = []
for i in range(size):
    a.append(list())
    for j in range(size):
        a[i].append(float(input()))
    b.append(float(input()))
times = 0
a = np.array(a)
b = np.array(b)

D = np.eye(size)*a
U = np.triu(a) - D
L = (np.triu(a.T) - D).T
x = np.zeros(size)
while True:
    times += 1
    tmp = x
    x = (np.eye(size) - np.linalg.inv(D)@a)@x + np.linalg.inv(D)@b
    if max(abs(tmp - x)) < 0.001:
        break
    print("{}\n".format(x))
print(x)
print(times)

# 示例
# 10 -2 -1 3
# -2 10 -1 15
# -1 -2 5 10
# 1 2 3

高斯赛德尔迭代法

import numpy as np

size = int(input())
a = []
b = []
for i in range(size):
    a.append(list())
    for j in range(size):
        a[i].append(float(input()))
    b.append(float(input()))
times = 0
a = np.array(a)
b = np.array(b)

D = np.eye(size)*a
U = np.triu(a) - D
L = (np.triu(a.T) - D).T
x = np.zeros(size)
while True:
    times += 1
    tmp = x
    x = -np.linalg.inv(D + L)@U@x + np.linalg.inv(D + L)@b
    if max(abs(tmp - x)) < 0.001:
        break
    print("{}\n".format(x))
print(x)
print(times)

# 示例
# 10 -2 -1 3
# -2 10 -1 15
# -1 -2 5 10
# 1 2 3

超松驰迭代法

import numpy as np
A = np.array([[10, -2, -1], [-2, 10, -1], [-1, -2, 5]])
b = np.array([[3], [15], [10]])
omega_list = [1 + 0.005 * i for i in range(100)]
length = len(A)
count = []
tmp = 0
x_new = 0
for omega in omega_list:
    times = 0
    x_0 = np.zeros((length, 1))
    while True:
        x_hold = x_0.copy()
        x_new = x_0.copy()
        for i in range(length):
            x_new[i][0] = x_0[i][0] + omega * (b[i][0] - sum([A[i][j] * x_new[j][0] for j in range(i)]) - sum(
                [A[i][j] * x_0[j][0] for j in range(i, length)])) / A[i][i]
            tmp = x_0
            x_0 = x_new.copy()
        times += 1
        if max(abs(tmp - x_new)) < 0.001:
            break
    count.append(times)
print(x_new)
print(count)
print()

最后说一句，python真香。-_-