# Math 2270
#  Maple Project 2 
# September, 1997
# 
#      We will study the geometric meaning of linear (hence matrix) maps
# f(x)=Ax, from R^n to R^m.  As does the book, we will focus on maps
# from R^2 to R^2 to illustrate more general properties.  Of course,
# computer graphics are most concerned with those maps and with R^3
# rotation maps and projection maps from R^3 to R^2.  
#      A Maple text version of this project may be found at the web site
# http://www.math.utah.edu/~korevaar/2270proj2.txt.     I recommend
# saving the file from the website onto you home directory, and then
# opening it from Maple, as a ``Maple text'' document.  You can then do
# the assignment by inserting the proper commands as well as italicized
# textual comments in answer to the various questions.  In order to
# execute multiline commands from the version you get from the web, you
# should use the ``Edit'' option in your maple window to ``join''
# execution groups which you have highlighted with your mouse.  
#     The very definition of a linear map, namely that f(u+v)=f(u)+f(v)
# and f(su)=sf(u), for all vectors u,v and scalars s, has the following
# geometric onsequences:
# (1)  Lines are mapped to lines, and parallel lines are mapped to
# parallel lines:  This is because a line can be described as a set L={u
# + tv| t in R}, where u is a point on the line and v is a direction
# vector.  Therefore the image set f(L):={f(u+tv)}={f(u)+tf(v)} is also
# a line, going through f(u) with direction f(v).  (If the directions
# f(v) turns out to be 0, then the line degenerates into a point.) 
# Therefore, line segments are mapped to line segments,  polygons are
# mapped to polygons, and regions bounded by polygons are mapped to
# regions bounded by polygons; if we know where the vertices go, we know
# everything.  Let's use these facts to draw the images of some
# polygonal regions under a matrix map.
#      First load the linear algebra and plotting libraries
> with(plots):with(linalg):
> verts0:=[[0,0],[1,0],[1,1],[0,1]]; 
>       #corners of unit square
> unitsquare:=polygonplot(verts0, color=`yellow`):
>       #this command make a polygon and colors
>       #the region inside it yellow.  Make sure
>       #to end this command with a colon!
> display(unitsquare);  #Now semicolon! this command
>       #shows the square
# When you executed the sequence of commands above you should have
# gotten a picture of a yellow unit square.  Now we will use a linear
# map with matrix A defined below to map this square to a parallelegram
> A:=matrix([[3,2],[-1,1]]);  #the matrix of our
>        #random linear transformation

> f:=x->evalm(A&*x);  #our linear map
# 1a) For our map f(x)=Ax defined above, what are the images of the
# vectors e1=[1,0] and e2=[0,1]?  How do you find these images from the
# rows or columns of A?
We use the ``map'' command below to see where the vertices of the unit
square are sent by f.  The syntax is to put the mapping function in as
the first argument, and the list of input vectors as the second
argument. The result of the command  will be the list of output
vectors.  Check (not to hand in) that this is what happens below to
the four points in the list verts0.   
> verts1:=map(f,verts0);
> image1:=polygonplot(verts1, color=`red`):
> display({unitsquare,image1});
# 1b)  Explain where f(e1) and f(e2) are represented in the picture you
# just made.  By the way, if you click on the plot you can choose to
# have the projection ``constrained '', in which case Maple will use the
# same scales for the vertical and horizontal axis.  Otherwise it scales
# the x and y-axes to make the picture fit nicely onto your screen.
# 
# 1c) If we compute f(f(x)):=f^2(x), then what will the matrix be for
# the resulting linear transformation? After answering that question,
# make the vertices for f(f(unitsquare)) as follows, and draw the
# corresponding image polygons.
> verts2:=map(f,verts1);
> image2:=polygonplot(verts2,`color`=blue):
> display({unitsquare,image1,image2});

# 1d)  Explain what the columns of the matrix for f^2 have to do with
# the picture you just made above.
# 1e)  What is the inverse function to f?  Hint : what is its matrix? 
# Verify that the inverse mapping takes image1 back to the unit square.

