Симп: Расчет собственных значений занимает очень много времени для (9,9) матрицы - PullRequest
0 голосов
/ 19 сентября 2019

Я пытаюсь вычислить собственные значения матрицы F, которая содержит 9 переменных, которые являются декартовыми координатами трех векторов.Время выполнения F.eigenvals() занимает не менее 15 минут, и я не набрался терпения, чтобы еще дольше ждать завершения программы.Я подозреваю, что я делаю что-то ужасно неправильно, что вызывает очень долгий расчет.Может быть, это смешивание объектов numpy и sympy, которое вызывает проблему.

У меня есть три атома воды с их соответствующими массами, из которых я вычисляю так называемую матрицу масс M, используя функцию Mmat.

import sympy as sy
import numpy as np
from scipy.constants import u as u

def Mmat(*args):

    N = len(args)
    args = np.array([[arg*u]*3 for arg in args]).reshape(N*3)
    m = args**(-1/2)
    I = np.identity(3*N)

    return m*I

mO, mH1, mH2 = 16, 1, 1
M = Mmat(mH1,mH2,mO)

Эту матрицу масс необходимо умножить по формуле F = M*H*M на матрицу Гессе H, чтобы получить (9,9) матрицу F, из которой собственные значенияизвлечены.Сам гессиан определяется частными производными потенциала V для декартовых координат трех атомов O, H1 и H2.

Расчет матрицы F сам занимает всего секунду, ноприменение метода .eigenvals() занимает вечность.Вот функция для вычисления H:

def dist(v1,v2):

    return ((v2-v1)**2).sum()**(-1/2)

def Hessian():

    kOH, kHH, dr1, dr2, dr3 = sy.symbols('kOH kHH dr1 dr2 dr3', real=True, positive=True)

    V = 1/2*kOH*(dr1)**2 +1/2*kOH*(dr2)**2 +1/2*kHH*(dr3)**2

    x1,x2,x3,x4,x5,x6,x7,x8,x9 = sy.symbols('x1 x2 x3 x4 x5 x6 x7 x8 x9', real=True)
    dOH10, dOH20, dH1H10 = sy.symbols('dOH10 dOH20 dH1H10', real=True, positive=True)

    vO = np.array([x1,x2,x3])
    vH1 = np.array([x4,x5,x6])
    vH2 = np.array([x7,x8,x9])

    r1 = dist(vO,vH1) - dOH10
    r2 = dist(vO,vH2) - dOH20
    r3 = dist(vH1,vH2) - dH1H10

    V = V.subs(dr1,r1).subs(dr2,r2).subs(dr3,r3)

    H = sy.hessian(V,[x1,x2,x3,x4,x5,x6,x7,x8,x9])

    return H

F = M*Hessian()*M

Реальная проблема возникает только в конце:

print(F.eigenvals())

Я смешиваю классы объектов неправильно?Это тот факт, что в матрице масс есть большие числа (6.022e26)?Это нормально, если матрица такого размера работает так медленно?

Заранее спасибо.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...