Computing Jacobian and Hessian

from sympy import * import pprint pp = pprint.PrettyPrinter(indent=4) def f(x, y, z): return x**3 - y**2 + z**2/(2**0.5) + (x*y - z*x) / y*z def g(x, y, z): return z*2 - (x*y)**0.5 / x*y*z**2 def h(x, y, z): return 14*x + 7*z*y / (x**3*y) x, y, z = symbols('x y z') flabels = 'f g h'.split() ff, fg, fh = f(x, y, z), g(x, y, z), h(x, y, z) fs, syms = [ff, fg, fh], [x, y, z] jac = [[[] for i in range(3)] for j in range(3)] # Hessians exist for each function hessians = [[[[] for i in range(3)] for j in range(3)] for k in range(3)] for i in range(3): for j in range(3): jac[i][j] = diff(fs[i], syms[j]) for k in range(3): hessians[i][j][k] = diff(jac[i][j], syms[k]) print('The Jacobian:') pp.pprint(jac) print('\nThe Hessians (for each function):') for i, hes in enumerate(hessians): print('of %s:' % (flabels[i])) pp.pprint(hes) print('\n') """ The Jacobian: [ [ 3*x**2 + z*(y - z)/y, x*z/y - 2*y - z*(x*y - x*z)/y**2, -x*z/y + 1.41421356237309*z + (x*y - x*z)/y], [ 0.5*y*z**2*(x*y)**0.5/x**2, -1.5*z**2*(x*y)**0.5/x, 2 - 2*y*z*(x*y)**0.5/x], [14 - 21*z/x**4, 0, 7/x**3]] The Hessians (for each function): of f: [ [6*x, z/y - z*(y - z)/y**2, -z/y + (y - z)/y], [ z/y - z*(y - z)/y**2, -2*x*z/y**2 - 2 + 2*z*(x*y - x*z)/y**3, x/y + x*z/y**2 - (x*y - x*z)/y**2], [ -z/y + (y - z)/y, x/y + x*z/y**2 - (x*y - x*z)/y**2, -2*x/y + 1.41421356237309]] of g: [ [ -0.75*y*z**2*(x*y)**0.5/x**3, 0.75*z**2*(x*y)**0.5/x**2, 1.0*y*z*(x*y)**0.5/x**2], [ 0.75*z**2*(x*y)**0.5/x**2, -0.75*z**2*(x*y)**0.5/(x*y), -3.0*z*(x*y)**0.5/x], [1.0*y*z*(x*y)**0.5/x**2, -3.0*z*(x*y)**0.5/x, -2*y*(x*y)**0.5/x]] of h: [[84*z/x**5, 0, -21/x**4], [0, 0, 0], [-21/x**4, 0, 0]] """

If you notice a mistake, feel free to contribute.