Scipy solves eigenproblem of a matrix incorrectly

user_3602168 0
Viewed 168

I am facing a very weird behaviour of SciPy when solving eigenproblem of a singular matrix, namely the computed eigenvalues are incorrect if I generate the matrix by some functions (matrix der in the code below). However, if I manually type the matrix in (der2) the diagonalisation seems to yield proper results. This can be also checked by subtracting both matrices, which is all done in the code below.

the code is

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

and the results are: eigenvalues of matrix der:

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

eigenvalues of 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]

see that der2 has one eigenvalue that is numerically zero, as should be since the matrix der has a null eigenvector which is simply [1,1,1,1,1] The biggest element of der-der2 is of order 10E-08. I suspect, that there is some type conversion problem, but have no idea where it comes from.

scipy linear-algebra matrix python
1 Answers

This isn't a full answer but I don't have enough reputation to comment.

You might have more luck posting on the math page since I think this is a numerical analysis problem.

If I do s,v = sp.linalg.eig(der) then all the u seem to satisfy the eignevalue equation in the sense that der@u - s*u is close to zero.

Another interesting thing is that the eigenvalues of der all have equal modulus, and are are have arguments equally spaced around the circle of this radius. Like scaled 5-th roots of -1.

In fact, if you take and equally weighted linear combination of the eigenvectors from the scipy output (i.e. right multiply v by the all ones vector) you will get close to a constant vector corresponding to the zero eigenvalue.

That said, I'm not very sure what is going on.

Also, your function ChebCollDer will output a regular ndarray of float instead of float128 since you do not initialize it at float128. That means that even though you do some calculations in float128, they get cast to the default scipy float before they are stored in d.

user_39012160
Related