Sympy：解决 Hessian 矩阵中的条目以获得更好的可读性？

Question

我对 sympy 很陌生，我不知道如何以格式良好的方式生成输出。现在我已经计算了潜在函数的 Hessian 矩阵：

V = 1/2*kOH*(r1)**2 +1/2*kOH*(r2)**2 +1/2*kHH*(r3)**2

三个谐振子项的通称形式为：

1/2*k*r**2.

所有变量都是正实数。

我的问题是，当我打印我的矩阵时，条目尚未解决，仅以功能方式显示。我希望条目在已经执行偏导数之后的形式中，而不仅仅是显示矩阵中每个点需要执行的推导。

def Hessian():

    '''
    sympy calc of hessian Matrix H for IR normal modes analysis
    from a potential V.

    Must be multiplicable with 9x9 matrix (somehow) in 
    the equation: F = M**(-1/2) * H * M**(-1/2)
    Here, F is the mass weighted Hessian, whose Eigenvalues
    contain the frequencies of the normal modes of water.
    M comes from the multiplication of the 3N-Dimensional
    mass-vector m with a 3N-dimensional identity matrix:
    M = m*I, I.shape = 3*N, 3*N, N = number of atoms in water. 
    '''

    kOH, kHH, r1, r2, r3 = sy.symbols('kOH kHH r1 r2 r3', real=True, positive=True)

    V = sy.Function('V')(1/2*kOH*(r1)**2 +1/2*kOH*(r2)**2 +1/2*kHH*(r3)**2)

    f = sy.hessian(V,[r1, r2, r3])

    sy.pprint(f)


Hessian()

补充：这实际上不是计算方面的一部分，因此也不是问题的一部分，但是如果有人知道他们在科学方面的知识：你能告诉我如何(3,3) 取决于三个距离的潜在依赖的 Hessian 应该与 (9,9) 质量矩阵相乘？函数的注释包含科学背景，如果你感兴趣的话。

Answer 1

你遇到的基本问题是：

In [39]: f = Function('f')                                                                                                        

In [40]: f(x)                                                                                                                     
Out[40]: f(x)

In [41]: f(x).diff(x)                                                                                                             
Out[41]: 
d       
──(f(x))
dx      

In [42]: f(x).diff(x).subs(x, 2*y)                                                                                                
Out[42]: 
⎛d       ⎞│     
⎜──(f(x))⎟│     
⎝dx      ⎠│x=2⋅y

理想情况下，SymPy 会将最后一个结果表示为 f'(2y) 之类的东西，但 SymPy 没有办法直接表示此类对象。理想情况下会有一个微分运算符 D，这样说 D(f)(x) 就等于 f(x).diff(x)。这样你就可以将其表示为 D(f)(2*y) 当然可以显示为 f'(2y).

当然，如果您在这里用一个函数代替 f，那么导数可以求值：

In [45]: f(x).diff(x).subs(x, 2*y).subs(f, Lambda(t, t**3))                                                                       
Out[45]: 
⎛d ⎛ 3⎞⎞│     
⎜──⎝x ⎠⎟│     
⎝dx    ⎠│x=2⋅y

In [46]: _.doit()                                                                                                                 
Out[46]: 
    2
12⋅y

要回答您的其他问题，显然您不能将 9x9 矩阵与 3x3 矩阵相乘。 F 的等式意味着 H 和 M 都是正方形且大小相同。要么你的质量矩阵真的只是 3x3，要么你的潜在函数实际上是 9 个坐标的函数。假设 r1 是原子 1 和原子 2 之间的距离，那么可能 r1 = sqrt((x1 - x2)**2 + (y1 - y2)**2 + (z1 - z2)**2) 在这种情况下你应该计算你的 Hessian wrt x1 等而不是 r1.

Sympy：解决 Hessian 矩阵中的条目以获得更好的可读性？

Sympy: Solve entries in Hessian Matrix for better readability?

python

sympy

hessian-matrix