# Solve() not working for Expenditure Minimization Problem in Sage

Though I've figured out how to go through a standard utility maximization problem analytically in sage Im having difficulty solving its dual problem, expenditure minimization. The code I have is the following:

x1, 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);
solve([dLdx == 0, dLdy == 0, dLdl == 0], x1, x2, l)


The output I get however is:

[-b*l*x1^a*x2^(b - 1) + p2 == 0, -a*l*x1^(a - 1)*x2^b + p1 == 0, -x1^a*x2^b + R == 0]


which are only the first order conditions and not solving for $x_1$ and $x_2$.

This code is exactly like my utility maximization code here however I only moved where m and U are.

why is this the case?

Note: if I were to consider the following code:

x1, x2, l, p1, p2, a, b,  R= var('x1, x2, l, p1, p2,  R')
U = x1^1/3*x2^2/3
m = p1*x1+p2*x2;
L = m+ l * (R-U);
dLdx = L.diff(x1);
dLdy = L.diff(x2);
dLdl = L.diff(l);
solve([dLdx == 0, dLdy == 0, dLdl == 0], x1, x2, l)


I obtain the following analytical solution:

[[x1 == 1/2*18^(1/3)*R/(R*p1/p2)^(2/3), x2 == 18^(1/3)*(R*p1/p2)^(1/3), l==1/2*18^(1/3)*p1/(R*p1/p2)^(2/3)],
[x1 == 1/4*18^(1/3)*R*(-I*sqrt(3) - 1)/(R*p1/p2)^(2/3), x2 == -1/2*18^(1/3)*
(R*p1/p2)^(1/3)*(I*sqrt(3) + 1), l == 1/4*18^(1/3)*p1*(-I*sqrt(3) - 1)/(R*p1/p2)^(2/3)],
[x1 == -1/4*18^(1/3)*R*(-I*sqrt(3) + 1)/(R*p1/p2)^(2/3), x2 == 1/2*18^(1/3)*(R*p1/p2)^(1/3)*(I*sqrt(3) - 1), l == -1/4*18^(1/3)*p1*(- I*sqrt(3) + 1)/(R*p1/p2)^(2/3)]]


My guess is that there is a bug with regards to the rules of multiplying undefined exponents.

edit retag close merge delete

Sort by » oldest newest most voted

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.lhs():S3.rhs().log().log_expand().factor().exp().log_expand()}
Sol.update({S2.lhs():S2.subs(Sol).rhs().log().log_expand().factor().exp()})
Sol.update({S1.lhs():S1.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:u.sage() for u in MSol}
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,

more

I'm having a bit of difficulty reading the unformatted code in the middle there can you just type set it?

2

My bad (a forgotten newline...). The markup used by ask.sagemath.org is not the most streamlined on the market. I'd rather use Github-flavored markdown (for which pandoc or org-mode offer convenient solutions...).

Thanks to Juanjo's answer to my question here I've managed to get a script which gives us an analytical solution for our hicksian demand equations. This is:

x1, 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);
step1=solve([dLdx == 0, dLdy == 0, dLdl == 0], p1, p2, R)
step2=solve((p1/p2).subs(step1)==p1/p2, x1)
good2step3 = solve((U/R).subs(step2).log().log_expand()==0, x2)


where good1step4 is our hicksiand demand for $x_1$ and good2step3 is our hicksian demand for $x_2$.

Its interesting how you need a different procedure when your constraint is non-linear.

more