# 2) Translations of objects are mapped to translations of the image
# objects and scalings of objects are mapped to scalings of the image
# objects:  
#      First let's review (or explain) what we mean by an ``object'' and
# by translations and scalings of an object.  An object  is some set S
# of points s.   By the  image of S we mean the collection of image
# points f(s) .  We write f(S) for this image (like we did for the line
# L and its image f(L) iin problem 1).   Similarly if b is a translation
# vector, then the object S+b means all points of the form s+b where s
# is in S.  If c is a scalar, then the scaled (or dilated) set cS means
# all points of the form cs, where s is in S.
#      Now, if we apply the linear map f to the translated object S+b we
# get all points of the form f(s+b)=f(s)+b, where s is in S (since f is
# linear), i.e. exactly the set f(S)+f(b), which is the translation of
# the image f(S) by the vector f(b).  Similarly, if we apply f to the
# scaled set cS be get f(cS) to be the set of all points f(cs)=cf(s)
# (since f is linear), i.e. the scaling cf(S) of the image set f(S).
#       Here's how to translate and scale the unit square from #1: 
# Let's translate it by the vector b=[2,3]:
> b:=[2,3];
> trans:=x->evalm(x+b);
> verts3:=map(trans,verts0);
> image3:=polygonplot(verts3,color=`yellow`):
>       #colon!
> display(unitsquare,image3);
# Here's how to scale it by a factor of 0.2:
> c:=0.2;
> shrink:=x->evalm(c*x);
> verts4:=map(shrink,verts0);
> image4:=polygonplot(verts4,color=`red`):
> display({image4,unitsquare});
# (By the way, when I do the command above the little square gets hidden
# behind the big one, except for its outline.  One way to see all of it
# is to uncolor the big one, but you don't have to do this.)
# 2a)  Describe where you expect the translated unit square, image3
# above, to be  mapped  by our linear map f from problem 1.  Then use
# the ``map'' command and the ``trans'' command, as well as polygonplot
# and display, to make a picture of the images of the unit square and
# its translation when they are mapped by f .  Verify that the images
# differ by the expected translation.
# 2b)  Describe what you expect the image of the scaled down square
# above to be when you apply f , and then draw a picture illustrating
# this.
# 
# 
# 
# 3)  Special linear transformations in R^2:
# 3a) Rotations:  The matrix for rotating by an angle theta is:
> Rot:=theta->matrix([[cos(theta),-sin(theta)],[sin(theta),cos(theta)]])
> ;

                             [cos(theta)    -sin(theta)]
             Rot := theta -> [                         ]
                             [sin(theta)    cos(theta) ]

# so to rotate the vector (2,3) by Pi/3, we would command
> theta:=Pi/3;
> evalm(Rot(theta)&*[2,3]);
> 
# You are to draw a picture of the unit square rotated by Pi/3 radians.
# 3b)  Reflections:  Define matrices which give reflections across the
# x-axis, across the y-axis, and across the line y=x.  Illustrate what
# each of these linear maps does to the unit square.  Finally, verify
# with Maple that  reflecting across the line y=x is the composition of
# first rotating Pi/4 in the clockwise direction, then reflecting across
# the x-axis, then rotating (back ) Pi/4 in the counterclockwise
# direction.
# 3c)  Translations and scalings revisited:  Explain why a translation
# by a non-zero map is NOT LINEAR.  (Show that the definition of linear
# does not hold.)  On the other hand, explain why  scaling IS a linear
# map and exhibit its matrix.  (This is a book homework problem
# disguised in a computer project.)
# 3d)  You can scale by different factors in the x and y-directions. 
# For example, what matrix expands the x-direction by 3 and the
# y-direction by 2?  Make a picture of what this scaling does to the
# unit square.
# 3e)  Projections:  Write down the matrix which projects points onto
# the x-axis.  Illustrate what this projection map does to the unit
# square.  Is there an inverse map? 
# 3f) Shears:  A shear of strength k in the x-direction is defined by
# the matrix
> Shear:=k->matrix([[1,k],[0,1]]);

                                      [1    k]
                        Shear := k -> [      ]
                                      [0    1]

# For k=1 explore what the shear map does to the unit cube.  What
# happens if you apply the shear map twice?  Three times?  What is the
# inverse map of a strength k shear?
# 3g)  Movies:  As you discovered in (3f), composing a shear with itself
# is another shear of twice the strength.  Make the following movie;
# it's a tiny indication of how computers and sequences of
# transformations can be used to make moving pictures. I stole the code
# from chapter 5 of the text ``Multivariable Mathematics with Maple'',
# by Carlson and Johnson.  All you have to do is execute it.  You
# should, however, pause to think about what the commands are doing.  
> iter:= proc(g,s,n)   #take a set of points s and
>      # apply the map g to it n times in a row.
>      # Save the (n+1) sets which you generate in
>      # an array called d.  Then convert d into a
>      # list.   
> local d,i;   #local variables won't pollute your use
>      #of the same letters elsewhere 
> d:=array(0..n);
> d[0]:=s;
> for i from 1 to n do
>   d[i]:=map(g,d[i-1]);  #apply the map g to the
>      #set you have at stage i-1 to get the set
>      #at stage i.
> od;
> convert(d,list);  #because you use a ``list'' data
>      #structure in the next routine.
> end:
> film:=proc(d)   #make a sequence of polygonplots
>      #from a sequence of vertex sets like the
>      #one you can generate above using iter.  In the
>      #end you can view your film with the display
>      #command.
> local i,j,n,F;
> n:=nops(d);
> F:=array(1..n);
> for i from 1 to n do
>    F[i]:=polygonplot(d[i]);
> od;
> convert(F,list);
> end:
# Now we use the two subroutines above to make a movie of a square
# getting sheared.
> k:=0.1;     #make a film of iterated shears
>      #where the shear strength of one step is 0.1
> g:=x->evalm(Shear(k)&*x);  #the mapping function

                               k := .1


                   g := x -> evalm(Shear(k) `&*` x)

> steps:=iter(g,verts0,20);  #make a list of 20 sets
>      #of vertices, starting with the unit cube vertex
>      #set and applying the shear 20 times in succession. 

> display(film(steps), insequence=true);
>      #a movie of a square getting squashed sideways.
>      #Use your menu options to view the movie.
>      #The command insequence=true means you want
>      #to see the frames in sequence, rather than
>      #a display of all the frames at once.

#