x1c 0d1x求解线性方程组的方法本身工作正确。用solve测试。但是当试图求解大量矩阵时，某些解向量不与真实向量收敛，因此，无法获得误差的正态分布。上面的直方图显示，有几个值非常大-这些都是解错的线性方程组。这些情况是用同样的函数分别推导和求解的，但得到了正确的解
import numpy as np
import matplotlib.pyplot as plt
def calculate_relative_error(true_solution, computed_solution):
n = len(true_solution)
rmse = np.sqrt(np.sum((true_solution - computed_solution)**2))/n
sup_norm = np.max(np.abs(true_solution - computed_solution))
return rmse, sup_norm
def generate_random_matrix(n = 6 ):
A = (np.random.random((n, n)).astype(np.float64) * 2 - 1).astype(np.float64)
while np.linalg.det(A) == 0:
generate_random_matrix(n)
return A
def tridiagonal_matrix(n):
main_diag = (np.random.random(n).astype(np.float64) * 2 - 1).astype(np.float64)
sub_diag = (np.random.random(n-1).astype(np.float64) * 2 - 1).astype(np.float64)
super_diag = (np.random.random(n-1).astype(np.float64) * 2 - 1).astype(np.float64)
A = np.diag(main_diag) + np.diag(sub_diag, k=-1) + np.diag(super_diag, k=1)
while np.linalg.det(A) == 0:
tridiagonal_matrix(n=6)
return A
def gauss(A, b, pivoting):
n = len(b)
a = np.hstack((A, b[:, np.newaxis])).astype(np.float64)
for i in range(n):
if pivoting:
max_index = np.argmax(np.abs(a[i:, i])) + i
a[[i, max_index]] = a[[max_index, i]]
for j in range(i + 1, n):
factor = a[j, i] / a[i, i]
a[j, i:] -= factor * a[i, i:]
x = np.zeros(n, dtype=np.float64)
for i in range(n - 1, -1, -1):
x[i] = (a[i, -1] - np.dot(a[i, i + 1:-1], x[i + 1:])) / a[i, i]
return x
def thomas(A, b):
gamma = [-A[0][1] / A[0][0]]
beta = [b[0] / A[0][0]]
n = len(b)
x = [0] * n
for i in range(1, n):
if i != n - 1:
gamma.append(-A[i][i + 1] / (A[i][i - 1] * gamma[i - 1] + A[i][i]))
beta.append((b[i] - A[i][i - 1] * beta[i - 1]) / (A[i][i - 1] * gamma[i - 1] + A[i][i]))
x[n - 1] = beta[n - 1]
for i in range(n - 2, -1, -1):
x[i] = gamma[i] * x[i + 1] + beta[i]
return x
num_matrices = 1000
methods = ["gauss_no_pivoting", "thomas"]
for method in methods:
errors_rmse = []
errors_sup_norm = []
for _ in range(num_matrices):
A_random = generate_random_matrix(6)
A = A_random.copy()
A_tridiagonal = tridiagonal_matrix(6)
b = np.array([1, 1, 1, 1, 1, 1]).astype(np.float64)
true_solution = gauss(A_random, b, pivoting=True)
computed_solution = None
if method == "gauss_no_pivoting":
computed_solution = gauss(A_random.copy(), b.copy(), pivoting=False)
elif method == "thomas":
computed_solution = thomas(A_tridiagonal.copy(), b.copy())
rmse, sup_norm = calculate_relative_error(true_solution, computed_solution)
errors_rmse.append(rmse)
errors_sup_norm.append(sup_norm)
plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.hist(errors_rmse, bins=20, color='blue', edgecolor='black')
plt.title(f'{method} - RMSE Histogram')
plt.xlabel('RMSE')
plt.ylabel('Frequency')
plt.subplot(1, 2, 2)
plt.hist(errors_sup_norm, bins=20, color='green', edgecolor='black')
plt.title(f'{method} - Sup Norm Histogram')
plt.xlabel('Sup Norm')
plt.ylabel('Frequency')
plt.tight_layout()
plt.show()字符串
我不明白异常大的sup_norm和rmse值出现在哪里，这会阻止获得正态分布。正如您在直方图中看到的，所有真实值都接近第一列，因此它们在图形上的一个位置累积。我想消除这样的大错误
1条答案
按热度按时间oyjwcjzk1#
你在随机选择中没有正确的避免奇异矩阵。
如果您尝试绘制
np.linalg.det(A)与rmse的关系图，您将看到只有当行列式非常接近0时，RMSE才会出现高值添加
字符串
在
errors_rmse和errors_sup_norm的类似线之后，然后如果你把
型
x1c 0d1x的数据
或以对数刻度
表示
尝试限制随机选择的矩阵，以明确非奇异矩阵，那么你永远不会看到均方根误差
你试图防止奇异矩阵的两个错误是
Error 1：递归错误
以下代码无效：
型
它递归地回调
tridiagonal_matrix，但丢弃结果，所以你在最后返回的只是原始的A，即使它的行列式是0。注意，一个校正可以是
型
但是如果你认为在递归的噩梦中，它是从...
while开始的，这里的递归有点多余，如果你想使用递归，你可以这样做。型
在scheme或caml中，这是一种优雅的方法。但在python中，就不那么优雅了。因为python不是一种终端递归语言。这意味着如果在找到矩阵之前持续太久，它会因为堆栈溢出（epsilon错误）而失败。
所以最好放弃递归而保留while
型
如果您真的不想使用
while True: ... if ...: break（有些人会武断地拒绝任何while True），您可以型
错误2：永远不要用==比较浮点数
不是
==0不是一个充分条件。首先，浮点数不是精确的数字。参见臭名昭著的0.1+0.2!=0.3。其次，即使与精确数字（如0）进行比较，由于数字错误的积累，大多数0也永远不会是精确的0。
为了在线性代数上下文中重用我之前的0.1+0.2-0.3示例，
np.linalg.det(0.1*np.identity(6)+0.2*np.identity(6)-0.3*np.identity(6))不是0。所以你总是需要一个容差，在这种情况下，因为你显然只是想检查方法的数值误差（而不是初始矩阵的数值误差），容差没有理由非常宽松。
所以我会
型
但如果这对您的需要来说“太容易”，您可以减少
atol用这种方法保证矩阵不奇异，得到了高斯方法所期望的结果
的
并证明您的托马斯实现存在问题（您的初始直方图似乎显示了几乎总是有效的东西。但这是因为由于奇异矩阵导致的巨大RMSE的情况非常少。但如果您删除奇异矩阵，因此放大“0”条，您会看到一般情况下RMSE总是太大。“rmse = 0”只是一个随机rmse值，不可能比其他更可能）
的
tl;干
1.这种分布是由于奇异矩阵
1.当matrix为单数时，您的“重试”方式将被忽略，因为您总是返回第一次尝试
1.你测试一个矩阵是否为奇异矩阵的方法（det==0）太严格了。在浮点数的数值计算中，总是有一个公差
1.一旦真正避免了奇异矩阵，你会看到分布集中在非常低的高斯误差值上。即使是非常罕见的极端值，这个分布仍然非常低（最多10阶）。
1.另一方面，一旦托马斯真正避免了奇异矩阵，您会看到该分布确实不包含错误的疯狂值，但RMSE>1是一个非常典型的情况，这只是证明您的托马斯实现不起作用（它不是“在一些奇怪的情况下不起作用”;它永远不会工作，而几个0 -很明显，如果你只放大第一个酒吧，你会看到这个0酒吧是由很多0.2，0.3，0.4，0.8，. -极少数“=0在数值误差范围内”就像一天两次停止的时钟是正确的）。
因为这不是你的问题，我有点懒，我不会试图理解为什么你的托马斯实现不工作，但它不是。我猜你最初的目标是看看它是否工作，所以使命完成：D