from sympy import symbols, Function, Eq, Derivative, Rational, UnevaluatedExpr, solve
import xv.km

x, y, z, lamda = symbols('x, y, z, lamda')
k = symbols('k', constant = True)
f = Function('f')
F = Function('F')
g = Function('g')

eq1 = Eq(F(x, y, z, lamda), f(x, y, z) + lamda*(g(x, y, z) - k))
eq1

$\displaystyle F{\left(x,y,z,\lambda \right)} = \lambda \left(- k + g{\left(x,y,z \right)}\right) + f{\left(x,y,z \right)}$

F_x = Eq(Eq(Derivative(eq1.lhs, x), Derivative(eq1.rhs, x)), 0, evaluate = False)
F_x

$\displaystyle \frac{\partial}{\partial x} F{\left(x,y,z,\lambda \right)} = \frac{\partial}{\partial x} \left(\lambda \left(- k + g{\left(x,y,z \right)}\right) + f{\left(x,y,z \right)}\right) = 0$

F_y = Eq(Eq(Derivative(eq1.lhs, y), Derivative(eq1.rhs, y)), 0, evaluate = False)
F_y

$\displaystyle \frac{\partial}{\partial y} F{\left(x,y,z,\lambda \right)} = \frac{\partial}{\partial y} \left(\lambda \left(- k + g{\left(x,y,z \right)}\right) + f{\left(x,y,z \right)}\right) = 0$

F_z = Eq(Eq(Derivative(eq1.lhs, z), Derivative(eq1.rhs, z)), 0, evaluate = False)
F_z

$\displaystyle \frac{\partial}{\partial z} F{\left(x,y,z,\lambda \right)} = \frac{\partial}{\partial z} \left(\lambda \left(- k + g{\left(x,y,z \right)}\right) + f{\left(x,y,z \right)}\right) = 0$

F_lamda = Eq(Eq(Derivative(eq1.lhs, lamda), Derivative(eq1.rhs, lamda)), 0, evaluate = False)
F_lamda

$\displaystyle \frac{\partial}{\partial \lambda} F{\left(x,y,z,\lambda \right)} = \frac{\partial}{\partial \lambda} \left(\lambda \left(- k + g{\left(x,y,z \right)}\right) + f{\left(x,y,z \right)}\right) = 0$

Example

Let x, y be real numbers. Find the maximum value of x + y, such that

eq2 = Eq(x**3 + y**3 + UnevaluatedExpr(x+y)/4 - Rational(15, 2), 0)
eq2

$\displaystyle x^{3} + y^{3} - \frac{15}{2} + \frac{x + y}{4} = 0$

eq3 = Eq(F(x, y, lamda), f(x, y) + lamda*(g(x, y) - k))
eq3

$\displaystyle F{\left(x,y,\lambda \right)} = \lambda \left(- k + g{\left(x,y \right)}\right) + f{\left(x,y \right)}$

f_eq = Eq(f(x, y), x + y)
f_eq

$\displaystyle f{\left(x,y \right)} = x + y$

g_eq = Eq(g(x, y), eq2.lhs)
g_eq

$\displaystyle g{\left(x,y \right)} = x^{3} + y^{3} - \frac{15}{2} + \frac{x + y}{4}$

eq4 = eq3.subs({f_eq.lhs: f_eq.rhs, g_eq.lhs: g_eq.rhs, k: 0})
eq4

$\displaystyle F{\left(x,y,\lambda \right)} = \lambda \left(x^{3} + y^{3} - \frac{15}{2} + \frac{x + y}{4}\right) + x + y$

F_x = Eq(Derivative(eq4.rhs, x), 0, evaluate = False).doit()
F_x

$\displaystyle \lambda \left(3 x^{2} + \frac{1}{4}\right) + 1 = 0$

F_y = Eq(Derivative(eq4.rhs, y), 0, evaluate = False).doit()
F_y

$\displaystyle \lambda \left(3 y^{2} + \frac{1}{4}\right) + 1 = 0$

F_lamda = Eq(Derivative(eq4.rhs, lamda), 0, evaluate = False).doit()
F_lamda

$\displaystyle x^{3} + \frac{x}{4} + y^{3} + \frac{y}{4} - \frac{15}{2} = 0$

