Demo #12 
Nonlinear Systems & Chaos

In this demo, we begin by illustrating how one applies Mathematica's numerical capabilities to systems of nonlinear differential equations.  We illustrate using the nonlinear predator-prey system:

x'(t) = x(1-y)
y'(t) = y(x-1)

We begin by plotting the direction field obtained by eliminating t between the two ODEs:
	\(TraditionalForm\`\(dy/dt\)\/\(dx/dt\) = \(dy\/dx = \(y(x - 1)\)\/\(x(1 - y)\)\)\)
using PlotVectorField[{dx,dy}, ...] as in previous demos

<< Graphics`PlotField`

plt1 = PlotVectorField[ {x(1-y),y(x-1)}, {x,-0.2,4}, {y,0,4}, Frame->True]

This slope field suggests that the trajectories are closed; i.e., the critical point at (1,1) is a center. To confirm, we find numerical solutions for a series of initial conditions and overlay the plots. A convenient way to do this is to set up a table of solution curves:

curve = Table[ n, {n,2,4} ]

This table will ultimately hold plots of our solution curves. Next, we define the ODE system

deq1 = x'[t] == x[t] ( 1 - y[t] );
deq2 = y'[t] == y[t] ( x[t] - 1 );

and then set up a "Do Loop" to make the plots:

Do[
	Clear[x,y];
		soln = NDSolve[ {deq1, deq2, x[0]==1, y[0]==n},
							{x[t],y[t]}, {t,0,15} ];
		x = soln[[1,1,2]]; 
		y = soln[[1,2,2]];
		curve[[n-1]] = ParametricPlot[ {x,y}, {t,0,15} ],
	{n,2,4} ];

Notice that we have varied the initial condition on y from 2 to 3 to 4.

Finally, we overlay the numerical solutions on top of the slope field:

Show[ plt1, curve[[1]], curve[[2]], curve[[3]] ]

These closed periodic solutions clearly encircle the critical point (1,1), which is thus a center.

Next, we examine how self-limiting (logistic) growth of the prey species in a predator-prey model
damps periodic oscillations about the stable equilibrium point.

First, we define the model:

deq1 = x'[t] == x[t] (1 - y[t]) + ϵ x[t] (1 - x[t]);
deq2 = y'[t] == y[t] ( x[t] - 1);

Next, we pick a small positive value of epsilon and start a trajectory at x[0]=2, y[0]=2:

ϵ = 0.1;
ic1 = x[0] == 2;
ic2 = y[0] == 2;

The following command generates the numerical solution:

soln = NDSolve[ {deq1, deq2, ic1, ic2},{x[t],y[t]}, {t,0,75.} ];

Plotting the trajectory in the phase plane:

ParametricPlot[ Evaluate[ {x[t],y[t]} /. soln ], {t,0,50} ]

So the trajectory spirals in to the critical point at (1,1) as we expected. Next, lets look at the x(t) and y(t) trajectories:

Plot[ Evaluate[ x[t] /. soln ], {t,0,75} ]

Plot[ Evaluate[ y[t] /. soln], {t,0,75} ]

These damped oscillations are indicative of the spiraling approach to the critical point in the phase plane.

Chaos in Nonlinear Systems

Next, we turn to explore the subject of chaos, within the forced Duffing equation:

\(TraditionalForm\`x \*"""\  + \ x'\  - x\  + \ x\^3\  = \ 
        f\_0\)cos(t)

We will employ the initial conditions x(0) = 1, x' (0) = 0 and then systematically increase the amplitude of the forcing term, (Subscript[f, 0]). As a first step, lets define the first-order system corresponding to this second-order ODE and the IC's:

deq1 = x'[t] == y[t];
deq2 = y'[t] == - y[t] + x[t] - x[t]^3 + f0 Cos[t];
ic1 = x[0] == 1;
ic2 = y[0] == 0;

Let's now begin with a small value of f0 and then work our way up towards stronger forcing:

f0 = 0.6;
soln1 = NDSolve[ {deq1, deq2, ic1, ic2},
				{x[t],y[t]}, {t,0,200}, MaxSteps -> 5000 ];

Plot[ Evaluate[ x[t] /. soln1 ], {t,0,200} ]

Thus, there is an initial transient, then periodic behavior. To confirm, we plot the last 50 time units:

Plot[ Evaluate[ x[t] /. soln1 ], {t,150,200} ]

Plot[Evaluate[y[t] /. soln1], {t, 150, 200}]

Now let's make a plot of the x-y trajectories

plt2 = 
      ParametricPlot[
        Evaluate[{x[t], y[t]} /. soln1], {t, 150, 200}]

Next, we will repeat the above with a larger force, f0 = 0.7.

f0 = 0.7;
soln2 = NDSolve[ {deq1, deq2, ic1, ic2},
				{x[t],y[t]}, {t,0,200}, MaxSteps -> 5000 ];

Plot[ Evaluate[ x[t] /. soln2 ], {t,150,200} ]

Now, the dynamical trajectory no longer contains a single harmonic (frequency), and the period of the repeat has "doubled" to 4 Pi. This is known as "period-doubling". The y trajectory looks similar and the phase plot looks like:

plt3 = ParametricPlot[ Evaluate[ {x[t],y[t]} /. soln2 ], {t,150,200} ]

Evidently the trajectory now goes around both of these orbits, resulting in the doubled period.

This strange period doubling continues, with an 8 Pi period when f0 = 0.75. Things get even more complicated
after that, as we shall see!

Next, let's try f0 = 0.8:

f0 = 0.8;
soln3 = NDSolve[ {deq1, deq2, ic1, ic2}, 
		{x[t],y[t]}, {t,0,200}, MaxSteps -> 5000 ]

Plot[ Evaluate[ x[t] /. soln3 ], {t,100, 200} ]

Plot[Evaluate[y[t] /. soln3], {t, 100, 200}]

plt4 = ParametricPlot[ Evaluate[ {x[t],y[t]} /. soln3 ], {t,50,200} ]

Obviously, additional period doubling has occured and the phase space trajectory has gotten very messy, i.e. "chaotic". Successive period doubling is one established route to chaos, although much remains to be understood about the complex dynamical processes that are associated with the transition to chaos.

To finish this demo, let's do one more calculation with f0 = 1.0:

f0 = 1.0;
soln4 = NDSolve[ {deq1, deq2, ic1, ic2},
				{x[t],y[t]}, {t,0,200}, MaxSteps -> 5000 ]

Plot[ Evaluate[ x[t] /. soln4 ], {t, 150, 200} ]

plt5 = ParametricPlot[ Evaluate[ {x[t],y[t]} /. soln4 ], {t,50,200} ]

Now, the bizzare thing is that although we had nearly attained chaos at f0 = 0.8, by increasing the forcing strength to f0 = 1.0, the dynamics has returned to a simpler 6 Pi periodic motion that now encloses BOTH critical points
at (1,0) and (-1,0). Increasing f0 to about 1.10, however, returns us to chaos again!

We can increase out understanding of these solutions by plotting the direction field of the (homogeneous) ODEs we started with
	x'  =  y
y'  =  -y  +  x  - (x^3)
Note that there are critical points at (-1,0), (0,0), and (1,0) where x' = y' = 0.  Thus we plot the direction field with the critical points shows as blue circles.

plt6  =  
      Show[PlotVectorField[{-y1 + x1 - (x1^3), y1}, {x1, -1.5, 
            1.5}, {y1, -0.75, 0.75}, Frame -> True], 
        ListPlot[{{-1, 0}, {0, 0}, {1, 0}, {1, 0}}, 
          PlotStyle -> {AbsolutePointSize[6], Hue[0.6]}]]

Now, let's look at our first solution:

Show[plt6, plt2]

The single-orbit solution circles the critical point at (1,0) owing to our choice of initial conditions.  Let's look at the next solution:

Show[plt6, plt3]

Here, we see how the period doubled solution fits on the direction field.  Let's try the chaotic solution:

Show[plt6, plt4]

An finally, we plot our the solution at the highest forcing where the orbits become simpler and once again close, meaning that the soution is periodic.

Show[plt6, plt5]

The message from all of this is that the dynamics of systems of nonlinear differential equations is a very rich and interesting area.Many things remain to be understood and it is an active area of current research.