Я пытался получить собственные значения и собственные векторы произведения экспонент матрицы с помощью функции linalg.eig()
из NumPy и столкнулся с большими числовыми ошибками.
Произошло следующее произведение:
def matrix(t): # matrix with n rows and n columns
...
S = scipy.linalg.expm( -1j * matrix(0) )
for i in range(N):
S_prev = S
S_next = scipy.linalg.expm( -1j * matrix(0 + i * delta_t) )
S = numpy.matmul(S_next, S_prev)
Я должен заметить, что я могу не просто написать S = scipy.linalg.expm( -1j * (matrix(0) + matrix(delta_t) + ... )
, потому что матрицы в показателе степени с различными t
не коммутируют.
Все собственные значения S
должны иметь абсолютное значение 1, поскольку в показателе степени есть мнимая единица 1j
, а matrix(t)
- эрмитово (все ее собственные значения действительны). Я проверяю этот факт:
vals, funcs = numpy.linalg.eig(S)
print(numpy.absolute(vals))
и вижу, что первое и последнее значения далеки от 1:
[2.01401514e+07 1.75005607e+03 1.18825323e+03 7.62941962e+02
4.93897530e+02 3.56924391e+02 3.10300966e+02 2.54560554e+02
1.84934167e+02 1.34093245e+02 1.25981874e+02 9.58105563e+01
5.12598735e+01 3.87933877e+01 2.88031162e+01 1.59672649e+01
1.18759460e+01 9.29488363e+00 8.71429012e+00 8.43350092e+00
8.43464605e+00 4.51898108e+00 2.62851258e+00 1.71905799e+00
1.25618921e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00
1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00
9.88249716e-01 1.00000000e+00 1.00000000e+00 1.01188999e+00
...
1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00
9.99999997e-01 1.00000000e+00 2.21288824e-01 1.00254283e+00
1.00036087e+00 9.99990372e-01 9.99542072e-01 1.18574662e-01
1.18558625e-01 1.14754055e-01 1.07586057e-01 8.42037792e-02
6.26281297e-02 3.47184828e-02 2.57776048e-02 1.95084377e-02
7.93765396e-03 2.37035650e-08 1.04372635e-02 7.45749680e-03
5.40732797e-03 3.92833442e-03 2.80172313e-03 3.22267950e-03
2.02470192e-03 1.31073074e-03 8.41571767e-04 5.71394099e-04]
Также я смотрю на собственные векторы хороших собственных значений (абсолютное значение которых равно 1) и вижу совершенно неправильное поведение, такое как высокочастотные колебания:
num = numpy.arange(n)
matplotlib.pyplot.plot(num, np.absolute(funcs[200])) # good eigenvalue number is 200
matplotlib.pyplot.ylabel("Vector Component Value")
matplotlib.pyplot.xlabel("Vector Component Number")
matplotlib.pyplot.show()
График
(Я беру абсолютное значение компонентов собственного вектора, потому что оно пропорционально плотности вероятности в моей задаче)
Можно ли исправить ошибки?