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