Visualizing linear transformations

Joshua R. Davis, Carleton College, Math 134, Fall 2024

Introduction

Linear transformations are arguably the point of linear algebra. They are heavily used throughout math, science, and technology. We will explore many of their applications throughout the course.

For example, structural geology is the study of how rocks deform through fracturing, folding, etc. Observing, describing, and modeling rock deformation can help us understand basic scientific questions about plate tectonics and the Earth’s history. Studying rock deformation can also help us with practical problems such as assessing earthquake hazards, designing dams, and finding natural resources. Rock deformations are often modeled using linear transformations.

In this Mathematica notebook, we practice visualizing linear transformations. If the geology content helps you understand them, then that’s great. But feel free to ignore the geology, if you just want to focus on the math itself.

Infrastructure

Run the huge cell of code below, but don’t worry about how it works. (Remember that, to run a cell, you click anywhere in the cell, and then press Shift-Enter.)

Clear[transformedCurvesPoints, transformer];
(* Helper function for transformer. *)
transformedCurvesPoints[bound_, curves_, points_, matrix_] :=
	Module[{objects},
		objects = {
			RGBColor[0.0, 0.0, 0.0],
			PointSize[Large],
			Point[Map[Function[point, Chop[point.Transpose[matrix]]], points]],
			AbsoluteThickness[2],
			Line[Map[Function[curve, Chop[curve.Transpose[matrix]]], curves]]};
		If[Length[matrix] == 2,
			Graphics[Join[{
				RGBColor[1.0, 0.0, 0.0],
				Line[{{bound, 0}, {0, 0}}],
				Line[{{-bound, 0}, {-bound/2, 0}}],
				RGBColor[0.0, 1.0, 0.0],
				Line[{{0, bound}, {0, 0}}],
				Line[{{0, -bound}, {0, -bound/2}}]},
				objects]],
			Graphics3D[Join[{
				RGBColor[1.0, 0.0, 0.0],
				Line[{{bound, 0, 0}, {0, 0, 0}}],
				Line[{{-bound, 0, 0}, {-bound/2, 0, 0}}],
				RGBColor[0.0, 1.0, 0.0],
				Line[{{0, bound, 0}, {0, 0, 0}}],
				Line[{{0, -bound, 0}, {0, -bound/2, 0}}],
				RGBColor[0.0, 0.0, 1.0],
				Line[{{0, 0, bound}, {0, 0, 0}}],
				Line[{{0, 0, -bound}, {0, 0, -bound/2}}]},
				objects]]]];
(* Returns a manipulator for a 2D or 3D linear transformation. *)
(* bound determines the length of the coordinate axes and often the view. *)
(* curves and points are some curves and points to be transformed. *)
(* matrix is a t-dependent matrix used to transform the curves and points. *)
transformer[bound_, curves_, points_, matrix_] :=
	Manipulate[transformedCurvesPoints[bound, curves, points, matrix], {t, 0, 1}];
(* Simple 2D and 3D geometry to transform. *)
box2D = {{{1, 1}, {-1, 1}, {-1, -1}, {1, -1}, {1, 1}}, 
{{1, 0}, {0, 0}}, {{0, 1}, {0, 0}}};
circlePoints[{a_, b_}, r_, n_] := 
Table[{r Cos[t] + a, r Sin[t] + b}, {t, 0, 2 Pi, 2 Pi / n}];
circle2D = circlePoints[{0, 0}, 1, 100];
rock2D = Append[box2D, circle2D];
box3D = {{{-1,-1,-1},{-1,-1,1}},{{-1,1,-1},{-1,1,1}},{{1,-1,-1},{1,-1,1}},{{1,1,-1},{1,1,1}},{{-1,-1,-1},{-1,1,-1}},{{-1,-1,1},{-1,1,1}},{{1,-1,-1},{1,1,-1}},{{1,-1,1},{1,1,1}},{{-1,-1,-1},{1,-1,-1}},{{-1,-1,1},{1,-1,1}},{{-1,1,-1},{1,1,-1}},{{-1,1,1},{1,1,1}},
{{0,0,0},{0,0,1}},{{0,1,0},{0,1,1}},{{1,0,0},{1,0,1}},{{0,0,0},{0,1,0}},{{0,0,1},{0,1,1}},{{1,0,0},{1,1,0}},{{0,0,0},{1,0,0}},{{0,0,1},{1,0,1}},{{0,1,0},{1,1,0}}};

2D basic examples

Here is a linear transformation that stretches the plane by a factor of 3 in both directions. Inspect the matrix and play with the animation, to get an idea for what the linear transformation does.

matrix = {{3, 0}, {0, 3}};
MatrixForm[matrix]
transformer[3, rock2D, {}, {{3^t, 0}, {0, 3^t}}]

Inspect the matrix and play with the animation. Then describe this linear transformation in words.

