# MAPLE ASSIGNMENT 11
# The following illustrates the principles of unconstrained optimization
# for a 2 variable function. In this assignment you will work mainly
# with expressions rather than functions as Maple seems to work better
# with expressions.
> f:=x^3-3*x*y+y^3;

                               3            3
                         f := x  - 3 x y + y

> with(linalg):
Warning, new definition for norm
Warning, new definition for trace
> gradf:=grad(f,[x,y]);

                           [   2                  2]
                  gradf := [3 x  - 3 y, -3 x + 3 y ]

> extpt:=solve({gradf[1],gradf[2]},{x,y});

extpt := {x = 0, y = 0}, {x = 1, y = 1},

                  2                              2
    {x = RootOf(_Z  + _Z + 1), y = -1 - RootOf(_Z  + _Z + 1)}

# The symbolic solve command finds the 2 extreme points and 1 garbage
# solution. The numeric solver command fsolve could also be used but it
# finds only 1 extreme point. The subs command can now be used to  check
# that Maple's solutions are correct.
> subs(extpt[1],gradf[1]);

                                  0

> subs(extpt[1],gradf[2]);

                                  0

# Change 1 to 2 to check the other extreme point. Now check the second
# order conditions to see what type each extreme point is.
> Hf:=hessian(f,[x,y]);

                                [6 x    -3 ]
                          Hf := [          ]
                                [-3     6 y]

> subs(extpt[1],Hf[1,1]);

                                  6

> subs(extpt[1],det(Hf));

                                  27

# Change 1 to 2 again to find the type of the second extreme point.
# Finally plot the graph of f and manipulate it to convince yourself
# that Maple has given the right answers. 
> plot3d(f(x,y),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);

                               1/2    1/2 3/2
                        f := (x    + y   )

> gradf:=grad(f,[x,y]);

                 [      1/2    1/2 1/2        1/2    1/2 1/2]
                 [    (x    + y   )         (x    + y   )   ]
        gradf := [3/4 ----------------, 3/4 ----------------]
                 [           1/2                   1/2      ]
                 [          x                     y         ]

> Hf:=hessian(f,[x,y]);

Hf :=

    [                                1/2    1/2 1/2
    [             1                (x    + y   )
    [3/16 ------------------ - 3/8 ---------------- ,
    [       1/2    1/2 1/2                3/2
    [     (x    + y   )    x             x

                                   ]
                     1             ]
    3/16 --------------------------]
           1/2    1/2 1/2  1/2  1/2]
         (x    + y   )    y    x   ]

    [
    [                 1
    [3/16 -------------------------- ,
    [       1/2    1/2 1/2  1/2  1/2
    [     (x    + y   )    y    x

                                    1/2    1/2 1/2]
                 1                (x    + y   )   ]
    3/16 ------------------ - 3/8 ----------------]
           1/2    1/2 1/2                3/2      ]
         (x    + y   )    y             y         ]

> simplify(Hf[1,1]);

                                 1/2      1/2
                                x    + 2 y
                     - 3/16 ---------------------
                             3/2   1/2    1/2 1/2
                            x    (x    + y   )

> det(Hf);

                     2  3/2      3/2  2    5/2        5/2
                  3 x  y    + 3 x    y  + x    y + x y
            9/128 ---------------------------------------
                           1/2    1/2 2  5/2  5/2
                         (x    + y   )  x    y

> 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;

                           1/2    1/2 3/2
                 PR := 2 (x    + y   )    - x - 3/2 y

> subs({x=1,y=2},PR);

                             1/2     3/2
                         2 (2    + 1)    - 4

