我有这个代码:
import numpy as np
from numpy import exp
def foo(a,b):
return np.array([[exp(-2*b*a), 0,-exp(b*a)/3 + exp(-2*b*a)/3],
[0, exp(-2*b*a), 0],
[-exp(b*a)/3 + exp(-2*b*a)/3,0, 10*exp(4*b*a)/9 - 2*exp(b*a)/9 + exp(-2*b*a)/9]])
a_s = np.linspace(0, 10, 100)
b_s = np.linspace(0.1, 5, 100)
store_det = np.zeros((len(a_s), len(b_s)))
for i, b in enumerate(b_s):
for j, a in enumerate(a_s):
store_det[i,j] = np.linalg.det(foo(a,b))
print(f'max: {np.max(store_det)} \t min: {np.min(store_det)}')
store_det_by_eigen = np.zeros((len(a_s), len(b_s)))
for i, b in enumerate(b_s):
for j, a in enumerate(a_s):
w,v = np.linalg.eig(foo(a,b))
store_det_by_eigen[i,j] = w[0]*w[1]*w[2]
print(f'max: {np.max(store_det_by_eigen)} \t min: {np.min(store_det_by_eigen)}')
对于a
和b
的每一个真实的值,这个方程的行列式应该是1。
考虑到这会产生一个3x3数组,我想计算特征值并在计算中使用它们。特征值的乘积应该等于行列式。然而,例如,对于a=10
和b=5
,乘积不再是1。
我做错什么了吗?我该怎么办?
3条答案
按热度按时间6yoyoihd1#
在@ev-br之上,我想强调的是,如果特征值算法是基于幂运算的,那么它们在数值上可能会更加不稳定(对于不同的特征值幅度)。老实说,我不知道numpy用的是哪一个,但当你的行列式变成一个小特征值和一个大特征值的乘积时,它似乎就崩溃了。
或者,您可以使用高精度库,如mpmath。例如:
sauutmhj2#
当a,B = 10,5时,矩阵的元素数量级为exp(50)和exp(-50)。毫不奇怪,计算会失去精度,因为你根本无法区分1+exp(-50)和1(即使exp(-50)不同于零)。
FWIW,尝试扩展精度很少是一个答案。你可能应该重新考虑你的方法,从你需要这些大/小条目的原因和你实际计算的内容开始。
1tu0hz3e3#
因此,对于一种情况,特征值为:
乘积接近1:
对于引用的问题案例:
不是那么近,但仍然在附近。这并不奇怪,因为eig值要么很小,要么很大。
看看这两个数组:
中间的一行(或一列)几乎全是0。