matrix = {{2, 0}, {0, -1}};
MatrixForm[matrix]
transformer[3, rock2D, {}, {{2^t, 0}, {0, 1 - 2 t}}]

This is arguably the most important example! (But not the most interesting.)

matrix = {{1, 0}, {0, 1}};
MatrixForm[matrix]
transformer[3, rock2D, {}, {{1, 0}, {0, 1}}]

What about this one?

matrix = {{-1, 0}, {0, 0}};
MatrixForm[matrix]
transformer[3, rock2D, {}, {{1 - 2 t, 0}, {0, 1 - t}}]

Here is a subtler example. Enter any nonzero vector y. What does the transformation have to do with that y?

y = {1,1};
matrix = Outer[Times, y, y] / Dot[y, y];
MatrixForm[matrix]
transformer[3, rock2D, {}, IdentityMatrix[2] + t( matrix - IdentityMatrix[2])]

2D examples that preserve area and orientation

This “simple shearing” is a rudimentary model for what happens to rock between two tectonic plates that are grinding past each other (such as the Pacific and North America, along the San Andreas fault). Try changing gamma and re-running the cell, to explore the meaning of gamma.

Clear[t];
gamma = 0.5;
vgt = {{0, gamma}, {0, 0}};
MatrixForm[MatrixExp[vgt]]
transformer[3, rock2D, {}, MatrixExp[t vgt]]

This “pure shearing” is a rudimentary model for what happens to rock between two tectonic plates that are ramming into each other (such as India and Asia). Change log(k) to explore its meaning.

Clear[t];
logK = 1.0;
vgt = {{logK, 0}, {0, -logK}};
MatrixForm[MatrixExp[vgt]]
transformer[3, rock2D, {}, MatrixExp[t vgt]]

This is pure rotation. For example, if a rock’s magnetic minerals are not aligned with the Earth’s current magnetic field, then maybe the rock has rotated since it formed (or the Earth’s magnetic field has changed, or both --- geology is complicated). Change theta to explore.

Clear[t];
theta = Pi;
vgt = {{0, -theta}, {theta, 0}};
MatrixForm[MatrixExp[vgt]]
transformer[3, rock2D, {}, MatrixExp[t vgt]]

But don’t take the interpolation seriously

Here are three examples. They all involve the same linear transformation, as you can see by comparing the matrices that are printed out. But they differ in how they interpolate from the start time to the end time.

This example rotates the plane through an angle of pi (meaning, 180 degrees).

Clear[t];
theta = Pi;
vgt = {{0, -theta}, {theta, 0}};
MatrixForm[MatrixExp[vgt]]
transformer[3, rock2D, {}, MatrixExp[t vgt]]

This example rotates the plane through an angle of 3 pi (meaning, 540 degrees). Compared to the previous example, it does an extra rotation of the plane, that is undetectable in the final result. (Change theta to any other odd multiple of pi, to make another example.)

Clear[t];
theta = 3 Pi;
vgt = {{0, -theta}, {theta, 0}};
MatrixForm[MatrixExp[vgt]]
transformer[3, rock2D, {}, MatrixExp[t vgt]]

This example produces the same final result without rotating at all. Instead, it “turns the plane inside out”.

Clear[t];
MatrixForm[{{-1, 0}, {0, -1}}]
transformer[3, rock2D, {}, {{1 - 2 t, 0}, {0, 1 - 2 t}}]

When we talk about a linear transformation of the plane (or any other space), we are talking about the net effect between where points started and where they ended up. There are infinitely many possible “paths” between the starting and ending configuration. Those paths are a topic for math beyond linear algebra.

In other words, I have provided animations, to help you see which points end up where. But you should not take the details of the animations too seriously, because there are infinitely many valid ways to run each animation.

Making your own example

If you’ve tried to edit the code above, to experiment with your own examples, then you might have found it difficult. Here are some instructions that might help.

First, you need to make a 2x2 matrix, whose entries depend on the variable t. The matrix must equal the identity matrix when t = 0. Also, it should equal your desired matrix when t = 1. Here’s an example showing a matrix, what it is when t = 0, and what it is when t = 1.

Clear[t];
mat = {{1, Sin[Pi t]}, {1 - Cos[Pi t], 1}};
MatrixForm[mat]
MatrixForm[mat /. t -> 0]
MatrixForm[mat /. t -> 1]

Then you put that matrix into my transformer function, like this:

transformer[3, rock2D, {}, mat]

To emphasize the point raised earlier, about there being infinitely many possible interpolations, here is the same linear transformation with a simpler interpolation.

Clear[t];
mat = {{1, 0}, {2 t, 1}};
MatrixForm[mat]
MatrixForm[mat /. t -> 0]
MatrixForm[mat /. t -> 1]
transformer[3, rock2D, {}, mat]

3D basic examples

