1. Harvesting the Seeds

Perhaps the first thing one wants to know about an autonomous system dx/dt = f(x, y), dy/dt = g(x, y) is where the equilibrium points are. For such a point, dx/dt and dy/dt are both zero, so the problem is the same as solving a pair of simultaneous equations. Standard root-finding methods can do this almost instantaneously provided they are given a starting value not too far away from the root. Here is an example. First we define the system that will be our main example, and then we appeal to FindRoot.

In[]:=

The root is returned in the form of a substitution rule. We can define rx and ry to be the x and y coordinates of the root, and then we plug into the given functions as a check (% refers to the preceding output).

In[]:=

So far, so good. But our problem here is a little different in that we want to know all the equilibrium points in a given rectangle. Here is where the power of Mathematica, and its openness in allowing us access to the data in a contour plot comes into play. Consider the following algorithm:

1. Generate the curve(s) representing f = 0 by using a contour plot for the 0 contour.

2. Evaluate g along the curves, and check for sign changes of g.

3. Since a sign change of g along f = 0 means that we are quite close to a root,
we can use the gathered values as seeds to FindRoot.

In fact, this procedure works very well and is not too difficult to program. The function Equilibria which accompanies this article implements it in detail. Since the main outline is so simple, we illustrate the programming methodology by following through an example. Graphics[ContourPlot[***]] returns a complicated graphics object that has imbedded in it the several Line objects making up the contours. Cases[***, _Line, Infinity] picks them out no matter where they are. First then strips off the Line header because we are really interested in these objects as lists. We reintroudce the Line headers in order to confirm visually that this scheme has worked.

In[]:=

Good. We now have the five curves as lists in contourData. Let us now focus on the central curve, given by contourData[[5]]. We apply g and look at the signs of the result.

In[]:=

Multiplying this vector by a shift of it causes the sign-changes to show up as -1s. There are six of them.

In[]:=

We now find the places in the plane where these sign-changes occurred.

In[]:=

This gives us six approximations. We define a function, root, that feeds a single seed to the rootfinder; mapping this function to the set of seeds produces six fairly exact roots.

In[]:=

The following figure, which shows the contours together with circles around the data we just found, shows that the approach has worked beautifully and we have found all six solutions on the central contour.

In[]:=

f[x_, y_] :=  x - y^2 Cos[y]
  g[x_, y_] := -y +  x Sin[x]
  contourData = Map[First,
    Cases[Graphics[
       ContourPlot[f[x, y], {x, -10, 10}, {y, -10, 10}, 
          Contours -> {0}, PlotPoints -> 50]],
       _Line, Infinity]];
  signs = Sign[Apply[g, contourData[[5]], 2]] // Short
  Rest[%] * Rest[RotateRight[%]] // Short
  Position[%, -1]
  seeds = contourData[[5, Flatten[%]]]
  root[{seedx_, seedy_}] := 
    FindRoot[{f[x, y] == 0, g[x, y] == 0}, {x, seedx}, {y, seedy}]
  Map[root, seeds]
  Show[ContourPlot[#[x, y], {x, -10, 10}, {y, -10, 10},
         Contours->{0}, ContourShading->False,
         PlotPoints ->80] & /@ {f, g},
       Graphics[{Circle[{x, y}, 0.32] /. %} ]];
  Show[ContourPlot[#[x, y], {x, -10, 10}, {y, -10, 10},
  Contours->{0}, ContourShading->False,
  DisplayFunction->Identity, PlotPoints->100, 
  ContourStyle->AbsoluteThickness[.6]] & /@ {f, g},
  
  Graphics[{AbsoluteThickness[.75],
      Circle[{x, y}, 0.32] /. %%} ],
     DisplayFunction->$DisplayFunction,
     PlotRange->All, AspectRatio->Automatic];

If your output will be a sequence of graphs, select a display option.
Display as:

As mentioned, Equilibria implements this in general. It finds all 16 equilibria points for our example in short order.

In[]:=

There are some options to Equilibria. Setting SeedsOnly to True causes just the seeds to be returned. The EquilibriaPlotPoints option controls the number of plot points in the contour plots that get the algorithm started; for difficult examples this might have to be increased from the default of 49. And there is a ShowEigenvalues option that produces a table of the equilibrium points together with the eigenvalues of the Jacobian at the points.

EXERCISE. Modify these ideas to get a routine that finds all complex roots of a complex function f(z) in a given rectangle. Just apply the same method, looking at f's imaginary part along a contour plot of the real part. FindRoot, does in fact work on complex functions, provided I appears in either the function or the initial value.

Up to New Visualization Ideas for Differential Equations