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。