# 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[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,

more

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

( 2020-07-16 16:32:30 +0100 )edit
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...).

( 2020-07-16 17:34:32 +0100 )edit

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