#      MATH 2270
# PROJECT 5
# Conics and Quadrics
# December 1998
# 
# This project is at the web address
# http://www.math.utah.edu/~korevaar/2270proj5.txt  Save it to your
# directory and open it from Maple as a Maple Text document.
# 
# Conic sections:
> with(linalg):with(plots):#for computations and pictures
> with(student):#to do algebra computations like completing the square 
Warning, new definition for norm
Warning, new definition for trace
# 
# Let's do the problem we did in class on Wednesday, using Maple:  This
# was problem #27 on page 501, which you should look at now.  
#      First let's see what conic section we have, by finding the
# eigenvalues of the quadratic form part of the equation:
> A:=matrix([[5,6],[6,0]]);

                                 [5    6]
                            A := [      ]
                                 [6    0]

> eigenvalues(A); #compute eigenvalues

                                -4, 9

# Since the eigenvalues are of opposite sign the conic must be a
# hyperbola (possibly degenerate).  We can verify this immediately by
# using the implicitplot command: 
# 
> implicitplot(5*x^2 + 12*x*y -12*sqrt(13)*x = 36,
>              x=-10..10,y=-10..10, color=`black`);

#    Yes, it was a hyperbola.
#    Now, what about all that work we did to explicitly pick new
# coordinates in which there was no cross term?  If we really need to do
# that, we can try to let Maple help, and although the commands seem a
# little cumbersome at first, when we're done we'll have a template that
# will work for any quadratic equation with only minor modifictations.  
#      Let's write the quadratic equation in standard form.  First we
# need the matrix B for the linear term. 
> B:=matrix([[-12*sqrt(13),0]]);
> 

                       B := [-12 sqrt(13)    0]

> fmat:=v->evalm(transpose(v)&*A&*v
>                  + B&*v ) ;
>      #this function takes a vector v and computes 
>      # transpose(v)Av + Bv ,
>      #which will be a one by one matrix
> f:=v->fmat(v)[1]; #this extracts the entry of the one
>      #by one matrix fmat(v)
>            

   fmat := v -> evalm(((transpose(v) `&*` A) `&*` v) + (B `&*` v))


                         f := v -> fmat(v)[1]

> fmat(vector([x,y])); #should be a one by one matrix
> f(vector([x,y]));  #should get its entry, which is what we want 
> 

               [-12 sqrt(13) x + (5 x + 6 y) x + 6 x y]


                -12 sqrt(13) x + (5 x + 6 y) x + 6 x y

> eqtn:=f([x,y])=36;

         eqtn := -12 sqrt(13) x + (5 x + 6 y) x + 6 x y = 36

> implicitplot(eqtn,x=-10..10,y=-10..10, color=`black`);
>     #this should give you the same picture as above
> 

# Now let's go about the change of variables:
> data:=eigenvectors(A);#get eigenvectors

           data := [-4, 1, {[1, -3/2]}], [9, 1, {[3/2, 1]}]

# You pick things out of the object above systematically, using inidices
# to work through the nesting of brackets:
> data[1];#first piece of data
> data[1][1];#eigenvalue
> data[1][2];#algebraic multiplicity
> data[1][3];#basis for eigenspace
> data[1][3][1];#actual eigenvector
# So that's how to extract the eigenvectors:
> v1:=data[1][3][1];#first eigenvector
> v2:=data[2][3][1];#second eigenvector
> u1:=v1/norm(v1,2);#normalized
> u2:=v2/norm(v2,2);#normalized

> P:=augment(u1,u2);#our orthogonal matrix

> f(evalm(P&*[u,v]));#do the change of variables

> simplify(%);#simplify it!

> completesquare(%,u); #complete the square in u

> completesquare(%,v);
> Eqtn:=%=36;

> standardform:=Eqtn/rhs(Eqtn);#rhs picks off the
>      #right hand side of an equation 

# That's what we got in class.
# 
Let's collect everything into one template .  This time it has  the
data from problem 30 on page 501, but obviously you could plug data
from any quadratic expression into it.   You can see from where the
colons and semicolons are that the output will be  the eigenvalues, 
the transition matrix,  a plot, and an equation just short of standard
form.  (You will get an amusing plot at first but can improve it by
picking a denser grid size.)  
> restart;with(linalg):with(plots):with(student):
Warning, new definition for norm
Warning, new definition for trace
> A:=matrix([[8,-8],[-8,8]]):
> B:=matrix([[33*sqrt(2),-31*sqrt(2)]]):
> C:=70:#the constant term if the rhs is zero
> fmat:=v->evalm(transpose(v)&*A&*v
>                  + B&*v + C):
> f:=v->fmat(v)[1]:
> eigenvals(A); #show the eigenvalues
> data:=eigenvectors(A):
> v1:=data[1][3][1]:#first eigenvector
> v2:=data[2][3][1]:#second eigenvector
> u1:=v1/norm(v1,2):#normalized
> u2:=v2/norm(v2,2):#normalized
> P:=augment(u1,u2);#show our orthogonal matrix
> implicitplot(f([x,y])=0,x=-10..10,y=-10..10,
>      grid=[25,25],color=`black`);
>      #increase grid size for better picture, but too big
>      #takes too long 
> f(evalm(P&*[u,v])):#do the change of variables
> simplify(%):#simplify it!
> completesquare(%,u): #complete the square in u
> completesquare(%,v):#and v
> Eqtn:=%=0;
> 
# 1)  Check your answers to the book homework problems #25,26, by using
# the template above, with the correct data from those problems.
# 
# Quadric Surfaces:
# 2)  Partially do your book homework from section 8.11 by having Maple
# find the eigenvalues of the quadratic part of the equations so that
# you can compute the inertia and classify the surface.  Also verify
# that the plot is what you claimed by drawing a picture of it with
# Maple.  By the way, your book HW problems are 8.10,  hand in
# #10,14,15,23,25,26; recommended 1,4
#      Hints:  To make 3-d plots you can use implicitplot3d.  If you
# make your grid too fine the plot will take a long time to make, so try
# the default grid at first and adjust if necessary.  You might also
# have to adjust your limits in the plot to get a better picture.  You
# can manipulate 3d plots with your mouse.  There are a lot of
# interesting plot options which you write into the command or access
# from the plotting toolbar.  One that helps me see things is to use the
# ``boxed'' axes option.  
#      The 3d plotting routines have been known to crash Maple, so save
# your file often.  Here's a plot to play with, #1 on page 511:
> implicitplot3d(x^2 + y^2 +2*z^2 - 2*x*y -4*x*z -4*y*z + 4*x = 8,
>               x=-5..5,y=-5..5,z=-5..5, axes=`boxed`);

# 3) (optional) Extra credit/for fun?:  Create a template to do the
# reduction into standard form for quadric surfaces, like the one
# written above for conic sections, and test it on several problems from
# the book.  By the way, if you do this you might want to keep in mind
# that if your eigenspace is more that one-dimensional Maple might not
# return an orthogonal basis for it, so you might have to Gram-Shmidt
# your basis in that case.  And in that case the command ``eigenvector''
# will return a different-looking object too.  Of course, ``most''
# quadratic forms don't have that ``problem'' and it is O.K. for this
# project if you don't worry about it.  Something you could add to your
# template is a little subroutine which outputs the type of quadric
# surface, based on an analysis of the eigenvalues.