21
$\begingroup$

In Lay's Linear Algebra and Its Applications textbook, he defines the matrix $$A=\begin{bmatrix} 4 & 11 & 14\\8 & 7 & -2 \end{bmatrix}$$ and claims that the transformation $T(x)=A x$ maps the sphere to an ellipse.

enter image description here

I'm still struggling with the mathematics on why this is true, but I would like to be able to come up with a way in Mathematica to demonstrate the transformation of the unit sphere to this ellipse to my students.

Update: There is some wonderful support on this page. Thanks to all, but I'd also like to see if it is possible to use Mathematica to produce the equation of the ellipse that is the border of the region.

$\endgroup$
21
$\begingroup$

Another visualization you might find helpful:

sp = MeshCoordinates[
   DiscretizeGraphics[Sphere[], MaxCellMeasure -> 0.02]];
el = Map[(A.#)~Join~{-10} &, sp];

I used DiscretizeGraphics and MeshCoordinates to get the points. From here I apply the mapping, but add a "z" coordinate for display purposes so you can see where each point on the sphere gets mapped to on the ellipse.

Show[Graphics3D[{Opacity[0.5, Blue], Sphere[], Arrowheads[0.02], 
   MapThread[Arrow[{##}] &, {sp, el}] }], ListPlot3D[el], 
 ImageSize -> Large, Axes -> True]

Mathematica graphics

$\endgroup$
20
$\begingroup$
x = Cos[ϕ] Sin[θ];
y = Sin[ϕ] Sin[θ];
z = Cos[θ];
A = {{4, 11, 14}, {8, 7, -2}};
ParametricPlot[
 A.{x, y, z}, {θ, 0, Pi}, {ϕ, 0, 2 Pi}, 
 MeshStyle -> {Red, Blue}, ColorFunction -> "MintColors"]

enter image description here

$\endgroup$
  • 3
    $\begingroup$ Very nice: shows what happens to the latitudes and longitudes, facilitating OP's interpretation of the actual transformation. $\endgroup$ – march Sep 30 '16 at 19:51
  • $\begingroup$ @yarchik Is all the code present here for drawing this image? I am using Mathematica 11 and I'm not getting all of those mesh lines. $\endgroup$ – David Sep 30 '16 at 21:15
  • $\begingroup$ @David, yes, I have MA9 $\endgroup$ – yarchik Sep 30 '16 at 21:22
  • 1
    $\begingroup$ On Mathematica 10 I get the lines only after I add Mesh -> 20 to the ParametricPlot. @David $\endgroup$ – Rahul Sep 30 '16 at 22:01
  • $\begingroup$ I'd like to thank all of my colleagues for all of their wonderful assistance. This page has some absolutely wonderful answers. $\endgroup$ – David Oct 1 '16 at 16:10
13
$\begingroup$

Visualization

From a GeometricTransformation on the unit sphere:

unitsphere = Sphere[];
transformTo2D = {{4, 11, 14}, {8, 7, -2}};
transformIn3D = Append[transformTo2D, {0, 0, 0}];

g = Graphics3D[
  {Black, GeometricTransformation[unitsphere, transformIn3D]}, 
  Boxed -> False
]

enter image description here

We can visualize the resulting ellipse together with the initial sphere with:

Graphics3D[
   {Red, unitsphere, Opacity[0.5], Black, 
    GeometricTransformation[unitsphere, transformIn3D]
   }, Boxed -> False
]

enter image description here

Equation of the ellipse

I'd also like to see if it is possible to use Mathematica to produce the equation of the ellipse that is the border of the region.

Parameters of the ellipse.

With mat the transformation matrix and p a point on the sphere,

mat = {{4, 11, 14}, {8, 7, -2}};
p = {x, y, z};

the length a of the semi-major axis of the ellipse and the coordinates of one of the corresponding points on the boundary are given by:

solmaj = Maximize[{Norm[mat.p], y >= 0}, Element[p, Sphere[]]]
(* {6 Sqrt[10], {x -> 1/3, y -> 2/3, z -> 2/3}} *)

a = solmaj [[1]]
(* 6 Sqrt[10] *)

pmaj = mat.solmaj [[2, All, 2]]
(* {18, 6} *)

Similarly, the length b of the semi-minor axis and the coordinates of one of the corresponding points are given by:

solmin = Maximize[{Norm[mat.p], Projection[mat.p, pmaj] == 0}, Element[p, Sphere[]]]
(* {3 Sqrt[10], {x -> -(2/3), y -> -(1/3), z -> 2/3}} *)

b = solmin[[1]]
(* 3 Sqrt[10] *)

pmin = mat.solmin[[2, All, 2]]
(* {3, -9} *)

The angle between the semi-major axis and the $x$-axis is:

α = VectorAngle[pmaj, {1, 0}]
(* ArcCos[3/Sqrt[10]] *)

Equations.

From the above definitions, we can plot the ellipse with the parametric equation:

ParametricPlot[
   {a Cos[t] Cos[α] - b Sin[t] Sin[α], 
    a Cos[t] Sin[α] + b Sin[t] Cos[α]
   }, {t, 0, 2 π},
   Epilog -> {PointSize[Large], DotDashed,
     Red, Point[{pmaj, -pmaj}], Line[{pmaj, -pmaj}], 
     Brown, Point[{pmin, -pmin}], Line[{pmin, -pmin}]
   }
]

enter image description here

and write the canonical form of the equation

$$ \left(\frac{x \cos{\alpha} + y \sin{\alpha}}{a}\right)^2 + \left(\frac{x \sin{\alpha} - y \cos{\alpha}}{b}\right)^2 = 1 $$

as

((x Cos[α] + y Sin[α])/ a)^2 + ((x Sin[α] - y Cos[α])/b)^2 // Simplify

$$ \frac{13 x^2 - 18 x y + 37 y^2}{3600} = 1 $$

$\endgroup$
  • $\begingroup$ May I ask how you got the canonical form? $\endgroup$ – David Oct 1 '16 at 5:15
  • $\begingroup$ @David See update. $\endgroup$ – user31159 Oct 1 '16 at 11:01
  • $\begingroup$ This will be really helpful. Thanks. $\endgroup$ – David Oct 1 '16 at 15:57
  • $\begingroup$ This is going to be a tremendous learning experience. Great answer! $\endgroup$ – David Oct 1 '16 at 16:09
  • $\begingroup$ @David Glad this was helpful. This was a nice question with many interesting answers! $\endgroup$ – user31159 Oct 1 '16 at 16:33
13
$\begingroup$

Generating equidistributed points on a sphere according to this post:

pts = With[{points = 200, samples = 40000, iterations = 20}, 
   Nest[With[{randoms = Join[#, RandomPoint[Sphere[], samples]]}, 
      Normalize @ Mean @ randoms[[#]] & /@ 
       Values @ PositionIndex @ Nearest[#, randoms]] &, 
    RandomPoint[Sphere[], points], iterations]];

ListPointPlot3D[#, BoxRatios -> 1]& @ pts

enter image description here

A = {{4, 11, 14}, {8, 7, -2}};

ellipse = Dot[A, #]& /@ pts;
(* which is the same as pts.Transpose[A] *)

ell = ListPlot[ellipse, Frame -> True]

enter image description here


The points on a unit sphere above are equidistributed. To make a uniform random distribution of points one needs to generate the spherical angles $b$ and $l$ according to

$$b=\arccos(2u-1)$$ $$l=2\pi v$$

where $u,v$ are $\mathcal{U}[0,1)$ iid:

ang = Table[{ArcCos[2 RandomReal[] - 1], 2 \[Pi] RandomReal[]}, {i, 1000}];

Then transform it to Cartesian coordinates $(x,y,z)$:

sph = {Sin[#1] Cos[#2], Sin[#1] Sin[#2], Cos[#1]}& @@@ ang;
ListPointPlot3D[#, BoxRatios -> 1]& @ sph

enter image description here

(EDIT: As Michael E2 pointed out in the comment, sph can be obtained with a built-in function: sph = RandomPoint[Sphere[], 1000].)

ellipse = Dot[A, #]& /@ sph;
pl2 = RegionPlot[ConvexHullMesh[ellipse], AspectRatio -> 1/GoldenRatio]

enter image description here

Thanks to Bob Hanlon for pointing out the RegionPlot.

$\endgroup$
  • 4
    $\begingroup$ RegionPlot[ConvexHullMesh[ellipse], AspectRatio -> 1/GoldenRatio] $\endgroup$ – Bob Hanlon Sep 30 '16 at 18:21
  • 1
    $\begingroup$ You could also use Show[ConvexHullMesh[ellipse], Frame -> True, AspectRatio -> 1/GoldenRatio] $\endgroup$ – Bob Hanlon Sep 30 '16 at 20:56
  • 2
    $\begingroup$ Note that in V10 we got RandomPoint[], so one may use pts = RandomPoint[Sphere[], 200]. Also note ellipse == pts.Transpose[A]. $\endgroup$ – Michael E2 Oct 1 '16 at 19:14
13
$\begingroup$

Another way:

mat = {{4, 11, 14}, {8, 7, -2}};
sph[{x_, y_, z_}] := x^2 + y^2 + z^2;

ContourPlot[sph[PseudoInverse[mat].{u, v}] == 1,
 {u, -20, 20}, {v, -20, 20},
 Epilog -> {Red, PointSize@Large, Point[{{18, 6}, {3, -9}}], 
   Line[{{{18, 6}, -{18, 6}}, {{3, -9}, -{3, -9}}}]}, 
 PlotLabel -> Simplify[sph[PseudoInverse[mat].{u, v}] == 1]]

Mathematica graphics

(If a set of points $\{p\}$ satisfies $F(p)=0$, then the transformed set $\{q=T(p)\}$ satisfies $F(T^{-1}(q))=0$.)


Aside: You can get the semi-major axes of the transformed ellipse from SingularValueDecomposition.

{U, Σ, V} = SingularValueDecomposition[mat];

U.Σ // MatrixForm

Mathematica graphics

The zero vector could be interpreted as one of the axes of the sphere collapsing under the transformation.

The axes of the sphere represented by the third component v are mapped to the above by the transformation:

mat.V // MatrixForm

Mathematica graphics

Some referencess:

David Austin, "We Recommend a Singular Value Decomposition", Feature Column, AMS (undated?, online)

D. Kalman, "A singularly valuable decomposition: The svd of a matrix", The College Mathematics Journal 27 (1996), no. 1, 2-23; revised 2002


More on the SVD:

If $(x,y)$ denotes the inner product, then $$\|Ax\|^2=(Ax,Ax)=(x,A^TAx) \le \lambda \cdot(x,x) = \lambda\,,$$ where $\lambda$ is the greatest eigenvalue of $A^TA$, which has nonnegative real eigenvalues. Consequently $\sqrt{\lambda}$ is the greatest singular value of $A$ and the maximum of $\|Ax\|$, which is achieved when $x$ is a unit eigenvector of $A^TA$ corresponding to the eigenvalue $\lambda$.


SVD Demo: An ${\bf R}^3 \rightarrow {\bf R}^2$ version of the illustration in the Austin article (ref. above).

Note that in the SVD below, V is an orthogonal matrix and represents a rigid motion that aligns the standard coordinate vectors with the eigenvectors of $A^TA$. In this case it is a reflection and -V is a rotation that accomplishes the same alignment. The demo allows the user to rotate the unit sphere from the pole aligned with the z-axis to a sphere with a pole aligned with the eigenvector corresponding to the singular value $0$. In fact, it may be rotated so that the pole may aligned with any of the three eigenvectors.

{U, Σ, V} = SingularValueDecomposition[mat]; (* from above *)
rot[θ_] = RotationMatrix[θ, First@Pick[##, 1] & @@ Reverse@Eigensystem[-V]];

axesplot[
   xform_,          (* transformation =  mat  or  IdentityMatrix[3]  *)
   angle_           (* rotation angle for  rot[]  *)
   ] := With[{
    axes = Transpose[rot[angle]], 
    colors = ColorData[97, "ColorList"][[{1, 3, 4}]]},
   {Red, PointSize@Large, Thick,
    MapThread[                         (* rotated axes *)
     Function[{axis, color}, {color, Point[axis], Line[{-axis, axis}]}],
     {Transpose[xform.axes], colors}],
    Gray, Thickness[Medium],           (* axes transformed by V *)
    InfiniteLine[{{0, 0, 0}, #}.Transpose[xform]] & /@ Transpose[V]}
   ];
sphParam[θ_, ϕ_] = CoordinateTransform["Spherical" -> "Cartesian", {1, ϕ, θ}];

Manipulate[
 GraphicsRow[{
   (* sphere *)
   Show[
    ParametricPlot3D[rot[t].sphParam[θ, ϕ],
     {θ, 0, 2 Pi}, {ϕ, 0, Pi}, PlotStyle -> None, Mesh -> 15],
    Graphics3D[axesplot[IdentityMatrix[3], -t]],
    AxesLabel -> {"x", "y", "z"}, Ticks -> None, 
    ViewPoint -> Dynamic@vp, ViewVertical -> Dynamic@vv, (* preserves view as t changes *)
    SphericalRegion -> True
    ],
   (* mat.sphere *)
   ParametricPlot[mat.rot[t].sphParam[θ, ϕ],
    {θ, 0, 2 Pi}, {ϕ, 0, Pi}, PlotStyle -> None, 
    Mesh -> 15, BoundaryStyle -> {Gray, Thin}, 
    Epilog -> axesplot[mat, -t], Axes -> None
    ]
   }],
 {t, 0., 2. Pi, TrackingFunction ->  (* sets "stops" at critical angles *)
   (t = Nearest[Pi Range[1, 5, 2]/3, #, {1, 0.07}] /. {{} -> #, {tt_} :> tt}; &),
  Appearance -> "Labeled"},
 (* ViewPoint, ViewVertical variables, without controls *)
 {{vp, {1.3, -2.4, 2}}, None}, {{vv, {0, 0, 1}}, None}
 ]

Mathematica graphics

Aligned with the eigenvectors:

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.