This example is a 3D version of the simple shearing above.

Clear[t];
gamma = 0.5;
vgt = {{0, gamma, 0}, {0, 0, 0}, {0, 0, 0}};
MatrixForm[MatrixExp[vgt]]
transformer[3, box3D, {}, MatrixExp[t vgt]]

This next example is like a 3D version of the pure shearing above.

Clear[t];
vgt = {{-Log[1.5], 0, 0}, {0, -Log[1.5], 0}, {0, 0, Log[2.25]}};
MatrixForm[MatrixExp[vgt]]
transformer[3, box3D, {}, MatrixExp[t vgt]]

This final example is like the previous one, but with extraneous rotation thrown into the interpolation. As you can see by comparing matrices, the two examples are actually showing us the same linear transformation.

Clear[t];
vgt = {{-Log[1.5], -2 Pi k, 0}, {2 Pi k, -Log[1.5], 0}, {0, 0, Log[2.25]}} /. k -> 1;
MatrixForm[MatrixExp[vgt]]
transformer[3, box3D, {}, MatrixExp[t vgt]]

More infrastructure

Now we build up some more complicated infrastructure for more complicated visualizations. Run this giant cell, without worrying about how it works.

Clear[layerCurve, zoneCurves, curvesGraphic, translation, scaling, deformation,isError, passError, euclideanRungeKutta, flow];
cubeCurves = 
	{{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 0, 0}},
	{{0, 0, 1}, {1, 0, 1}, {1, 1, 1}, {0, 1, 1}, {0, 0, 1}},
	{{0, 0, 0}, {0, 0, 1}}, {{1, 0, 0}, {1, 0, 1}}, {{1, 1, 0}, {1, 1, 1}}, {{0, 1, 0}, {0, 1, 1}}};
layerCurve[numSamples_] :=
	Join[
		{{0, 0, 0}},
		Table[{1, s / numSamples, 0}, {s, 0, numSamples}],
		Table[{0, 1 - s / numSamples, 0}, {s, 0, numSamples}]];
zoneCurves[numLayers_, numSamples_] :=
	Module[{layer},
		layer = layerCurve[numSamples];
		Join[
			{{{0, 0, 0}, {0, 0, 1}}, {{1, 0, 0}, {1, 0, 1}}, {{0, 1, 0}, {0, 1, 1}}, {{1, 1, 0}, {1, 1, 1}}},
			Flatten[Table[translation[{layer}, {0, 0, l / numLayers}], {l, 0, numLayers}], 1]]];
curvesGraphic[curves_] :=
	Graphics3D[Join[{
			Thickness[Large],
			Red, Line[{{0, 0, 0}, {0.5, 0, 0}}],
			Green, Line[{{0, 0, 0}, {0, 0.5, 0}}],
			Blue, Line[{{0, 0, 0}, {0, 0, 0.5}}],
			Thickness[Small],
			Black},
		Map[Line, curves]], Boxed->False];
translation[curves_, d_] :=
	Map[Function[point, point + d], curves, {2}];
scaling[curves_, abc_] :=
	Map[Function[point, abc point], curves, {2}];
deformation[curves_, f_] :=
	Map[Function[point, f.point], curves, {2}];
isError[value_] :=
	And[ListQ[value], StringQ[value[[1]]]];
passError[function_, arguments_, header_] :=
	Module[{value}, value = Catch[Quiet[Check[Apply[function, arguments], {"message:", function, arguments}]]];
		If[isError[value],
			Throw[Append[header, value]],
			value]];
euclideanRungeKutta[g_, x_, n_, t0_, t1_] :=
	Module[{h, last, t, i, k1, k2, k3, k4},
		h = (t1 - t0) / n;
		last = g;
		t = t0;
		For[i = 0, i < n, i += 1,
			(* Compute x[last], etc. in a way that lets errors be reported. *)
			k1 = passError[x, {t, last}, {"euclideanRungeKutta: k1: g, n, i:", g, n, i}];
			k2 = passError[x, {t + 0.5 h, last + 0.5 h k1}, {"euclideanRungeKutta: k2: g, n, i:", g, n, i}];
			k3 = passError[x, {t + 0.5 h, last + 0.5 h k2}, {"euclideanRungeKutta: k3: g, n, i:", g, n, i}];
			k4 = passError[x, {t + h, last + h k3}, {"euclideanRungeKutta: k4: g, n, i:", g, n, i}];
			last += (h / 6.0) (k1 + 2.0 k2 + 2.0 k3 + k4);
			t += h];
		last];
flow[curves_, v_, t1_] :=
	Map[Function[point, euclideanRungeKutta[point, v, 10, 0, t1]], curves, {2}];
