Students and Mathematica:
Swirl Discoveries. A 2D Linear Differential Equations Investigation

The coefficient matrix of 2-dimensional linear differential equation systems has complex eigenvalues and eigenvectors if and only if there are lines of attraction.  The proof arose from and makes use of a graphical field plot representation of such systems.  The use of Mathematica for both algebraic calculation and generation of these field plots allowed the proof to arise quite naturally.

By Nathan Linger

This paper gives a proof of a result I discovered in a Differential Equations class at the University of Illinois at Urbana-Champaign with Professor Jerry Uhl.  The result is that the coefficient matrix of 2-dimensional linear differential equation systems has complex eigenvalues and eigenvectors if and only if there are lines of attraction.  The proof arose from and makes use of a graphical field plot representation of such systems.  The use of Mathematica for both algebraic calculation and generation of these field plots allowed the proof to arise quite naturally.

Introduction

My differential equations class was studying 2 dimensional linear differential equations of the form

({{x' = a x + b y}, {y' = c x + d y}})

To visualize such a system we plotted the vector field

F[x, y] = {x', y'}

Comparisons between the resulting plots of vector fields reveals properties of the system. Here are two such systems and their corresponding field plots.

x' = 
      2  x - 3 y 
y' = - 4 y

PlotVectorField[(({{2, -3}, {0, -4}})).{x, y},{x,  -1,  1}, {y,  -1,  1},PlotPoints -> 10];

All trajectories in this system are repelled from one of two trajectories onto one of two other trajectories.  We learned that these characteristic trajectories lie on the eigenvectors of the system's coefficient matrix.

x' = 
      x + 3 y 
y' = -2 x + y

PlotVectorField[(({{1, 3}, {-2, 1}})).{x, y},{x,  -1,  1}, {y,  -1,  1},PlotPoints -> 10];

All trajectories in this system end up rotating outwards forever.  We learned that systems of this sort have coefficient matrices with non-real eigenvalues.  These systems were called pure swirlers.

Questions arose about other features of the field plot representations.  In the case of a pure swirler, how can one determine from the coefficient matrix whether its trajectories will swirl clockwise or counter-clockwise.  This is the question that I began investigating.

Definition of Swirl

It is useful to interpret x and y as the coordinates of a particle's position.  Then solutions of the system are possible particle trajectories.  We can already determine a particle's velocity from its position.  Another useful quantity is its acceleration.  This can easily be derived by taking derivatives of the original equations to obtain

({{x'' = a x' + b y'}, {y'' = c x' + d y'}})

So acceleration is obtained by applying the same coefficient matrix to the velocity.

velocity=(({{a, b}, {c, d}})).{x, y}
acceleration=(({{a, b}, {c, d}})).velocity

What we want to define is which way a particle is "turning" at a given point in its trajectory.  Is it  being pulled to the left or right of its current path?  This is the same question as whether or not its acceleration vector is on the left or right of its velocity vector.  The way to determine this is with a cross product.

For a cross product, we need 3 dimensions instead of just two.  So we append a z-coordinate of 0 to each vector.

velocity = Append[velocity, 0]
acceleration = Append[acceleration, 0]

Now we are ready to define Swirl as the z-coordinate of the cross product of the velocity and acceleration of a point at position {x,y}.

CartesianSwirl[(({{a_, b_}, {c_, d_}}))][x_, y_]=Cross[velocity, 
          acceleration][[3]]

It is also useful to think of Swirl in polar coordinates.

PolarSwirl[(({{a_, b_}, {c_, d_}}))][r_, θ_]=CartesianSwirl[(({{a, b}, {c, d}}))][r Cos[θ], 
            r Sin[θ]]//TrigReduce//Simplify

Interpretation of Swirl

Swirl is simply a measurement of which way the imaginary particle at that point is turning.  If it is veering to the left (turning counter-clockwise), Swirl is positive.  If it is veering to the right (turning clockwise), Swirl is negative.  If it is not veering to either the left or the right (it is going in a perfectly straight line), Swirl is zero.

So the system is a pure swirler iff Swirl is never zero.  The direction of pure swirler can be obtained by sampling the sign of swirl at some point.

To drive this point home, let's take a look at some field plots of some systems and their corresponding Swirl measurements (in polar coordinates).

FieldAndSwirlPlot[(({{1, 3}, {-2, 1}}))];

Here we have a clockwise pure swirler.  The polar swirl plot show this: Swirl is negative (clockwise) for all angles.

FieldAndSwirlPlot[(({{-1, -2}, {3, -2}}))];

Here we have a counter-clockwise pure swirler.  The polar swirl plot show this: Swirl is positive (counter-clockwise) for all angles.

FieldAndSwirlPlot[(({{2, 2}, {2, -1}}))];

Here we have a "partitioned" swirler which turns counter-clockwise at some places and clockwise at others.  The polar swirl plot show this: Swirl is positive (counter-clockwise) for some ranges of angles and negative (clockwise) for others.

