python – Scipy错误地解决了矩阵的特征问题

2023年4月15日 155次阅读

当我解决奇异矩阵的特征问题时,我正面临SciPy非常奇怪的行为,即如果我通过某些函数(下面的代码中的矩阵)生成矩阵,则计算的特征值是不正确的.但是,如果我在(der2)中手动键入矩阵,则对角化似乎会产生适当的结果.这也可以通过减去两个矩阵来检查,这些都是在下面的代码中完成的.

代码是

import numpy as np
import scipy as sp
from scipy.linalg import eigvals

def cbar(k, n):
    """
    cbar function for coefficients
    """
    if k==0 or k==n-1:
        return np.float128(2.)
    else:
        return np.float128(1.)    

def ChebCollDer(x):
    """
    ChebCollDer  Chebyshev collocation differentiation matrix.
    """
    xx=np.array(x)  
    n=xx.size
    d=np.zeros((n,n))
    for i in range(n):
        for k in range(n):
            if i!=k:
                d[i,k]=cbar(i, n)*np.float128(sp.power(-1., i+k))/(xx[i]-xx[k])/cbar(k,n)
    for i in range(n):
        tmp=-sp.sum(d[i,:])
        d[i,i]=tmp

    return d
nn=5
xx=(np.cos(sp.pi*np.linspace(0,1.0,nn)))/2.
der=ChebCollDer(xx)
print eigvals(der)
der2=[[ 11.0 ,-13.656854,  4.0 ,-2.3431458 , 1.0],
 [ 3.4142136, -1.4142136, -2.8284271,  1.4142136, -0.58578644],
 [-1.0,  2.8284271,  1.110223e-16, -2.8284271,  1.0]
 ,[0.58578644, -1.4142136,  2.8284271,  1.4142136, -3.4142136],
 [-1.0,2.3431458,-4.0 ,13.656854,-11.0]]
print eigvals(der2)

print der-der2

结果是：
矩阵的特征值：

[ 0.00389434+0.00282825j  0.00389434-0.00282825j -0.00148641+0.00457958j -0.00148641-0.00457958j -0.00481586+0.j        ]

der2的特征值：

[  9.71161644e-02+0.j          -9.71161490e-02+0.j          -3.08158279e-08+0.j 7.69629619e-09+0.09711159j   7.69629619e-09-0.09711159j]

看到der2有一个数值为零的特征值,应该是因为矩阵der有一个零特征向量,它只是[1,1,1,1,1]
der-der2的最大元素是10E-08阶.
我怀疑,有一些类型转换问题,但不知道它来自哪里.

最佳答案 这不是一个完整的答案,但我没有足够的声誉来评论.

您可能会在数学页面上发布更多运气,因为我认为这是一个数值分析问题.

如果我这样做,v = sp.linalg.eig(der)那么所有你似乎都满足eignevalue方程,因为der @ u – s * u接近于零.

另一个有趣的事情是,der的特征值都具有相等的模数,并且具有围绕该半径的圆等距间隔的参数.像缩放的第5根-1.

实际上,如果从scipy输出中获取特征向量的同等加权线性组合(即右乘v乘所有向量),您将接近对应于零特征值的常数向量.

那就是说,我不太确定发生了什么.

此外,你的函数ChebCollDer将输出一个常规的浮点数而不是float128,因为你没有在float128初始化它.这意味着即使你在float128中进行了一些计算,它们也会在存储到d之前被转换为默认的scipy浮点数.

点赞