Revision history [back]

It cannot be done "automagically" (neither by Maxima's solver (default) nor Sympy's), but "manually" solving each equation at a time, high school-style, works.

Elements of the problem :

1, x2, l, p1, p2, a, b,  R= var('x1, x2, l, p1, p2, a, b,  R')
U = x1^a*x2^b
m = p1*x1+p2*x2
L = m+ l * (R-U)
dLdx = L.diff(x1)
dLdy = L.diff(x2)
dLdl = L.diff(l)

Equations to solve :

E1=dLdx==0
E2=dLdy==0
E3=dLdl==0

Solve one-by-one, thus eliminating all unknowns (with a little manual help):

S1=((E1-p1).log().log_expand()-I*pi).solve(x1)
S2=((E2.subs(S1)-p2).log().log_expand()-I*pi).solve(x2)
S3=((E3.subs(S1).subs(S2)-R).log().log_expand()-I*pi).solve(l)

Substitute each unknown in partial solutions (and clean up a bit) : Sol={S3[0].lhs():S3[0].rhs().log().log_expand().factor().exp().log_expand()} Sol.update({S2[0].lhs():S2[0].subs(Sol).rhs().log().log_expand().factor().exp()}) Sol.update({S1[0].lhs():S1[0].subs(Sol).rhs().log().log_expand().factor().exp()})

So we get :

sage: Sol
{l: e^(-(a*log(R) + b*log(R) + a*log(a) + b*log(b) - a*log(p1) - b*log(p2) - log(R))/(a + b)),
 x2: e^(-(a*log(a) - a*log(b) - a*log(p1) + a*log(p2) - log(R))/(a + b)),
 x1: e^((b*log(a) - b*log(b) - b*log(p1) + b*log(p2) + log(R))/(a + b))}

If you happen to have access to Mathematica, you can check your solution against it :

MSol=mathematica.Solve([E1,E2,E3],[x1,x2,l])
sage: MSol
{{x1 -> E^((b*Log[a] - b*Log[b] - b*Log[p1] + b*Log[p2] + Log[R])/(a + b)), 
  x2 -> E^((-(a*Log[a]) + a*Log[b] + a*Log[p1] - a*Log[p2] + Log[R])/
     (a + b)), l -> E^((-(a*Log[a]) - b*Log[b] + a*Log[p1] + b*Log[p2] + 
      Log[R] - a*Log[R] - b*Log[R])/(a + b))}}

Which can be translated as :

MST={u[1]:u[2].sage() for u in MSol[1]}
sage: MST
{x1: e^((b*log(a) - b*log(b) - b*log(p1) + b*log(p2) + log(R))/(a + b)),
 x2: e^(-(a*log(a) - a*log(b) - a*log(p1) + a*log(p2) - log(R))/(a + b)),
 l: e^(-(a*log(R) + b*log(R) + a*log(a) + b*log(b) - a*log(p1) - b*log(p2) - log(R))/(a + b))}

Same difference...

HTH,