A Connection to Eigenvectors

Let's investigate the points where Swirl is zero.  Here are the constraints on x and y that make Swirl zero:

Solve[CartesianSwirl[(({{a, b}, {c, d}}))][x, y]==0,{x, y}]

This is a familiar looking result.

Eigenvectors[(({{a, b}, {c, d}}))]

Let's set y=1 in the original solution

Solve[CartesianSwirl[(({{a, b}, {c, d}}))][x, 1]==0,{x}]

The {x,y} coordinates where swirl is 0 are two lines through the origin having the same slope as their corresponding eigenvectors.  So these eigenvector lines partition the plane into sections of equivalent direction (sign) of Swirl.

This is why all trajectories coming from systems with non-real eigenvectors swirl in the same direction (clockwise or counterclockwise) at all points.  If they didn't, i.e. there is one trajectory turning clockwise and another turning counterclockwise, then swirl changes from negative to positive.  Since Swirl is a continuous function it must be zero at some point.  This is impossible since that "point" is lies on a non-real eigenvector.

A Closed Form for Direction of Pure Swirlers

How can we decide which direction it will swirl (counter-clockwise or clockwise) using only the coefficient matrix?  This problem is equivalent to finding the sign of swirl given only the coefficient matrix (if we know it will have the same sign at all points).

A good starting point is to find the extrema of PolarSwirl (when r is fixed and t varies).  So solve for θ after setting the derivative of PolarSwirl (with respect to t) to 0.

D[PolarSwirl[(({{a, b}, {c, d}}))][r, θ],θ]==0

This can be simplified to the equivalent expression

((-a + d)/(b + c)) = Tan[2 θ]

Tsols = 
      Solve[((-a + d)/(b + c)) == 
          Tan[2 θ], θ]

Since Tan[2  θ] has a period of (π/2), any addition of (π/2) is also a solution.

Tsols  =  
      Append[Tsols, {θ -> ((π/2) + θ /. 
                  Tsols)[[1]]}]

Plugging these two values back into the function gives us the extrema of PolarSwirl.  Since the signs of coefficients are unknown, the signs of the extrema are also unknown.

MaxAndMin = (FullSimplify //@ (PolarSwirl[(({{a, b}, {c, d}}))][r, θ]/.Tsols))/. \(\@\(1 + a_\^2\/b_\^2\) -> \@\(b\^2 + a\^2\)\/b\)

These two extrema values are of the form X ± Y where

X=(1/2)(r^2)(c - b)Det[(({{a, b}, {c, d}}))];

and

Y=(1/2)(r^2)Det[(({{a, b}, {c, d}}))](Sqrt[(((b + c))^2) + (((a - d))^2)]);

Expand /@ MaxAndMin === Expand /@ {X + Y, X - Y}

Now we have what we were looking for.  The sign of Swirl in a pure swirler is the sign of X±Y, or just the sign of X.  X is easiest to compute, so we'll use it.  The (1/2) and (r^2) terms don't affect the sign at all, so we are left with the condition (c - b)Det[({{a, b}, {c, d}})] >0 for counter-clockwise swirl and (c - b)Det[({{a, b}, {c, d}})] <0 for clockwise swirl (in the case of a pure swirler).

Evidence of Consistency

Toward the end of the previous section we found that the extrema of PolarSwirl are given by:

\(\(MaxAndMin = 
	{\(-\(1\/2\)\)\ \((\(-b\) + 
                c + \@\(\((b + c)\)\^2 + \((a - d)\)\^2\))\)\ \((b\ c - 
                a\ d)\)\ r\^2, 
          1\/2\ \((b - 
                c + \@\(\((b + c)\)\^2 + \((a - d)\)\^2\))\)\ \((b\ c - 
                a\ d)\)\ r\^2};\)\)

Now all "pure swirlers" will have PolarSwirl extrema of the same sign, since Swirl never equals zero in "pure swirler" systems.  We can test whether the extrema have the same sign by testing whether or not their product  is positive.

\(SameSignCondition\  = 

      Factor[Expand[
		Times\  @@ \ MaxAndMin
		]]\  > \ 0\)

We can divide out all positive values to obtain

(a^2) + 4 b c - 2 a d + (d^2) < 0

This is the exact same condition for the eigenvectors of the coefficient matrix to be complex.  A nice check of consistency.

Eigenvalues[(({{a, b}, {c, d}}))]

References

Davis, B. and Jerry Uhl, J. .Differential Equations and Mathematica. MathEverywhere. Columbus, Ohio 1998

About the author

Nathan Linger
Computer Science / Mathematics Undergraduate
University of Illinois at Urbana-Champaign
rlinger@uiuc.edu

Electronic Subscriptions

Included in the distribution for each electronic subscription is the file swirl.nb containing Mathematica code for the material described in this article.