Nonlinear system analysis

In the following we do a complete analysis of a nonlinear system.  This is broken into a number of different parts:

This, loosely, gives a sense of how we can analyze such systems.  I've chosen to do all of this in Mathematica, but much of it can be done by hand, and all of it could be done in other packages like Maple if we were so inclined.  The important thing here is not how I did things in Mathematica, but rather what the methods and steps are, and what they tell us about the system.

General analysis without solving much of anything

Let's consider the nonlinear system
which might model some sort of populations, or might have just been pulled out of thin air.

direction field and equilibrium points

To get a quick sense of what this does, let's first generate a direction field.  I'm using Mathematica, so this involves the commands



We see that the points (0,0), (2,0), somewhere around (0, 2.8), somewhere around (1.8, 1.25), and maybe also somewhere around (0.5, 2) are interesting.  We suspect that these are the critical points (equilibrium solutions, because it's autonomous) of the system.  To see, we solve the system.  Numerical (decimal) values are easier to understand than roots of roots, so I'll use NSolve rather than Solve:


So we were right: the critical points (equilibrium solutions) are
    (0, 0), (2, 0), (0, 3), (1.7221, 1.1021), and (0.5807, 1.8979).
(There are a couple of other points, but they are complex-valued, so they don't have any physical significance.)  To get a better idea of what is going on around these points, we need to linearize around them and find the behavior there.  Before doing that, though, let's look at the system in a different way.

back to top

nullclines and equilibrium points

As we saw in lab 5, we can also investigate the behavior of the system by looking at nullclines.  These are just where the [Graphics:mmaimages/nlsys2_gr_8.gif] and [Graphics:mmaimages/nlsys2_gr_9.gif] equations are individually zero, and give where trajectories in the phase plane must be moving only vertically (for the [Graphics:mmaimages/nlsys2_gr_10.gif] null clines) or horizontally (for the [Graphics:mmaimages/nlsys2_gr_11.gif] ).

The [Graphics:mmaimages/nlsys2_gr_12.gif] nullclines are when
or, when
    [Graphics:mmaimages/nlsys2_gr_14.gif] or [Graphics:mmaimages/nlsys2_gr_15.gif]
which is,
    [Graphics:mmaimages/nlsys2_gr_16.gif] or  [Graphics:mmaimages/nlsys2_gr_17.gif].

The [Graphics:mmaimages/nlsys2_gr_18.gif] nullclines are when
which is, when
    [Graphics:mmaimages/nlsys2_gr_20.gif] or  [Graphics:mmaimages/nlsys2_gr_21.gif]
    [Graphics:mmaimages/nlsys2_gr_22.gif] or  [Graphics:mmaimages/nlsys2_gr_23.gif].

We can plot these to see what they say:



where the blue curve and the [Graphics:mmaimages/nlsys2_gr_26.gif] axis are the [Graphics:mmaimages/nlsys2_gr_27.gif] nullclines and the red curves (including the [Graphics:mmaimages/nlsys2_gr_28.gif] axis) are the [Graphics:mmaimages/nlsys2_gr_29.gif] nullclines.  Note that the equilibrium points are



So the equlibrium points are, as we expect, at the intersection of the nullclines.  It's interesting to superimpose this plot with the direction field



Note that the direction field arrows do cross the nullclines exactly horizontally (for the [Graphics:mmaimages/nlsys2_gr_35.gif] nullclines) or vertically (for the [Graphics:mmaimages/nlsys2_gr_36.gif] nullclines), and that the equilibrium points are where we'd expect.

Important Points:
back to top

Linearization at the different equilibrium points

Next, we might want to get a better picture of the behavior of the system at each of the equilibrium points.  Let's take these in turn and linearize the system around them.  To do this, we find the Jacobian and evaluate it at each equilibrium point.  Letting
    [Graphics:mmaimages/nlsys2_gr_37.gif]  and  [Graphics:mmaimages/nlsys2_gr_38.gif]
we need [Graphics:mmaimages/nlsys2_gr_39.gif], [Graphics:mmaimages/nlsys2_gr_40.gif], [Graphics:mmaimages/nlsys2_gr_41.gif] and [Graphics:mmaimages/nlsys2_gr_42.gif] to evaluate the jacobian.  As long as I'm in Mathematica, I'll go ahead and do this with it (fx, fy, gx and gy are, respectively, the derivatives we want):


Now, let's plug in the various points.

(Let's make all of the arrows that we draw around the equilibrium points be of length 0.375:)


linearization at (0, 0)

At (0,0), the coefficient matrix is


So the linearized system here is [Graphics:mmaimages/nlsys2_gr_47.gif]  We can find the eigenvalues and eigenvectors of the coefficient matrix [Graphics:mmaimages/nlsys2_gr_48.gif] easily with Mathematica:


The eigenvalues are 2 and 3, and the corresponding eigenvectors are [Graphics:mmaimages/nlsys2_gr_51.gif] and [Graphics:mmaimages/nlsys2_gr_52.gif]  Thus, around (0,0), the behavior is given by
We can illustrate this with arrows along the eigenlines (which go in the direction of the eigenvectors):


back to top

linearization at (2,0)

Continuing as before, at (2,0), the coefficient matrix of the linearized system is


So that the solution of the system is determined by the eigenvalues and eigenvectors of this matrix,


So the solution near (2,0) is given by [Graphics:mmaimages/nlsys2_gr_61.gif]  The eigenlines in this case therefore appear something like



So the point (2,0) appears to be a saddle point, with trajectories converging towards it in the [Graphics:mmaimages/nlsys2_gr_64.gif] direction and leaving along a line that is almost vertical.

back to top

linearization at (0,3)




Note that the point (0,3) is stable -- all trajectories actually approach it!  Can you see this from the direction field?

back to top

linearization at (1.7221, 1.1021)




And this, too, is a stable point.

back to top

linearization at (0.5807, 1.8979)

One more time:



And this is again an unstable saddle.

back to top

putting all of the linearizations together

Each of the linearized problems give us the behavior near one of the equilibrium points.  To see how all of this fits together, let's add them all together with the equilibrium points and nullclines (let's also recolor the equilibrium points as green for unstable and blue for stable):



Could we at this point guess the behavior of different trajectories in the plane?  Sure!  Note that because there are two attracting equilibrium points, however, there will be some trajectories which we won't be able to predict.  To put everything that we've done so far together, let's add the direction field to the above plot:



Be sure that you can tell how each of the different parts of our analysis (finding the direction field, finding equilibrium points, finding nullclines, and the solution of the linearizations at each equilibrium points) tell us similar and different things about the behavior of the system.

Important Points:
back to top

Numerical solution for trajectories

Because this is a nonlinear system, we can't actually solve it by hand.  So if we want to see what trajectories look like, we have to solve it numerically.  Mathematica has a pretty good numerical solver, so we'll just use that a couple of times here (for different initial conditions) and plot the resulting trajectories.


This is the list of different initial conditions that we'll use: each pair, e.g., {0.1,0.05}, is a pair of values for {[Graphics:mmaimages/nlsys2_gr_89.gif]}.  So we have 8 different sets of initial conditions in our iclist.


This is just a convenient way of solving the differential equation (with Mathematica's numerical solver, NDSolve) and plotting the resulting trajectory.  By defining a function like this, I save a lot of typing.


Ok, drum roll!  Let's see them:


So: if we start at (0.1,0.05) (the initial condition closest to the origin: it's marked below with a red dot)



We start at (0.1,0.05), and then move up towards the equilbrium point at (0.58,1.89) -- but that's a saddle, so we then swoop away from it and head to the stable equlibrium point at (1.72,1.1).

We can similarly discuss the other 7 trajectories.  Let's put them all together with the linearized solutions and equilibrium points:



And then, just for good measure, let's also add in the direction field and nullclines:



And how beautiful it is...

back to top

Converted by Mathematica      November 29, 2001   (with substantial subsequent editing)