如何从Sympy表达式中提取线性二次项并求解M/V/K/A/B/C矩阵

我有一个依赖变量数组X(t)和dX(t)的通用表达式E，希望找到简便方法从中提取线性项和二次项，现有代码仅对简单表达式生效。

我的目标是得到矩阵M、V、K、A、B和常量C，满足如下关系：

E = dX @ M @ dX + dX @ V @ X + X @ K @ X + A @ dX + B @ X + C

示例场景 #

import sympy as sp

t = sp.symbols("t")
x = sp.Function("x")(t)
y = sp.Function("y")(t)
dx = sp.diff(x, t)
dy = sp.diff(y, t)
X = sp.Array([x, y])
dX = sp.Array([dx, dy])

E = 7*(dx**2) + 2*dx*dy + 6*(dy**2)  # M = [[7, 1],
                                     #      [1, 6]]
E += 4*dx*x + 6*dx*y + 3*dy*x + 5*dy*y  # V = [[4, 6],
                                        #      [3, 5]]
E += (-3)*x**2 + (-4)*x*y + (-5)*y**2  # K = [[-3, -2],
                                       #      [-2, -5]]
E += 1*dx + (-7)*dy  # A = [1, -7]
E += (-5)*x + 9*y  # B = [-5, 9]
E += 80  # C = 80

现有实现与问题 #

目前实现思路是手动求导并逐次减去对应项，以下代码在表达式仅含常量系数时运行正常：

diff = sp.derive_by_array
def dot(A, B):
    # 实现A和B的矩阵乘法
    ndimA = len(sp.shape(A))
    C = sp.tensorproduct(A, B)
    D = sp.tensorcontraction(C, (ndimA-1, ndimA))
    return D

M = diff( diff(E, dX), dX)/2
E -= dot(dX, dot(M, dX))
V = diff( diff(E, X), dX)
E -= dot(dX, dot(V, X))
K = diff( diff(E, X), X)/2
E -= dot(X, dot(K, X))
A = diff(E, dX)
E -= dot(A, dX)
B = diff(E, X)
E -= dot(B, X)
C = sp.expand(E)

但当传入非常量项时结果错误，比如E = cos(x)的情况：

# 预期输出
M = [[0, 0],
     [0, 0]]
V = [[0, 0],
     [0, 0]]
K = [[0, 0],
     [0, 0]]
A = [0, 0]
B = [0, 0]
C = cos(x)

# 实际输出
M = [[0, 0],
     [0, 0]]
V = [[0, 0],
     [0, 0]]
K = [[-cos(x)/2, 0],
     [0, 0]]
A = [0, 0]
B = [-x(t)**2*sin(x(t))/2 + x(t)*cos(x(t)) - sin(x(t)), 0]
C = x(t)**3*sin(x(t))/2 - x(t)**2*cos(x(t))/2 + x(t)*sin(x(t)) + cos(x(t))

我尝试过分解E的表达式结构，但是希望存在更简便的实现方案。

内容的提问来源于stack exchange，提问作者Carlos Adir