# MAPLE ASSINMENT 11
# CURVATURE OF SURFACES
# The domains of the following functions are restricted to {(x,y)/x>0 &
# y>0}.
> U:=x^(1/4)*y^(1/2);
> plot3d(U,x=0..5,y=0..10-2*x);
# Manipulate the plot to convince yourself that this surface is concave.
# Then go back and change (1/4,1/2) to (1,3/4)
# or to anything where the sum of the 2 powers is > 1. The resulting
# surface is neither concave or convex. It is however quasiconcave which
# may be clearer to you if you again make some changes in the appearance
# of the plot. The next commands show  with calculus that the surface
# is not concave.
> with(linalg):
> HU:=hessian(U,[x,y]);
> det(HU);
# The 0 diagonal element in the Hessian rules out strict concavity.  The
# negative det. of the Hessian rules out concavity completely.  Now go
# back to the exponents (1/4,1/2) in U and try the above again.  a11 in
# the Hessian is negative and the det. of the Hessian is
# positive so the function is strictly concave. The next surface looks
# like it might be strictly convex but it is a little difficult to tell.
# From the Hessian matrix you can see
# that it is. 
> V:=sqrt(x+y)/(x*y);
> plot3d(V,x=1..5,y=2..12-2*x);
> HV:=hessian(V,[x,y]);
> simplify(HV[1,1]);
> det(HV);
# Maple appears not to have heard of quasiconcavity but a little trick
# can be used 
# to construct the bordered Hessian. This is done next for the
# utility function whose graph you plotted earlier. If you do not
# understand the trick, ask your TA to explain it.
> U:=x*y^(3/4);
> V:=lambda*U;
> HV:=hessian(V,[lambda,x,y]);
> HUb:=subs(lambda=1,%);
> det(HUb);
# Since the determinant of the bordered Hessian is positive, U is
# strictly quasiconcave. This can also be 
# seen by looking at the curvature of the level curves.
> with(plots):
> contourplot(U,x=0..5,y=0..5);
# UNCONSTRAINED OPTIMIZATION OF A 2 VARIABLE FUNCTION
# The following example illustrates the principles of unconstrained
# optimization
# for a 2 variable function. 
> f:=x^3-3*x*y+y^3;
> with(linalg):
> gradf:=grad(f,[x,y]);
> statpt:=solve({gradf[1],gradf[2]},{x,y});
# The symbolic solve command finds the 2 stationary points and 1 garbage
# solution. The numeric solver command fsolve could also be used but it
# finds only 1 stationary point. The subs command can now be used to 
# check
# that Maple's solutions are correct.
> subs(statpt[1],gradf[1]);
> subs(statpt[1],gradf[2]);
# Change 1 to 2 to check the other stationary point. Now check the
# second order conditions to see what type each stationary point is.
> Hf:=hessian(f,[x,y]);
> subs(statpt[1],Hf[1,1]);
> subs(statpt[1],det(Hf));
# The negative determinant of the Hessian means the stationary point is
# a saddle point. 
# Change 1 to 2 again to find the type of the second stationary point.
# Finally plot the graph of  f and manipulate it to convince yourself
# that Maple has given the right answers. 
> plot3d(f,x=-2..3,y=-2..3);
# THE PROFIT MAXIMIZING FIRM
# The commands below define a production function, check for positive
# marginal products, check for strict concavity and then plot the
# function. Manipulate the 3d plot so it looks like what you might
# expect.
> f:=(sqrt(x)+sqrt(y))^(3/2);
> with(linalg):
> gradf:=grad(f,[x,y]);
> Hf:=hessian(f,[x,y]);
> simplify(Hf[1,1]);
> det(Hf);
> with(plots): 
> contourplot(f,x=0..20,y=0..20);
> plot3d(f,x=0..20,y=0..20-x);
# The following calculations have p:=2;r:=1;w:=1.5;
> PR:=2*f-1*x-(3/2)*y;
> subs({x=1,y=2},PR);
> evalf(%);
> gradPR:=grad(PR,[x,y]);
> extpt:=solve({gradPR[1],gradPR[2]},{x,y});
> subs(extpt,gradPR[1]);
> evalf(%);
# Check that the second first order condition is also satisfied.
> plot3d(PR(x,y),x=0..30,y=0..30,orientation=[-15,45]);
# Profit functions can be very flat but there is indeed a maximum within
# the highest contour line which you should try to see better by
# manipulating the plot..
> contourplot(PR(x,y),x=0..50,y=0..30);
> PRmax:=subs(extpt,PR);
> evalf(%);
# You can try to do better by substituting your own numbers for the ?s
# below. Then check the second order conditions for a max and show that
# extpt is where an isoquant is tangent to an isocost line.
> evalf(subs({x=?,y=?},PR));
> HPR:=hessian(PR,[x,y]);
> evalf(subs(extpt,HPR[1,1]));
> evalf(subs(extpt,det(HPR)));
> extq:=evalf(subs(extpt,f));
> extC:=subs(extpt,1*x+(3/2)*y);
> implicitplot({f=extq,1*x+(3/2)*y=extC},x=0..30,y=0..30);
# Click on the tangency point to see if you get close to extpt.
# Now you can  get Maple to solve for the input demand functions where
# r has been normalized to 1.
> PR:=p*f-x-w*y;
> gradPR:=grad(PR,[x,y]);
> dem:=solve({gradPR[1],gradPR[2]},{x,y});
> subs(dem,gradPR[1]);
> simplify(%);
# This actually does simplify to 0 if you simplify further by hand but
# Maple doesn't recognize this. Now you can plot the demand for y in the
# conventional way as a function of its price w. To do this, it seems to
# be helpful to convert dem  to a function and then assign a fixed
# value to p when plotting. The output price in this context is called a
# shifter of the demand curve.
> yd:=subs(dem,y);
> ydem:=unapply(yd,p,w);
> plot(ydem(3,w),w=1..2);
> plot({ydem(2,w),ydem(3,w),ydem(4,w)},w=1..2);
# The demand curves for y slope down as expected. Which way does an
# increase in p shift the demand curve? Next you can plot the 2 variable
# demand surface. Does the shape of the demand surface conform to your
# expectations? (manipulate the plot to make it more understandable).
# You can ask your TA to explain this plot to you if you wish.  
> contourplot(ydem(p,w),p=3..6,w=1..2);
> plot3d(ydem(p,w),p=3..6,w=1..2);
# Now the firm's supply curve can be defined and plotted. Now w is a
# shifter of the supply curve.
> sup:=subs(dem,f);
> qsup:=unapply(sup,p,w);
> qsup(3,1);
> plot(qsup(p,1),p=3..6);
> plot({qsup(p,1),qsup(p,2),qsup(p,3)},p=3..6);
# Which of these supply curves corresponds to each of the w values?
> contourplot(qsup(p,w),p=3..6,w=1..2);
> plot3d(qsup(p,w),p=3..6,w=1..2);
# Does the supply surface conform to your expectations about supply?
# Finally you can define and plot the firm's maximized profit function.
> PRm:=subs(dem,PR);
> evalf(subs({p=3,w=1},PRm));
> PRmax:=unapply(PRm,p,w);
> plot(PRmax(p,1),p=3..6);
> plot({PRmax(p,1),PRmax(p,2),PRmax(p,3)},p=3..6);
> contourplot(PRmax(p,w),p=3..6,w=1..2);
> plot3d(PRmax(p,w),p=3..6,w=1..2,orientation=[-15,45]);
# Which way did an increase in w shift the maximized profit curve? Does
# the maximized profit surface conform to your expectations?
> 
