In the following we do a complete analysis of a nonlinear
system. This is broken into a number of different
parts:
Let's consider the nonlinear system
which might model some sort of populations, or might have just been pulled out of thin air.
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.
As we saw in lab 5, we can also investigate the behavior of the system by looking at nullclines. These are just where the and equations are individually zero, and give where trajectories in the phase plane must be moving only vertically (for the null clines) or horizontally (for the ).
The nullclines are when
,
or, when
or
which is,
or .
The nullclines are when
which is, when
or
so
or .
We can plot these to see what they say:
where the blue curve and the axis are the nullclines and the red curves (including the axis) are the 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 nullclines) or vertically (for the nullclines), and that the equilibrium points are where we'd expect.
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
and
we need , , and 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:)
At (0,0), the coefficient matrix is
So the linearized system here is We can find the eigenvalues and eigenvectors of the coefficient matrix easily with Mathematica:
The eigenvalues are 2 and 3, and the corresponding eigenvectors are and 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):
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 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 direction and leaving along a line that is almost vertical.
Similarly,
Note that the point (0,3) is stable -- all trajectories actually approach it! Can you see this from the direction field?
Again:
And this, too, is a stable point.
One more time:
And this is again an unstable saddle.
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.
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 {}. 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...