Существует функция hessian
для выражений и метод jacobian
для матриц.
Вот функция и переменные вашей задачи:
>>> from sympy.abc import x, y
>>> from sympy import ordered, Matrix, hessian
>>> eq = x**2/2 + 5*y**2 + 2*(x - 2)**4/3 + 8*(y + 1)**4
>>> v = list(ordered(eq.free_symbols)); v
[x, y]
Мы можем написать собственный помощник для градиента, который создаст матрицу и использует для нее метод jacobian
:
>>> gradient = lambda f, v: Matrix([f]).jacobian(v)
Тогда величины можно рассчитать как:
>>> gradient(eq, v)
Matrix([[x + 8*(x - 2)**3/3, 10*y + 32*(y + 1)**3]])
>>> hessian(eq, v)
Matrix([
[8*(x - 2)**2 + 1, 0],
[ 0, 96*(y + 1)**2 + 10]])