Constrained optimization with strange result.
cyrille piatecki
assume(a,'real')
assume(x,'real')
assume(y,'real')
assume(p_x,'real')
assume(p_y,'real')
assume(D,'real')
assume(Rev,'real')
assume(R,'real')
assume(l,'real')
assume(a<1)
assume(a>0)
assume(p_x>0)
assume(p_y>0)
assume(R>0)
U =(1/(a+1))*x^(a+1)+y
show(LatexExpr(r'\text{La fonction d}^\prime\text{utilité est }U(x,y) = '),U)
D= x*p_x + y*p_y
show(LatexExpr(r'\text{La Dépense } D = '),D)
Rev= R
show(LatexExpr(r'\text{Le Revenu } Rev = '),R)
L=U+l*(Rev-D)
show(LatexExpr(r'\text{Le lagrangien est } \mathcal{L}(x, y, λ) = '),L)
FOC = [diff(L,x),diff(L,y),diff(L,l)]
show(LatexExpr(r'\text{Les condition du premier ordre sont } \left\{\begin{array}{c}\mathcal{L}_x= 0\\\mathcal{L}_y= 0\\\mathcal{L}_λ= 0\end{array}\right. '))
show(LatexExpr(r'\text{soit }'))
show(LatexExpr(r'\mathcal{L}_x= 0 \Longleftrightarrow '),FOC[0]==0)
show(LatexExpr(r'\mathcal{L}_y= 0 \Longleftrightarrow '),FOC[1]==0)
show(LatexExpr(r'\mathcal{L}_λ= 0 \Longleftrightarrow '),FOC[2]==0)
sol=solve([FOC[0]==0,FOC[1]==0,FOC[2]==0],x,y,l)
show(sol)
Emmanuel Charpentier
Variables of the form z_xxxx are integer variables created by Maxima, which attempts to give you also the complex roots, if any, thus ignoring the assumptions on x, y and l. Note that :
sage: solve(FOC[0], x)
---------------------------------------------------------------------------
[ Snip… ]
TypeError: Computation failed since Maxima requested additional constraints; using the 'assume' command before evaluation *may* help (example of legal syntax is 'assume(l>0)', see `assume?` for more details)
Is l positive, negative or zero?
sage: with assuming(l>0): print(solve(FOC[0], x))
[
x == (l*p_x)^(1/a)
]
sage: with assuming(l<0): print(solve(FOC[0], x))
[
x^a == l*p_x
]
sage: with assuming(l<0): print(solve(FOC[0], x, to_poly_solve=True))
[x == (l*p_x)^(1/a)*e^(2*I*pi*z4353/a)]
Interestingly:
sage: solve([FOC[0]==0,FOC[1]==0,FOC[2]==0],x,y,l)
[[l == (1/p_y), x == (p_x/p_y)^(1/a)*e^(2*I*pi*z3540/a), y == -(p_x*(p_x/p_y)^(1/a)*e^(2*I*pi*z3540/a) - R)/p_y]]
sage: solve([FOC[0]==0,FOC[1]==0,FOC[2]==0],x,y,l, algorithm="sympy")
[{x: (p_x/p_y)^(1/a), l: 1/p_y, y: -(p_x*(p_x/p_y)^(1/a) - R)/p_y}]
sage: solve([FOC[0]==0,FOC[1]==0,FOC[2]==0],x,y,l, algorithm="fricas")
[[l == (1/p_y), x == (p_x/p_y)^(1/a)*e^(2*I*pi*z3891/a), y == -(p_x*(p_x/p_y)^(1/a)*e^(2*I*pi*z3891/a) - R)/p_y]]
sage: giac.solve(giac(FOC),giac([x,y,l])).sage()
[[(p_x/p_y)^(1/a), -(p_x*(p_x/p_y)^(1/a) - R)/p_y, 1/p_y]]
HTH,
cyrille....@univ-orleans.fr
De: Emmanuel Charpentier <emanuel.c...@gmail.com>
À: sage-support <sage-s...@googlegroups.com>
Envoyé: Mon, 29 Nov 2021 11:03:37 +0100 (CET)
Objet: [sage-support] Re: Constrained optimization with strange result.
Emmanuel Charpentier
Le lundi 29 novembre 2021 à 16:33:41 UTC+1, cyrille piatecki a écrit :
Thanks Emmanuel for your precious answer. But It generates some few new questions :- is there a place in the documentation where I can find the information on `solve()` and mainly its options ?
The documentation, of course…
- if I understand clearly z_{6497} is an integer but how to fix it to zero --- when the number change at each iteration
That’s why I used a methods sequence to designate it, rather than using its name…
- sympy seems to be the good approach
Beware : see below…
but it is not self evident that to call y one must typpeset sol2[0][x]
It is, because algorithm="sympy" will cause the results to be expressed as dictionaries and D[x] is the canonical way to get the value of the entry of dictionary D indexed by x. Basic Python…
- the giac way is certainly the better but it keeps no track of the variable's order.
Again, ask for a solution dictionary. As for algorithm="giac", I have seen it go pear-shaped a couple times…
Now for the various expression of solutions : consider :
print(table([[u,solve(FOC,[x,y,l], solution_dict=True, algorithm=u)] for u in ["maxima", "fricas", "sympy", "giac"]]))
maxima [{l: 1/p_y, x: (p_x/p_y)^(1/a)*e^(2*I*pi*z3541/a), y: -(p_x*(p_x/p_y)^(1/a)*e^(2*I*pi*z3541/a) - R)/p_y}]
fricas [{l: 1/p_y, x: (p_x/p_y)^(1/a)*e^(2*I*pi*z3892/a), y: -(p_x*(p_x/p_y)^(1/a)*e^(2*I*pi*z3892/a) - R)/p_y}]
sympy [{x: (p_x/p_y)^(1/a), l: 1/p_y, y: -(p_x*(p_x/p_y)^(1/a) - R)/p_y}]
giac [{x: (p_x/p_y)^(1/a), y: -(p_x*(p_x/p_y)^(1/a) - R)/p_y, l: 1/p_y}]
Both maxima and fricas try to explicitly express the set of solutions of the equation z^a==p_x/p_y, which is a set of a complexes if a is a positive integer. (I leave to you (as en exercise ;-) to determine what it means (if any…) if a is rational, algebraic or transcendental, real or complex…).
OTOH, both sympy and giac use the notation (p_x/p_y)^(1/a) to implicitly denote *the very same set of solutions to the very same equation. One could say that tey are glossing over whatever maxima and fricas insist on. Choose your poison…
HTH,
cyrille....@univ-orleans.fr
"maxima", "fricas", "sympy", "giac"De: Emmanuel Charpentier <emanuel.c...@gmail.com>
À: sage-support <sage-s...@googlegroups.com>
Objet: Re: [sage-support] Re: Constrained optimization with strange result.
Emmanuel Charpentier
Thanks Emmanuel for your so precious answer. But obviously the doc must be upgraded. From"maxima", "fricas", "sympy", "giac"only "sympy" has a reference.
(Almost) right : the documentation I pointed to doesn’t list anything but maxima. The up to date documentation is similar, but gives example of algorithm="sympy" uses…
Would you care to file a ticket against that ? A further problem being that the documentatin mentions the use of algorithm for single-equation solving, whereas your current example shows success for (an) equation system. This should be clarified…
I wonder If one could add a parameter asking, in all cases, only for Real solutions.