solns = solve((F_x, F_y, F_lamda), x, y, lamda, dict = True)
solns

[{x: 3/2, y: 3/2, lamda: -1/7},
 {x: -3/4 + sqrt(31)*I/4,
  y: -3/4 + sqrt(31)*I/4,
  lamda: 1/14 - 9*sqrt(31)*I/434},
 {x: -3/4 - sqrt(31)*I/4,
  y: -3/4 - sqrt(31)*I/4,
  lamda: 1/14 + 9*sqrt(31)*I/434}]

for sol in solns:
    for var, val in sol.items():
        display(Eq(var, val))

$\displaystyle x = \frac{3}{2}$

$\displaystyle y = \frac{3}{2}$

$\displaystyle \lambda = - \frac{1}{7}$

$\displaystyle x = - \frac{3}{4} + \frac{\sqrt{31} i}{4}$

$\displaystyle y = - \frac{3}{4} + \frac{\sqrt{31} i}{4}$

$\displaystyle \lambda = \frac{1}{14} - \frac{9 \sqrt{31} i}{434}$

$\displaystyle x = - \frac{3}{4} - \frac{\sqrt{31} i}{4}$

$\displaystyle y = - \frac{3}{4} - \frac{\sqrt{31} i}{4}$

$\displaystyle \lambda = \frac{1}{14} + \frac{9 \sqrt{31} i}{434}$

x_val = solve(F_x, x)[0]
x_val

$\displaystyle \frac{\sqrt{3} \sqrt{\frac{- \lambda - 4}{\lambda}}}{6}$

y_val = solve(F_y, y)[0]
y_val

$\displaystyle \frac{\sqrt{3} \sqrt{\frac{- \lambda - 4}{\lambda}}}{6}$

lamda_val = solve(F_lamda.subs({x: x_val, y: y_val}), lamda)[0]
lamda_val

$\displaystyle - \frac{1}{7}$

x_val.subs(lamda, lamda_val).simplify(force = True)

$\displaystyle \frac{3}{2}$

y_val.subs(lamda, lamda_val).simplify(force = True)

$\displaystyle \frac{3}{2}$

def find_lagrange_multipliers(f_eq, g_eq, k = 0):
    f_eq = Eq(f(x, y), f_eq)
    g_eq = Eq(g(x, y), g_eq)
    eq3 = Eq(F(x, y, lamda), f(x, y) + lamda*(g(x, y) - k))
    eq4 = eq3.subs({f_eq.lhs: f_eq.rhs, g_eq.lhs: g_eq.rhs,})
    F_x = Eq(Derivative(eq4.rhs, x), 0, evaluate = False).doit()
    F_y = Eq(Derivative(eq4.rhs, y), 0, evaluate = False).doit()
    F_lamda = Eq(Derivative(eq4.rhs, lamda), 0, evaluate = False).doit()
    solns = solve((F_x, F_y, F_lamda), x, y, lamda, dict = True)
    dict_solns = {x: [], y: []}
    for sol in solns:
        dict_solns[x].append(sol[x])
        dict_solns[y].append(sol[y])
        
    return dict_solns

f_eq = x + y
g_eq = x**3 + y**3 + UnevaluatedExpr(x+y)/4 - Rational(15, 2)

solns = find_lagrange_multipliers(f_eq, g_eq)
solns

{x: [3/2, -3/4 + sqrt(31)*I/4, -3/4 - sqrt(31)*I/4],
 y: [3/2, -3/4 + sqrt(31)*I/4, -3/4 - sqrt(31)*I/4]}

for i, j in solns.items():
    print('The values of', i)
    for val in j:
        display(Eq(i, val, evaluate = False))
    print('_'*50)

The values of x

$\displaystyle x = \frac{3}{2}$

$\displaystyle x = - \frac{3}{4} + \frac{\sqrt{31} i}{4}$

$\displaystyle x = - \frac{3}{4} - \frac{\sqrt{31} i}{4}$

__________________________________________________
The values of y

$\displaystyle y = \frac{3}{2}$

$\displaystyle y = - \frac{3}{4} + \frac{\sqrt{31} i}{4}$

$\displaystyle y = - \frac{3}{4} - \frac{\sqrt{31} i}{4}$

__________________________________________________