# MAPLE ASSIGNMENT 7
# EIGENVALUES & EIGENVECTORS
# Maple has commands for computing eigenvalues and eigenvectors as
# illustrated below. If A is not symmetric some eigenvalues may be
# complex numbers. You might like to edit A and try this out.
> A:=matrix(4,4,[1,-2,3,-4,-2,3,-4,5,3,-4,5,-6,-4,5,-6,7]);
> det(A);
> with(linalg):
> lambda:=evalf(Eigenvals(A));
# Note that two of these numbers are 0 and what you are seeing are
# rounding errors.
> evalf(Eigenvals(A,vecs));
> X:=evalm(vecs);
# Now to verify that Maple is getting the right answer
> evalm(A&*X);
> evalm(lambda[1]*col(X,1));
# Note that again there are slight rounding errors. Edit this command to
# verify that the other columns of X are also eigenvectors of A. Maple
# will also produce the characteristic polynomial of A and this clearly
# shows that 2 of the eigenvalues are zero.
> charpoly(A,x);
# DEFINITENESS
# Maple has commands for checking the definiteness of a symmetric matrix
# without having to compute the signs of the principal minors.
> definite(A, 'positive_def');
> definite(A, 'negative_def');
> definite(A, 'positive_semidef');
> definite(A, 'negative_semidef');
# You already knew these answers. You can edit the commands above to
# show that transpose(A&*A is pos. semidef. but not pos. def..
> B:=evalm(transpose(A)&*A);
# The following calculations illustrate that A is not pos or neg semi
# def by computing several values of transpose(X)&*A&*X. 
> X1:=vector([-2,-3,0,1]);
> evalm(X1&*A&*X1);
> X2:=vector([-8,3,1,-6]);
> evalm(X2&*A&*X2);
> X3:=vector([0,0,0,1]);
> evalm(X3&*A&*X3);
# GRAPHS OF 2 VARIABLE FUNCTIONS
# The first example of a 3d plot has nothing to do with economics but
# may be attractive.
> plot3d(sin(x)*cos(y),x=-Pi..Pi,y=-Pi..Pi);
# If you click on the plot, a number of menus will appear along the top
# of the window
# which allow you to change the appearance of the plot. The AXES menu
# allows you to put in axes in
# several different ways.  The STYLE menu allows you to change the way
# the plot is shaded.  Try
# some of these options out. When you click on a menu item, your plot
# disappears into a box.  Click on the REDRAW button to see the changed
# plot. (marked R). If you wish to print your plot click the PRINT
# command on the FILE menu.
> with(plots);
# 
# Now try plotting the production function I will be using as an
# example in class. A good choice of 
# style here may be "patch and contour". Putting in normal axes may also
# help
> plot3d(x^.25+x*y+y^.25,x=0..50,y=0..50);

# I did not like the direction of the axes in this plot.  The directions
# can be controlled by changing the orientation angles given at the left
# of the second toolbar when you click on the plot. In the plot you just
# drew you are looking at the origin at an
# angle of 45 degrees from the positive x-axis and an angle of 45
# degrees above the (x,y) plane. 
# These view angles can be changed by editing the numbers in the
# toolbar. For example if the first 45 is changed to -45 then  the x and
# y axes are rotated 90 degrees counterclockwise.  You 
# can experiment with other choices if you wish. 
> plot3d(x^.25+x*y+y^.25,x=0..50,y=0..50,orientation=[-45,45]);
# The domain over which the surface is plotted can be changed from a
# rectangle by making the limits for y functions of x. Here the domain
# is reduced
# to{(x,y)/0<=x,0<=y,x+y<=50}. I recommend changing the first
# orientation angle to -15 degrees to get a better looking plot as well
# as making other changes to STYLE and AXES to suit your own tastes.. 
> plot3d(x^.25+x*y+y^.25,x=0..50,y=0..50-x);

# PARTIAL DERIVATIVES
> f:=(x,y)->sin(x)*cos(y);

                     f := (x, y) -> sin(x) cos(y)

> g:=(x,y)->x^.25+x*y+y^.25;

                                   .25          .25
                   g := (x, y) -> x    + x y + y

# The partial derivatives of multiple variable functions can be
# obtaiined using the diff command.
> fx:=diff(f(x,y),x);

                         fx := cos(x) cos(y)

> fy:=diff(f(x,y),y);

                         fy := -sin(x) sin(y)

> gx:=diff(g(x,y),x);

                                  .25
                            gx := ---- + y
                                   .75
                                  x

> gy:=diff(g(x,y),y);

                                      .25
                            gy := x + ----
                                       .75
                                      y

