Phase portraits of numeric approximations (think nonlinear systems)

Unfortunately, many nonlinear systems of differential equations can't be solved (by Mathematica, at least) in any reasonable sort of manner. You can use the critical points of the system (we are talking mainly about 2-dimensional systems here) along with the eigenvalues of the linear approximaiton to the system and its phase portrait to analyze these systems.  While the material given above is theoretically sufficient to draw these things, it turns out there are some annoying little issues that can really cuase problems if you aren't careful (and sometimes even if you are).

So, let's start with the following nonlinear autonomous system:

or where and

Let's find our critical points:

Now, to classify these critical points, we need to evaluate the Jacobian at each of them and find its Eigenvalues (at the origin).  So, we need to define a Jacobian first:

At each critical point, you need to find the Jacobian and then find the eigenvalues of that matrix.  So, for example, at the first critical point:

So, telll me, what does it mean that these critical points are imaginary with negative real part?  (Chcck your answer down below when we graph the phase portrait.)  Now, it is actually pretty easy to do this for all of the critical points at once (each row is a pair of eigenvalues, listed in the same order as the critical points in the list above):

	7.3898418826991215`
	3.750122581916463`

Take a moment to go through and classify each critical point.  Now, let's graph the vector field (well, the direction field, to be precise) for this system:

(Notice that we choose our domain to plot over so we can see the critical points and the region near them.  If your vectors are too long or too short, you can change the value of the ScaleFactor option.)

Let's graph our critical points:

(Notice two things here:  We use ListPlot since we just want to plot some points, and we choose a PointSize of about 2% to make the points big enough to see easily.)

Sometimes, it is also handy to graph your "null-clines" as well:

Now, it's time to put them all together:

This is pretty cool, but it would also be really helpful to add in some trajectories as well.  Unfortunately, there is really not a nice automatic way to pick out the most useful trajectories to plot.  Basically, you have to choose these individually (unless you just get lucky).  Basically, you pick points in interesting regions of the graph to be initial conditions.  Of course, since that means those would be the starting points of the graph, you might want to actually plot them for some negative t values as well as postive (i.e., backtrack along the curves to see where they came from).  

Now, if the system can be solved symbolically in a reasonable amount of time, this isn't all that difficult (see the examples earlier in this notebook).  However, if we try to solve this one, we get:

So no dice there.  Clearly, we need to use NDSolve to find the solution.  So, for example, if we want to plot a curve with initial condition we could do:

This sounds good.  But what about through the point ?

Uh, I think we have a problem here.  If you read the error messages, it is warning you that it can't compute a good solution all the way out to (it actually doesn't get much beyond ), so it tells you it is extrapolating to get all the way to 2.  As you can probably guess from looking at the numbers on the axes, this was pretty wildly unsuccessful.  You could fix this by simply graphing from to :

This works okay if you are only graphing a few curves, but it would really be nice to supply a list of points and then have Mathematica graph the solution curves for all of them, without having to set separate domains for individual curves.  The best way to do this would be to just instruct Mathematica to not extrapolate approximate functions beyond their domain.  Unfortunately, I haven't been able to figure out a way to do that (except on a case by case basis).  

So, after spending way too much time on this over the weekend, I decided to just write my own function to take care of this.  This was really annoying.  The good news for you is that you don't really have to worry about it.  You can just use the function I wrote (think cut-and-paste) and everything works out okay.  (You owe me for this one...)  Now, since I really didn't want to spend even more of my weekend making this really fancy, this means that you have to use this function exactly as I show you; I haven't built in any fancy error-handling or usage notes or anything else.  If you prefer not to use it, you can just plot individual curves (naming each one of them) and then combine them all together using the Show command.  On the other hand, if you prefer to take advantage of my function, I list it below (don't forget to execute it); do not change anything in it (unless you really know Mathematica pretty well) or you might break it.  (Unless you are just interested, it really isn't at all important that you understand how it works.  Actually, if you need to, you could change the value for PlotPoints, but be careful.)

So, how do you use it?  Fortunately, it is very simple.  You just have to supply it with the system (without any sort of initial conditions), a list of points you want curves through, the domain you want t to vary over, and the domains for x and y.  Like this:

Now, let's put all our graphs together:

Let's skip the null-clines and add in some more points (you should think about what points to choose to make things clearer):

Or, if you're just lazy (like me), you could have Mathematica plug in a grid of points for you (and if you're lucky, it will look good):

(You need the Flatten command here to knock out some excess {}'s.)

Well, you get the idea.  The combination of the direction field with the solution curves actually gives you a pretty decent idea of what's going on here.  Of course, you still have to figure out what that means, but that depends on the specific application.