Clear[transpressionGraphicTF];
transpressionGraphicTF[h0_, gamma_, k_, t_] :=
	Module[{zoneA, zoneB, blockA, blockB, f, disp},
		(* Set up the geometry at time 0. *)
		zoneA = scaling[zoneCurves[4, 4], {2, h0, 2}];
		zoneB = translation[zoneA, {-1, 0, -1}];
		zoneA = translation[zoneA, {-1, -h0, -1}];
		blockA = scaling[cubeCurves, {2, 1, 2}];
		blockB = translation[blockA, {-1, h0, -1}];
		blockA = translation[blockA, {-1, -3h0, -1}];
		(* Compute the block displacement at time t. *)
		f = MatrixExp[t {{0, gamma, 0}, {0, Log[k], 0}, {0, 0, -Log[k]}}];
		disp = f.{0, h0, 0} - {0, h0, 0};
		(* Compute the geometry at time t. *)
		blockA = translation[blockA, -disp];
		blockB = translation[blockB, disp];
		zoneA = deformation[zoneA, MatrixExp[t {{0, gamma, 0}, {0, Log[k], 0}, {0, 0, -Log[k]}}]];
		zoneB = deformation[zoneB, MatrixExp[t {{0, gamma, 0}, {0, Log[k], 0}, {0, 0, -Log[k]}}]];
		curvesGraphic[Join[blockA, blockB, zoneA, zoneB]]];
Clear[velocity, transpressionGraphicRC];
velocity[h0_, gamma_, k_] :=
	Function[{t, x},
		Module[{h, u1, u2},
			h = h0 k^t;
			u1 = gamma h;
			u2 = Log[k] h;
			{u1 / h x[[2]],
			(3 u2) / (2 h) x[[2]] - u2 / (2 h^3) x[[2]]^3,
			(3 u2) / (2 h^3) x[[2]]^2 x[[3]] - (3 u2) / (2 h) x[[3]]}]];
transpressionGraphicRC[h0_, gamma_, k_, t_] :=
	Module[{zoneA, zoneB, blockA, blockB, f, disp},
		(* Set up the geometry at time 0. *)
		zoneA = scaling[zoneCurves[4, 16], {2, h0, 2}];
		zoneB = translation[zoneA, {-1, 0, -1}];
		zoneA = translation[zoneA, {-1, -h0, -1}];
		blockA = scaling[cubeCurves, {2, 1, 2}];
		blockB = translation[blockA, {-1, h0, -1}];
		blockA = translation[blockA, {-1, -3h0, -1}];
		(* Compute the block displacement at time t. *)
		f = MatrixExp[t {{0, gamma, 0}, {0, Log[k], 0}, {0, 0, -Log[k]}}];
		disp = f.{0, h0, 0} - {0, h0, 0};
		(* Compute the geometry at time t. *)
		blockA = translation[blockA, -disp];
		blockB = translation[blockB, disp];
		zoneA = flow[zoneA, velocity[h0, gamma, k], t];
		zoneB = flow[zoneB, velocity[h0, gamma, k], t];
		curvesGraphic[Join[blockA, blockB, zoneA, zoneB]]];

More complicated examples

This linear transformation is a combination of pure shearing and simple shearing. Deformable material is trapped between two strong, rigid blocks. The material is sheared in the x (red) direction and smushed in the y (green) direction. To compensate for the smushing, it extrudes in the z (blue) direction. This is a more accurate model of what is happening to the rock along the San Andreas fault.

Clear[t];
gamma = 1.0;
logK = -0.5;
MatrixForm[MatrixExp[{{0, gamma, 0}, {0, logK, 0}, {0, 0, -logK}}]]
Manipulate[transpressionGraphicTF[0.5, 1.0, Exp[logK], t], {t, 0, 1}]

The preceding model is still very simplistic. For example, you can see that the deforming rock slips in a frictionless way (apparently) against the rigid rock. Depending on the scientific goals of the model, that might be a problem.

The following model is arguably more realistic and definitely more complicated. It is a non-linear model constructed specifically for this kind of geologic situation. You can see that it’s non-linear, because lines bend into curves as the deformation progresses.

Clear[t];
gamma = 1.0;
logK = -0.5;
Manipulate[transpressionGraphicRC[0.5, gamma, Exp[logK], t], {t, 0, 1}]

The preceding non-linear model is still pretty simplistic. In practice, more sophisticated models are built up by gluing together many linear transformations defined on little patches of space. To produce physically reasonable results, the gluing is controlled by mathematical “laws” such as the Navier-Stokes equations. The final result is computed using an algorithm called the finite element method (FEM).

Linear algebra shows up in this work in multiple ways. For example, the deformation on each little patch of space is described by a matrix and hence is definitely linear algebra, right? But how does the FEM compute the numbers that make up those matrices? Well, the numbers are determined by the gluing laws, which turn out to be an enormous system of linear equations! So the FEM solves a giant linear algebra problem, to get numbers, which it then plugs into various matrices, to describe linear transformations.