Wednesday, October 21, 2009

Stability analysis

I haven't been posting for a few months because I've been suffering from another bout of depression (as I have many times before) which leaves me mostly bed-ridden and unable to concentrate. Thanks to Welbutrin and an SSRI I'm starting to get better although I still tire easily and spend many days partially in bed. That said, at least I can now concentrate long enough to get a little work done and use my brain for a few hours a day. This one blog entry has taken me about two weeks to write. Thanks to everyone who's been helpful and understanding while I've been sick.

To help with the amorphous computing projects, I've been trying to get my head around concepts of stability analysis in differential equations. In most of the textbooks I've read, stability is discussed with an analogy to balls rolling on curved surfaces such as the one below showing the ball (representing the state of some system) near an unstable equilibrium at the the top of a hill and thus tending to roll down into a stable equilibrium at the bottom (these analogies ignore the momentum -- its as if the ball, i.e. the system, is moving through viscous honey).

While these "ball rolling" analogies are helpful, they become detrimental to understanding when the functions have complex-valued curvature. This confused me for a long time until I realized that I should think about systems as fields using quiver plots like the following.

In a quiver plot, an arrow is drawn to indicate the solution to the function at each point and this permits viewing both the magnitude and direction of the field. In the above case, we see a 2D section of a counter-clockwise rotational field from the Lotka Volterra equation which is a simplistic model of a predator / prey relationship; the details of this system do not matter, suffice to say the model creates oscillating behavior as a rise in population of the prey creates a corresponding rise in the population of the predator which in turn reduces the prey which reduces the predator, etc.

If we were to start the system at some arbitrary point in this space, say 0.5 rabbits/acre and 0.4 foxes/acre and then we can integrate the system forward to get a trajectory representing the evolution of the system. In the following example, the system starts at the red star and then the blue arc is ticked off in equal intervals of time. As the rabbit population grows, so too does the fox population until there's a rapid turn around as the foxes over-hunt the rabbits causing the fox population to plummet. If we let it continue we see that it would get back to where we started and it would keep on cycling in this trajectory forever.

The idea of stability analysis is to examine such plots for patterns around the equilibria. The above case entails a fairly complex analysis, so let's start simpler and return to this case later.

Let's start instead with a simple linear system with a single stable equilibrium at the origin. No matter where we start the system, it will evolve into this equilibrium at the origin. This is a stable equilibrium as all the arrows point into the same place.

The above simple linear system is represented by:

where Z(t) is the 2-vector that is the state of the system and and M = [ -0.1 0; 0 -0.1 ]. Note that M in this case is a uniform negative scaling matrix.

Now consider the opposite: a uniform positive scaling matrix, we get the graph where all the arrows point outwards:

In this case, even if we start extremely close to the equilibrium, say at [0.01 0.01], the system diverges rapidly. This is an unstable equilibrium.

The following is the in-between case: stable in one direction, unstable in another. This is called a "saddle point" because of the shape of a saddle: curving up (from horse's nose to rear) and curving down (from side to side). This is M = [ +0.1 0; 0 -0.1 ] . In this case we have a change in convexity around the equilibrium so an approach from any direction towards the equilibrium causes a deflection. We can therefore see that only in the case that all directions point inwards will an equilibrium be stable.

So far, these systems are like the valley and the mountain analogies typically used. But now consider this linear system which isn't like those:

Here, M is a rotation matrix = [ 0 0.1; -0.1 0 ]. Here if we start anywhere we get a circular orbit. It's for this kind of linear system that the "ball-rolling" analogies fail. There is no "surface" shape that would make a ball roll around in an orbit maintaining whatever radii it started at. The "ball-rolling" analogies just don't make any sense but the quiver plot does.

Finally, for illustration, here's a case where there's rotation and a saddle point at the same time.

So, how can characterize these cases mathematically? Eigen decomposition.

First, we solve the system for the places where the d/dt is zero, i.e. the equilibria. For the above linear cases, this is extremely simple: the equilibrium is at the origin. Now, we look at the curvature of the system around that equilibrium. Again, for these linear cases, this is very simple: the curvature is the same everywhere. That is to say, no matter where you look, the percentage change of the arrows is the same.

Eigen decomposition tells us the components of curvature. Consider the following system which is unstable all directions away from the equilibrium but it is more curved (ie "more unstable") in one direction than another. In green are plotted the eigen vectors scaled by the eigen values. In this case, the eigen values are 0.1 and 0.3, both positive, telling is that this system is unstable in all directions. The eigen vector (ie. the direction of instability) doesn't really matter if we're just asking the simple question "is this stable or not?" but I've chosen to plot it anyways just for illustration -- as you can see, the eigen vectors point in the directions of maximum and minimum curvature.

The sign of the eigen values tells us about the local curvature. A system is only stable if all the eigen values are negative: i.e. all the arrows point inwards. If the system has rotation, then the eigen values will be complex valued. That's what complex value means: that there's rotation in the scale factor. For example, the following has eigenvalues 0.17 + 0.09i and 0.17 - 0.09i and you can see the rotation in the field. Note that the two eigen values are complex conjugates of each other. This is just saying that you can look at this as either clockwise or counter-clockwise rotation by just changing the sign. Complex eigen values will always come in conjugate pairs.

Another way to intuit the complex-valued eigen values is to consider the definition of eigen direction as those directions that don't change under the given transformation. Consider a rotation matrix. What direction is invariant under a rotation? Nothing in the plane of the rotation! Whatever vector is invariant must somehow be perpendicular to the plane -- and what do we call the direction perpendicular to any system? -- the (stupidly named) "imaginary" direction.

You can also see from this why the complex valued eigen values come in conjugate pairs -- both the out-of-paper and into-paper directions are valid rotational axes.

Now that we've seen positive, negative, and complex-valued eigenvalues for linear system, we're ready to look at a more complicated system. In a non-linear system there may be more than one equilibrium and the curvature may vary from place to place. However, the eigen-decomposition tool can only be used on linear systems so we look at the approximate linear system right next to the equilibrium. This is called "linearization".

For example, consider this plot for a dampening pendulum. There's two equilibria: one at the bottom (angle zero) and one at the top (angle pi). The one at the top is a shearing unstable point -- a tiny displacement either way will knock the pendulum away. By the way, the graph is periodic: read x axis 2*pi the same as zero. You can see that starting the pendulum higher near the top causes it to spin around a lot faster and more times as you'd expect.

It's taken me 2 weeks to write this part so I'll post it now and another entry later on linearization of complicated systems like this one and the stadium wave we're working on for the amorphous computing demo.

1 comment:

Stability Analysis said...

Fantastic work showing a simple explanation of the limits of the ball+hill analogy. It's very helpful introducing the concept, but when people get too fixated on thinking that way they can get the wrong idea about how the system should behave. Quiver plots are great for overcoming this!