Geoffroy - 1 year ago 112

Python Question

My goal is to create a Python3 program to operate differential operations on a vector V of size N. I did so, test it for basic operation and it works (differentiation, gradient...).

I tried to write with that basis more complex equations (Navier-Stokes, Orr-Sommerfeld,...) and I tried to validate my work by calculating the eigenvalues of these equations.

As these eigenvalues were completely unexpected, I simplify my problem and I am currently trying to calculate the eigenvalues only for the differentiation matrix (see below). But the results seem wrong...

Thanks in advance for your help, because I do not find any solution to my problem...

I use Chebyshev spectral method to operate the differentiation of vectors.

I use the following Chebyshev package (translated from Matlab to Python):

http://dip.sun.ac.za/%7Eweideman/research/differ.html

That package allow me to create a differentiation matrix DM, obtained with:

`nodes, DM = chebyshev.chebdiff(N, maximal_order)`

To obtain the 1st, 2nd, 3th... order differentiation, I write for example:

`dVdx1 = np.dot(DM[0,:,:], V)`

d2Vdx2 = np.dot(DM[1,:,:], V)

d3Vdx3 = np.dot(DM[2,:,:], V)

I tested that and it works.

I've build different operators based on that differentiation process.

I've tried to validate them by finding their eigenvalues. It didn't go well so I am just trying right now with DM only.

I do not manage to find the right eigenvalues of DM.

I've tried with different functions:

`numpy.linalg.eigvals(DM)`

scipy.linalg.eig(DM)

scipy.sparse.linalg.eigs(DM)

sympy.solve( (DM-x*np.eye).det(), x) [for snall size only]

I do not want to directly use the matrix DM, so I wrapped into a function which operates the differentiation (see code below) like that:

`dVdx1 = derivative(V)`

The reason why I do that comes from the global project itself.

This is useful for more complicated equations.

Creating such a function prevents me from using directly the matrix DM to find its eigenvalues (because DM stay inside the function).

For that reason, I use a scipy.sparse.LinearOperator to wrap my method derivative() and use it as an input of scipy.sparse.eig().

Here is the code to compute these eigenvalues:

`import numpy as np`

import scipy

import sympy

from scipy.sparse.linalg import aslinearoperator

from scipy.sparse.linalg import eigs

from scipy.sparse.linalg import LinearOperator

import chebyshev

N = 20 # should be 4, 20, 50, 100, 300

max_order = 4

option = 1

#option 1: building the differentiation matrix DM for a given order

if option == 1:

if 0:

# usage of package chebyshev, but I add a file with the matrix inside

nodes, DM = chebyshev.chebdiff(N, max_order)

order = 1

DM = DM[order-1,:,:]

#outfile = TemporaryFile()

np.save('DM20', DM)

if 1:

# loading the matrix from the file

# uncomment depending on N

#DM = np.load('DM4.npy')

DM = np.load('DM20.npy')

#DM = np.load('DM50.npy')

#DM = np.load('DM100.npy')

#DM = np.load('DM300.npy')

#option 2: building a random matrix

elif option == 2:

j = np.complex(0,1)

np.random.seed(0)

Real = np.random.random((N, N)) - 0.5

Im = np.random.random((N,N)) - 0.5

# If I want DM symmetric:

#Real = np.dot(Real, Real.T)

#Im = np.dot(Im, Im.T)

DM = Real + j*Im

# If I want DM singular:

#DM[0,:] = DM[1,:]

# Test DM symmetric

print('Is DM symmetric ? \n', (DM.transpose() == DM).all() )

# Test DM Hermitian

print('Is DM hermitian ? \n', (DM.transpose().real == DM.real).all() and

(DM.transpose().imag == -DM.imag).all() )

# building a linear operator which wrap matrix DM

def derivative(v):

return np.dot(DM, v)

linop_DM = LinearOperator( (N, N), matvec = derivative)

# building a linear operator directly from a matrix DM with asLinearOperator

aslinop_DM = aslinearoperator(DM)

# comparison of LinearOperator and direct Dot Product

V = np.random.random((N))

diff_lo = linop_DM.matvec(V)

diff_mat = np.dot(DM, V)

# diff_lo and diff_mat are equals

# FINDING EIGENVALUES

#number of eigenvalues to find

k = 1

if 1:

# SCIPY SPARSE LINALG LINEAR OPERATOR

vals_sparse, vecs = scipy.sparse.linalg.eigs(linop_DM, k, which='SR',

maxiter = 10000,

tol = 1E-3)

vals_sparse = np.sort(vals_sparse)

print('\nEigenvalues (scipy.sparse.linalg Linear Operator) : \n', vals_sparse)

if 1:

# SCIPY SPARSE ARRAY

vals_sparse2, vecs2 = scipy.sparse.linalg.eigs(DM, k, which='SR',

maxiter = 10000,

tol = 1E-3)

vals_sparse2 = np.sort(vals_sparse2)

print('\nEigenvalues (scipy.sparse.linalg with matrix DM) : \n', vals_sparse2)

if 1:

