Maxima: unable to solve Lagrangian

Question

I think this is may be more of a technical issue than a mathematical issue. I don't know...

I have the following function I want to maximize:

$$\arg\max_{x,y}\ U(x,y) = \frac{1-e^{-x}}{2} + \frac{1-e^{-y}}{2}$$

subject to $x+y \leq 50$ and $x,y\geq0$.

(My actual poblem is more complicated than this of course, but take this as a MWE please. :))

We can easily see this has a maximum. Using Maxima:

U(x,y) := (1-exp(-x))/2 + (1-exp(-y))/2;
wxplot2d([U(x,50-x)], [x,0,50]);

To solve this using Lagrangian optimzation along the $x+y=50$ frontier, I would do:

$$\begin{cases} \dfrac{\partial U}{\partial x} = \lambda \\[4pt] \dfrac{\partial U}{\partial y} = \lambda \\[4pt] x+y = 50 \end{cases} $$

We can easily solve this. $\frac{\partial U}{\partial x} = \frac{e^{-x}}{2}=\lambda \iff \frac{e^{-x}}{2} = \frac{e^{-50+x}}{2} \iff x=25$, as could be seen in the graph anyway.

But now I try to do this in Maxima:

U(x,y) := (1-exp(-x))/2 + (1-exp(-y))/2;
solve([diff(U(x,y),x)=lambda, diff(U(x,y),y)=lambda, x+y=50], [x,y,lambda]);

And it says there is no solution:

[]

What gives? I haven't done anything in optimization or Maxima for a long time... But I cannot for the life of me discover what the problem is...

EDIT: I would prefer to know where my Maxima usage is wrong (or my mathematics), but I will approve any answers that suggest me to use another symbolic solver.

Answer 1

Your solution is correct. It is also the same in Mathematica. One can still use "Reduce" in Mathematica as follows:

Reduce[{E^-x/2 - a == 0, E^-y/2 - a == 0, x + y == 50}, {x, y, a}]

with the following answer:

C[1] \[Element] Integers && (x == 25 + I \[Pi] + 2 I \[Pi] C[1] || x == 25 + 2 I \[Pi] C[1]) && y == 50 - x && a == E^-x/2

Since the problem is simple one can solve this problem without using any Lagrangian multiplier as well.

The function $U$ is increasing both in $x$, $y$ and in both. So this means for the maximization $x$ and $y$ should be as large as possible. So this means the constraint $x+y\leq 50$ should hold with equality, that is $x+y=50$. Then, one can simply insert $x:=-y+50$ in $U$, take the derivative w.r.t. $y$ and see that $y$ must be $25$. Accordingly $x=25$.

You may think about asking questions about programming in mathematica.stackexchange

Answer 2

In Maple:

U:= (1-exp(-x))/2 + (1-exp(-y))/2;
solve({diff(U,x)=lambda, diff(U,y)=lambda, x+y=50},{lambda,x,y});

$$ \left\{ \lambda=\frac{1}{2 {\rm e}^{25}},x=25,y=25 \right\} $$

Answer 3

For future reference, I was able to solve this using Sympy, though I had to get rid of the $\lambda$... (If someone knows why please let me know.)

from sympy import *
x, y = symbols('x y')

U = (1-exp(-x))/2 + (1-exp(-y))/2
dU_dx = diff(U, x)
dU_dy = diff(U, y)

print solve((Eq(dU_dx,dU_dy), Eq(x+y,50)), (x, y))