> evalf(");

                             3.502284402

> gradPR:=grad(PR,[x,y]);

             [      1/2    1/2 1/2            1/2    1/2 1/2      ]
             [    (x    + y   )             (x    + y   )         ]
   gradPR := [3/2 ---------------- - 1, 3/2 ---------------- - 3/2]
             [           1/2                       1/2            ]
             [          x                         y               ]

> extpt:=solve({gradPR[1],gradPR[2]},{x,y});

                                   225
                     extpt := {x = ---, y = 25/4}
                                   16

> subs(extpt,gradPR[1]);

                  1/2   1/2         1/2  1/2 1/2   1/2  1/2
    3/50 (1/16 225    16    + 1/4 25    4   )    25    4    - 3/2

> evalf(");

                                  0

# Check that the second first order conditon 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..30,y=0..30);
> PRmax:=subs(extpt,PR);

                           1/2   1/2         1/2  1/2 3/2   375
       PRmax := 2 (1/16 225    16    + 1/4 25    4   )    - ---
                                                            16

> evalf(");

                             7.812500000

# 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));

                              7.75804916

> HPR:=hessian(PR,[x,y]);

HPR :=

    [                               1/2    1/2 1/2
    [            1                (x    + y   )
    [3/8 ------------------ - 3/4 ---------------- ,
    [      1/2    1/2 1/2                3/2
    [    (x    + y   )    x             x

                                  ]
                    1             ]
    3/8 --------------------------]
          1/2    1/2 1/2  1/2  1/2]
        (x    + y   )    y    x   ]

    [
    [                1
    [3/8 -------------------------- ,
    [      1/2    1/2 1/2  1/2  1/2
    [    (x    + y   )    y    x

                                   1/2    1/2 1/2]
                1                (x    + y   )   ]
    3/8 ------------------ - 3/4 ----------------]
          1/2    1/2 1/2                3/2      ]
        (x    + y   )    y             y         ]

> assign(b);
> evalf(subs(extpt,HPR[1,1]));

                            -.02488888888

> evalf(subs(extpt,det(HPR)));

                            .002133333334

> extq:=evalf(subs(extpt,f));

                         extq := 15.62500000

> extC:=subs(extpt,1*x+(3/2)*y);

                                     375
                             extC := ---
                                     16

> with(plots):
> implicitplot({f=extq,1*x+(3/2)*y=extC},x=0..30,y=0..30);

# Click on the plot, then 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;

                            1/2    1/2 3/2
                  PR := p (x    + y   )    - x - w y

> gradPR:=grad(PR,[x,y]);

            [        1/2    1/2 1/2              1/2    1/2 1/2    ]
            [    p (x    + y   )             p (x    + y   )       ]
  gradPR := [3/4 ------------------ - 1, 3/4 ------------------ - w]
            [            1/2                         1/2           ]
            [           x                           y              ]

> dem:=solve({gradPR[1],gradPR[2]},{x,y});

                           4        2           4        2
                      81  p  (1 + w)       81  p  (1 + w)
          dem := {y = --- -----------, x = --- -----------}
                      256      4           256      2
                              w                    w

> subs(dem,gradPR[1]);

        /                   / 4        2\1/2
        |        1/2    1/2 |p  (1 + w) |
1/108 p |1/256 81    256    |-----------|
        |                   |     2     |
        \                   \    w      /

                          / 4        2\1/2\1/2
               1/2    1/2 |p  (1 + w) |   |      1/2    1/2   /
     + 1/256 81    256    |-----------|   |    81    256     /
                          |     4     |   |                 /
                          \    w      /   /

    / 4        2\1/2
    |p  (1 + w) |
    |-----------|    - 1
    |     2     |
    \    w      /

> simplify(");

      // 4        2\1/2   / 4        2\1/2\1/2   / 4        2\1/2
      ||p  (1 + w) |      |p  (1 + w) |   |      |p  (1 + w) |
    p ||-----------|    + |-----------|   |    - |-----------|
      ||     2     |      |     4     |   |      |     2     |
      \\    w      /      \    w      /   /      \    w      /
    -------------------------------------------------------------
                          / 4        2\1/2
                          |p  (1 + w) |
                          |-----------|
                          |     2     |
                          \    w      /

# 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 necessary to convert dem[2] 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);

                                   4        2
                              81  p  (1 + w)
                        yd := --- -----------
                              256      4
                                      w

> ydem:=unapply(yd,p,w);

                                         4        2
                                    81  p  (1 + w)
                  ydem := (p, w) -> --- -----------
                                    256      4
                                            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 surfaces. Does the shape of the demand surface conform to your
# expectations? (manipulate the plot to make it more understandable) 
> 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);

       /                   / 4        2\1/2
       |        1/2    1/2 |p  (1 + w) |
sup := |1/256 81    256    |-----------|
       |                   |     2     |
       \                   \    w      /

                          / 4        2\1/2\3/2
               1/2    1/2 |p  (1 + w) |   |
     + 1/256 81    256    |-----------|   |
                          |     4     |   |
                          \    w      /   /

> qsup:=unapply(sup,p,w);

                  /                   / 4        2\1/2
                  |        1/2    1/2 |p  (1 + w) |
qsup := (p, w) -> |1/256 81    256    |-----------|
                  |                   |     2     |
                  \                   \    w      /

                          / 4        2\1/2\3/2
               1/2    1/2 |p  (1 + w) |   |
     + 1/256 81    256    |-----------|   |
                          |     4     |   |
                          \    w      /   /

> qsup(3,1);

                          1/2    1/2    1/2    1/2 3/2
               1/16384 128    (81    256    324   )

> plot(qsup(p,1),p=3..6);

> plot({qsup(p,1),qsup(p,2),qsup(p,3)},p=3..6);

> 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);

         /                   / 4        2\1/2
         |        1/2    1/2 |p  (1 + w) |
PRm := p |1/256 81    256    |-----------|
         |                   |     2     |
         \                   \    w      /

                          / 4        2\1/2\3/2        4        2
               1/2    1/2 |p  (1 + w) |   |      81  p  (1 + w)
     + 1/256 81    256    |-----------|   |    - --- -----------
                          |     4     |   |      256      2
                          \    w      /   /              w

            4        2
       81  p  (1 + w)
     - --- -----------
       256      3
               w

> evalf(subs({p=3,w=1},PRm));

                              68.3437503

> PRmax:=unapply(PRm,p,w);

                     /                   / 4        2\1/2
                     |        1/2    1/2 |p  (1 + w) |
PRmax := (p, w) -> p |1/256 81    256    |-----------|
                     |                   |     2     |
                     \                   \    w      /

                          / 4        2\1/2\3/2        4        2
               1/2    1/2 |p  (1 + w) |   |      81  p  (1 + w)
     + 1/256 81    256    |-----------|   |    - --- -----------
                          |     4     |   |      256      2
                          \    w      /   /              w

            4        2
       81  p  (1 + w)
     - --- -----------
       256      3
               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? 