# SCIPY SPARSE AS LINEAR OPERATOR

vals_sparse3, vecs3 = scipy.sparse.linalg.eigs(aslinop_DM, k, which='SR',

maxiter = 10000,

tol = 1E-3)

vals_sparse3 = np.sort(vals_sparse3)

print('\nEigenvalues (scipy.sparse.linalg AS linear Operator) : \n', vals_sparse3)

if 0:

# NUMPY LINALG / SAME RESULT AS SCIPY LINALG

vals_np = np.linalg.eigvals(DM)

vals_np = np.sort(vals_np)

print('\nEigenvalues (numpy.linalg) : \n', vals_np)

if 1:

# SCIPY LINALG

vals_sp = scipy.linalg.eig(DM)

vals_sp = np.sort(vals_sp[0])

print('\nEigenvalues (scipy.linalg.eig) : \n', vals_sp)

if 0:

x = sympy.Symbol('x')

D = sympy.Matrix(DM)

print('\ndet D (sympy):', D.det() )

E = D - x*np.eye(DM.shape[0])

eig_sympy = sympy.solve(E.det(), x)

print('\nEigenvalues (sympy) : \n', eig_sympy)

Here are my results (for N=20):

`Is DM symmetric ?`

False

Is DM hermitian ?

False

Eigenvalues (scipy.sparse.linalg Linear Operator) :

[-2.5838015+0.j]

Eigenvalues (scipy.sparse.linalg with matrix DM) :

[-2.58059801+0.j]

Eigenvalues (scipy.sparse.linalg AS linear Operator) :

[-2.36137671+0.j]

Eigenvalues (scipy.linalg.eig) :

[-2.92933791+0.j -2.72062839-1.01741142j -2.72062839+1.01741142j

-2.15314244-1.84770128j -2.15314244+1.84770128j -1.36473659-2.38021351j

-1.36473659+2.38021351j -0.49536645-2.59716913j -0.49536645+2.59716913j

0.38136094-2.53335888j 0.38136094+2.53335888j 0.55256471-1.68108134j

0.55256471+1.68108134j 1.26425751-2.25101241j 1.26425751+2.25101241j

2.03390489-1.74122287j 2.03390489+1.74122287j 2.57770573-0.95982011j

2.57770573+0.95982011j 2.77749810+0.j ]

The values returned by scipy.sparse should be included in the ones found by scipy/numpy, which is not the case. (idem for sympy)

I've tried with different random matrices instead of DM (see option 2) (symmetric, non-symmetric, real, imaginary, etc...), which had small size N (4,5,6..) and also bigger ones (100,...).

That worked

By changing parameters like 'which' (LM, SM, LR...), 'tol' (10E-3, 10E-6..), 'maxiter', 'sigma' (0) in scipy.sparse... scipy.sparse.linalg.eigs always worked for random matrices but never for my matrix DM. In best cases, found eigenvalues are close to the ones found by scipy, but never match.

I really do not know what is so particular in my matrix.

I also dont know why using scipy.sparse.linagl.eig with a matrix, a LinearOperator or a AsLinearOperator gives different results.

I DO NOT KNOW HOW I COULD INCLUDE MY FILES CONTAINING MATRICES DM...

For N = 4 :

`[[ 3.16666667 -4. 1.33333333 -0.5 ]`

[ 1. -0.33333333 -1. 0.33333333]

[-0.33333333 1. 0.33333333 -1. ]

[ 0.5 -1.33333333 4. -3.16666667]]

Every idea is welcome.

May a moderator could tag my question with :

scipy.sparse.linalg.eigs / weideman / eigenvalues / scipy.eig /scipy.sparse.lingalg.linearOperator

Geoffroy.

Answer Source

I spoke with a few colleague and solve partly my problem. My conclusion is that my matrix is simply very ill conditioned...

In my project, I can simplify my matrix by imposing boundary condition as follow:

```
DM[0,:] = 0
DM[:,0] = 0
DM[N-1,:] = 0
DM[:,N-1] = 0
```

which produces a matrix similar to that for N=4:

```
[[ 0 0 0 0]
[ 0 -0.33333333 -1. 0]
[ 0 1. 0.33333333 0]
[ 0 0 0 0]]
```

By using such condition, I obtain eigenvalues for scipy.sparse.linalg.eigs which are equal to the one in scipy.linalg.eig. I also tried using Matlab, and it return the same values.

To continue my work, I actually needed to use the generalized eigenvalue problem in the standard form

λ B x= DM x

It seems that it does not work in my case because of my matrix B (which represents a Laplacian operator matrix). If you have a similar problem, I advise you to visit that question: http://scicomp.stackexchange.com/questions/10940/solving-a-generalised-eigenvalue-problem

(I think that) the matrix B needs to be positive definite to use scipy.sparse. A solution would be to change B, to use scipy.linalg.eig or to use Matlab. I will confirm that later.

I wrote a solution to the stack exchange question I post above which explains how I solve my problem. I appears that scipy.sparse.linalg.eigs has indeed a bug if matrix B is not positive definite, and will return bad eigenvalues.