Контекст:
Моя цель - создать программу Python3 для выполнения дифференциальных операций над вектором V размера N. Я сделал это, протестировал ее для базовой операции, и она работает (дифференциация, градиент...).
Я попытался написать на этой основе более сложные уравнения (Навье-Стокса, Орра-Зоммерфельда,...) и попытался проверить свою работу, вычислив собственные значения этих уравнений.
Поскольку эти собственные значения были совершенно неожиданными, я упрощаю свою задачу и в настоящее время пытаюсь вычислить собственные значения только для матрицы дифференцирования (см. ниже). Но результаты кажутся неправильными...
Заранее благодарю за помощь, так как не могу найти решение своей проблемы...
Определение ДМ:
Я использую спектральный метод Чебышева для дифференцирования векторов. Я использую следующий пакет Чебышева (переведенный с Matlab на Python): http://dip.sun.ac.za/%7Eweideman/research/differ.html
Этот пакет позволяет мне создать матрицу дифференциации DM, полученную с помощью:
nodes, DM = chebyshev.chebdiff(N, maximal_order)
Чтобы получить дифференциацию 1-го, 2-го, 3-го... порядка, я пишу, например:
dVdx1 = np.dot(DM[0,:,:], V)
d2Vdx2 = np.dot(DM[1,:,:], V)
d3Vdx3 = np.dot(DM[2,:,:], V)
Я проверил это, и это работает. Я создал различных операторов на основе этого процесса дифференциации. Я попытался проверить их, найдя их собственные значения. Это не пошло хорошо, поэтому я просто пытаюсь прямо сейчас только с DM. Мне не удается найти правильные собственные значения DM.
Я пробовал с разными функциями:
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]
Почему я использую scipy.sparse.LinearOperator:
Я не хочу напрямую использовать матрицу DM, поэтому я обернул функцию, которая выполняет дифференцирование (см. код ниже):
dVdx1 = derivative(V)
Причина, по которой я это делаю, связана с самим глобальным проектом. Это полезно для более сложных уравнений.
Создание такой функции не позволяет мне напрямую использовать матрицу DM для нахождения ее собственных значений (поскольку DM остается внутри функции). По этой причине я использую scipy.sparse.LinearOperator, чтобы обернуть мой метод производной() и использовать его в качестве ввода scipy.sparse.eig().
Код и результаты:
Вот код для вычисления этих собственных значений:
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)
Вот мои результаты (для 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 ]
Значения, возвращаемые scipy.sparse, должны быть включены в найденные scipy/numpy, что не так. (то же самое для Симпи)
Я пробовал использовать разные случайные матрицы вместо DM (см. вариант 2) (симметричные, несимметричные, действительные, мнимые и т. д.), которые имели как малый размер N (4,5,6..), так и больший размер единицы (100,...). Это сработало
Изменяя такие параметры, как «который» (LM, SM, LR...), «tol» (10E-3, 10E-6..), «maxiter», «sigma» (0) в scipy.sparse... scipy.sparse.linalg.eigs всегда работал для случайных матриц, но никогда для моей матрицы DM. В лучшем случае найденные собственные значения близки к найденным scipy, но никогда не совпадают.
Я действительно не знаю, что такого особенного в моей матрице. Я также не знаю, почему использование scipy.sparse.linagl.eig с матрицей, LinearOperator или AsLinearOperator дает разные результаты.
Я НЕ ЗНАЮ, КАК Я МОГУ ВКЛЮЧИТЬ СВОИ ФАЙЛЫ, СОДЕРЖАЩИЕ МАТРИЦЫ DM...
Для 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]]
Каждая идея приветствуется.
Может ли модератор пометить мой вопрос следующим образом: scipy.sparse.linalg.eigs/weideman/eigenvalues/scipy.eig/scipy.sparse.lingalg.linearOperator
Жоффруа.
