Function to calculate minima and maxima of a function using Lagrange Multipliers

[hamorabi]

Sometimes, we are interested in finding the minima and maximum of a function ff subjected to a constraint by one or more constraints, usually some curve(s) along the surface of f.

The function lagrangemult takes 3 parameters:
1. ff: this is the function f(x,y,z) such that f(x,y,z)=0
2. cnt: this is/are the constraint(s) to which ff is subjected to. The constraints should be entered separately in a list and defined as g(x,y,z)=0
3. vars: These are the variables used.

Here are a couple of examples, mostly taken from this page:

PHP Code:

lagrangemult(5x-3y,[x^2+y^2-136],[x,y])
lagrangemult(x^2+y^2,[x^2+y^2+2x-2y+1],[x,y])
lagrangemult(x*y*z,[x+y+z-1],[x,y,z]) 


With 2 constraints:

PHP Code:

lagrangemult(4y-2z,[2x-y-z-2, x^2+y^2-1],[x,y,z]) 


The result is presented in a matrix format. The first j columns correspond to the values of of the variables. For instance, column 1 is for the x-coordinate, column 2 for y, etc. The last column is the function ff evaluated at the point P(x,y,z). Each row corresponds to either a minima or a maxima. For example, if there are 4 rows, then there can be 3 local minima and 1 maxima, or some other mixture.


The function works well in most cases except in those where the solver fails to find a solution even though a solution exists. I have attached a screenshot that demonstrates this. You can clearly see that the function arrives at the same equations as in the example. However, for a reason which I don't understand at the moment, the solver fails. If someone has any idea why this is happening or what should be done to prevent it, I'd be happy to adjust the code.

This is the problematic case. Also see attached image:

PHP Code:

lagrangemult(x^2+2y^2+z,[x^2-z^2-2],[x,y,z]) 


Finally, the code:

PHP Code:

#cas
def lagrangemult(ff,cnt,vars):
    lx = len(vars)
    ly = lx + len(cnt)
    λλ = seq(expr("λ" + str(ii)),ii,0,len(cnt)-1)
    vars = concat(vars,λλ)
    hh = sum(λλ[ii]*cnt[ii],ii,0,len(cnt)-1) + ff
    ss = zeros(seq(equal(diff(hh,vars[ii]),0),ii,0,ly-1),vars)[0:ly,0:lx+1]
    for ii in range(0, rowDim(ss)):
        ss[ii][lx] = subst(ff,seq(equal(vars[jj],ss[ii][jj]),jj,0,lx-1))
    return (ss)
#end