# Second partial derivatives can also be computed.
> gxx:=diff(g(x,y),x,x);

                                     .1875
                            gxx := - -----
                                      1.75
                                     x

> gxy:=diff(g(x,y),x,y);

                               gxy := 1

> gyx:=diff(g(x,y),y,x);

                               gyx := 1

> gyy:=diff(g(x,y),y,y);

                                     .1875
                            gyy := - -----
                                      1.75
                                     y

# Note that gxx and gyy are always negative. This is the law of
# diminishing marginal products. Notice also that gxy = gyx.
# Now you can plot the graphs of the marginal product functions. Change
# the appearance of the plot in any way that like.  Can you see the law
# of diminishing marginal products in these graphs?
> plot3d(gx,x=0..1,y=0..1);

# Diminishing MP's are indicated if when you point the cursor at the
# y-axes and then move the cursor 
# from there parallel to the x-axis you cross ever lower and lower
# contour lines. Notice that I had to 
# reduce the domain for the plot substantially in order to look at this
# property as  gx(x,y) quickly 
# becomes a constant = y as a function of x alone. The next command
# plots gy.  Does this plot also 
# show DMP?
> plot3d(gy,x=0..1,y=0..1);

# You might like to change g(x,y) at this stage to something more well
# known among economists such 
# as a Cobb-Douglas g(x,y)=x^alpha*y^beta where you can make your own
# choices for alpha & 
# beta. Both should be < 1 and probably the sum <=1 is a good choice
# too. You can probably enlarge 
# the domains for the plots of the MPs and still be able to see DMP.
# 
# Partial derivatives of the first and second orders can also be
# calculated all at once.
> with(linalg):
> grad(f(x,y),[x,y]);

                   [cos(x) cos(y), -sin(x) sin(y)]

> grad(g(x,y),[x,y]);

                         [.25           .25 ]
                         [---- + y, x + ----]
                         [ .75           .75]
                         [x             y   ]

> hessian(f(x,y),[x,y]);

                  [-sin(x) cos(y)    -cos(x) sin(y)]
                  [                                ]
                  [-cos(x) sin(y)    -sin(x) cos(y)]

> hessian(g(x,y),[x,y]);

                         [  .1875           ]
                         [- -----       1   ]
                         [   1.75           ]
                         [  x               ]
                         [                  ]
                         [             .1875]
                         [   1       - -----]
                         [              1.75]
                         [             y    ]

# In these commands f(x,y) could be replaced by the expression f(x,y) is
# equal to.  The variables do 
# not have to be named x & y. It is possible to do a subset of the
# variables the function is defined over.
> grad(K^(1/2)*L^(1/3)*R^(1/6),[K,L,R]);

            [     1/3  1/6       1/2  1/6       1/2  1/3]
            [    L    R         K    R         K    L   ]
            [1/2 ---------, 1/3 ---------, 1/6 ---------]
            [       1/2            2/3            5/6   ]
            [      K              L              R      ]

> grad(K^(1/2)*L^(1/3)*R^(1/6),[K,R]);

                    [     1/3  1/6       1/2  1/3]
                    [    L    R         K    L   ]
                    [1/2 ---------, 1/6 ---------]
                    [       1/2            5/6   ]
                    [      K              R      ]

> hessian(K^(1/2)*L^(1/3)*R^(1/6),[K,L,R]);

       [       1/3  1/6            1/6                 1/3    ]
       [      L    R              R                   L       ]
       [- 1/4 ---------     1/6 ---------      1/12 --------- ]
       [         3/2             1/2  2/3            1/2  5/6 ]
       [        K               K    L              K    R    ]
       [                                                      ]
       [        1/6               1/2  1/6             1/2    ]
       [       R                 K    R               K       ]
       [ 1/6 ---------     - 2/9 ---------     1/18 --------- ]
       [      1/2  2/3              5/3              2/3  5/6 ]
       [     K    L                L                L    R    ]
       [                                                      ]
       [        1/3                1/2                1/2  1/3]
       [       L                  K                  K    L   ]
       [1/12 ---------     1/18 ---------     - 5/36 ---------]
       [      1/2  5/6           2/3  5/6               11/6  ]
       [     K    R             L    R                 R      ]

> hessian(K^(1/2)*L^(1/3)*R^(1/6),[K,R]);

                [       1/3  1/6             1/3    ]
                [      L    R               L       ]
                [- 1/4 ---------     1/12 --------- ]
                [         3/2              1/2  5/6 ]
                [        K                K    R    ]
                [                                   ]
                [        1/3                1/2  1/3]
                [       L                  K    L   ]
                [1/12 ---------     - 5/36 ---------]
                [      1/2  5/6               11/6  ]
                [     K    R                 R      ]

