Lagrange Multiplier

import numpy as np
import sympy as sp

Example

x, y, la= sp.symbols("x y la", real = True)

f = (x - 3)**2 + (y - 2)**2
f_num = sp.lambdify((x, y), f)
phi = (x / 3)**2 + (y / 2)**2 - 1
phi_num = sp.lambdify((x, y), phi)
L = f - la * phi
L_num = sp.lambdify((x, y, la), L)

L_grad = sp.Matrix([sp.diff(L, x), sp.diff(L, y), sp.diff(L, la)])
L_grad_num = sp.lambdify((x, y, la), list(L_grad)) # use list here to get a nicer output as flat array

H = sp.hessian(L, (x, y, la))
H_num = sp.lambdify((x, y, la), H)

Newton's Method

sol_old = np.array([2, 1, -2]) # for min
sol_new = np.array([2, 1, -2]) # for min

sol_old = np.array([-2, -1, -2]) # for max
sol_new = np.array([-2, -1, -2]) # for max
solutions = [sol_old]

for i in range(20):
    g_old = np.array(L_grad_num(*sol_old))
    H_old = np.array(H_num(*sol_old))
    sol_new = sol_old - np.linalg.solve(H_old, g_old) # better to avoid inverse and ask numpy to solve a system of linear equations directly
    solutions.append(sol_new)
    sol_old = sol_new
    
sol_new

array([-2.88143787, -0.55670291, 18.37032177])

np.sqrt(f_num(sol_new[0], sol_new[1]))

6.413114781967053

phi_num(sol_new[0], sol_new[1])

0.0