• • 编辑**：我正在寻找一个算法来解决一个大系统，可解，线性IVP，可以处理非常小的浮点值。解决特征向量和特征值是不可能的numpy.linalg.eig()作为返回值是复杂的，不应该，它也不支持numpy.float128，并且矩阵不是对称的，所以numpy.linalg.eigh()不起作用。Symy可以在无限长的时间内完成，但是在运行了5个小时之后我放弃了。scipy.integrate.solve_ivp()使用隐式方法（试过Radau和BDF），但是输出错误很大。有没有库、方法、算法或解决方案可以处理这么多非常小的数字？

你可以忽略剩下的。
我有一个150x150的稀疏（22500的~500个非零项）矩阵表示一阶系统，线性微分方程。我试图找到这个矩阵的特征值和特征向量，来构造一个函数，作为这个系统的解析解，这样我就可以给它一个时间，它就会给我每个变量的值。我过去曾对类似的40x40矩阵使用过这种方法，它非常（十，在某些情况下是数百次）比scipy.integrate.solve_ivp()更快，并且还使模型后分析更容易，因为我可以使用scipy.optimize.fmin()找到最大值和最大变化率，或者在inf上计算函数以查看如果放置时间足够长，则会沉淀。
然而，这一次，numpy.linalg.eig()似乎不喜欢我的矩阵，给了我复杂的值，我知道这些值是错误的，因为我正在建模的物理系统不可能有复杂的增长率或衰减率（或正弦解），它的变量的复杂值要小得多。我认为这是一个刚性或浮点舍入问题，其中底层LAPACK算法无法处理非常小的值（最小值约为3e-14，并且大多数非零值具有相似的标度）或某些值之间的差异（最大值约为4000，但大于1的值仅显示少数几次）。
我看到过类似用户的问题建议使用sympy来求解特征值，但是当它在5个小时后还没有解出我的矩阵时，我认为这对我的大型系统来说不是一个可行的解决方案。我也看到过使用numpy.real_if_close()来去除复数值的虚部的建议，但是我也不确定这是一个好的解决方案;numpy.linalg.eig()的几个特征值是0，这对我来说是错误的标志，但是另外几乎所有的实部都与虚部具有相同的尺度（非常小），这也使我质疑它们的有效性。我的矩阵是实的，但是不幸的是不对称，所以numpy.linalg.eigh()也不可行。
现在我可以运行scipy.integrate.solve_ivp()任意长的时间（几千个小时），这可能需要很长的时间来计算，然后使用scipy.optimize.curve_fit()来近似我想要的解析解，因为我对它们的形式有很好的了解。这并不理想，因为它使我的程序慢得多，我甚至不确定它是否能解决我在numpy.linalg.eig()中遇到的刚度和圆整问题;我怀疑Radau或BDF将能够导航刚度，但不是四舍五入。
有人有什么想法吗？有没有其他找到特征值的算法可以处理这个问题？numpy.linalg.eig()能和numpy.float128一起工作而不是numpy.float64，或者即使是额外的精度也无济于事？
我很乐意根据要求提供更多的细节。如果需要，我愿意改变语言。