Nonlinear Systems

Initialization

<<Graphics`PlotField`

Autonomous systems

Consider the autonomous system:

    x'[t] == y[t] + x[t]y[t]
    y'[t] == x[t] + y[t]^2

Find all critical points of this system and identify the type and stability of each. Plot the vector field corresponding to this system and a representative sample of trajectories on a single graph.

Solution

F[x_, y_] = y + x y;
G[x_, y_] = x + y^2;

Find the critical points:

cp = Solve[{F[x,y] == 0, G[x,y] == 0}]

The critical points are thus {-1,-1}, {-1,1}, and {0,0}. To analyze them we need to examine the matrix of partial derivatives evaluated at each.

A = {{D[F[x,y],x], D[F[x,y],y]},
     {D[G[x,y],x], D[G[x,y],y]}};
A//MatrixForm

At the critical point {-1,-1}:

A1 = A /. cp[[1]];
A1//MatrixForm

Eigensystem[A1]

Since the eigenvalues have the same signs, the point {-1,-1} is an improper node. It is asymptotically stable since the eigenvalues are negative.
At the point {-1,1}:

A2 = A /. cp[[2]];
A2//MatrixForm

Eigensystem[A2]

The eigenvalues have the same signs, so the point {-1,0} is again an improper node. It is  unstable since the eigenvalues are positive.
At the point {0,0}:

A2 = A /. cp[[3]];
A2//MatrixForm

Eigensystem[A2]

The eigenvalues have opposite signs, so the critical point is a saddle point, necessarily unstable.
We use the option ScaleFunction->(Log[50#+1]&) in plotting the vector field to help the vectors show up better near the critical point.

v = PlotVectorField[{F[x,y],G[x,y]},{x,-2,2},{y,-2,2},
ScaleFunction->(Log[50#+1]&), Frame->True]

From the vector field plot, it appears that initial points with y coordinates equal to -2, 0, or 2 would produce a representative set of trajectories. Both First[] and Flatten[] are present below only to get rid of extra braces. Evaluate[] keeps Mathematica from complaining about not being able to compile the functions.

solns = Flatten[Table[{x[t],y[t]} /. First[NDSolve[
        {x'[t] == y[t] + x[t]y[t],
         y'[t] == x[t] + y[t]^2,
         x[0] == i, y[0] == j}, {x[t],y[t]},{t,-3,3}]],
        {i,-2,2,0.5}, {j,-2,2,2}],1];

p = ParametricPlot[Evaluate[solns],{t,-3,3}]

Show[p,v,PlotRange->{{-2,2},{-2,2}}]

Pendulum

Consider the equations describing the motion of a damped pendulum:

    x'[t] == y[t]
    y'[t] == -2.0 Sin[x[t]] - 0.75 y[t]

Here x[t] is the angle of the pendulum from the down position and y[t] = x'[t] is its angular velocity. Plot the vector field corresponding to this system and a representative sample of trajectories on a single graph. Then plot the trajectories for initial position down, theta[0] == 0, and the initial velocity equal to b == 3.0, to 5.0 in steps of 0.5. Experiment to find a value for the initial velocity b such that the trajectory approaches the unstable equilibrium point {Pi,0} corresponding to the pendulum being stationary in an upright position.

Solution

PlotVectorField[{y,-2.0Sin[x]-0.75y},{x,-2,8},{y,-2,5},
    ScaleFunction->(Log[50#+1]&), Frame->True]

solns = Table[{x[t],y[t]} /. First[NDSolve[
        {x'[t] == y[t],
         y'[t] == -2.0Sin[x[t]]-0.75y[t],
         x[0] == 0, y[0] == b}, {x[t],y[t]}, {t,0,20}]],
        {b, 3.0, 5.0, 0.5}];

ParametricPlot[Evaluate[solns],{t,0,20},
PlotRange->{{-1,8},{-2,5}}]

Here the critical velocity appears to be something less than 4.5. For initial velocities below this, the pendulum swings back and forth with ever decreasing angles x and velocity y coordinates. For initial velocities above 4.5, the pendulum swings around once, bringing all the angles to a range between Pi and 3Pi. The damping causes the pendulum to slow down and again the velocity decreases to 0 while the the angles approach 2Pi. Let's plot some solutions near the critical velocity.

ParametricPlot[Evaluate[{x[t],y[t]} /. First[NDSolve[
        {x'[t] == y[t],y'[t] == -2.0Sin[x[t]]-0.75y[t],
         x[0] == 0, y[0] == 4.4},
        {x[t],y[t]},{t,0,20}]]],{t,0,20},
        PlotRange->{{-1,8},{-2,5}}]

ParametricPlot[Evaluate[{x[t],y[t]} /. First[NDSolve[
        {x'[t] == y[t],y'[t] == -2.0Sin[x[t]]-0.75y[t],
         x[0] == 0, y[0] == 4.47},
        {x[t],y[t]},{t,0,20}]]],{t,0,20},
        PlotRange->{{-1,8},{-2,5}}]

ParametricPlot[Evaluate[{x[t],y[t]} /. First[NDSolve[
        {x'[t] == y[t],y'[t] == -2.0Sin[x[t]]-0.75y[t],
         x[0] == 0, y[0] == 4.467},
        {x[t],y[t]},{t,0,20}]]],{t,0,20},PlotRange->{{-1,8},{-2,4}}]

ParametricPlot[Evaluate[{x[t],y[t]} /. First[NDSolve[
        {x'[t] == y[t],y'[t] == -2.0Sin[x[t]]-0.75y[t],
         x[0] == 0, y[0] == 4.468465},
        {x[t],y[t]},{t,0,20}]]],{t,0,20},
        PlotRange->{{-1,8},{-2,5}}]

It looks like an initial velocity of 4.468465 is close to the critical value.

Undamped and Overdamped Pendulums

Undamped

v = PlotVectorField[{y,-Sin[x]},{x,-2,8},{y,-4,4},
    ScaleFunction->(Log[50#+1]&), Frame->True]

solns = Flatten[Table[{x[t],y[t]} /. First[NDSolve[
        {x'[t] == y[t],
         y'[t] == -Sin[x[t]],
         x[0] == a, y[0] == b},{x[t],y[t]}, {t,-4,4}]],
        {a,0,2Pi,2Pi},{b, -3.0, 3.0,0.5}],1];

p = ParametricPlot[Evaluate[solns],{t,-4,4},
PlotRange->{{-3,10},{-4,4}}]

Show[v,p,PlotRange->{{-2,8},{-4,4}}]

Overdamped:

v = PlotVectorField[{y,-Sin[x]-3y},{x,-2,8},{y,-3,3},
    ScaleFunction->(Log[50#+1]&), Frame->True]

solns = Flatten[Table[{x[t],y[t]} /. First[NDSolve[
        {x'[t] == y[t],
         y'[t] == -Sin[x[t]]-3y[t],
         x[0] == a, y[0] == b},{x[t],y[t]}, {t,-4,4}]],
        {a,-2,8},{b, -3.0, 3.0, 6.0}],1];

p = ParametricPlot[Evaluate[solns],{t,-4,4},
    PlotRange->{{-2,8},{-3,3}}]

Show[p,v,PlotRange->{{-2,8},{-3,3}}]]

Competing Species

Consider the competing species model:

    x'[t] == x[t]((3/2) - (1/2)x[t] - y[t])
    y'[t] == y[t]((3/4) - (1/8)x[t] - y[t])

Find all critical points of this system and identify their type and stability. Plot the vector field corresponding to this system and a representative sample of trajectories on a single graph. Determine the limiting behaviour as t goes to infinity and interpret the results in terms of the populations of the two species.

Solution

Clear[x,y]

F[x_,y_] = x((3/2) - (1/2)x - y);
G[x_,y_] = y((3/4) - (1/8)x - y);

The critical points are:

cp = Solve[{F[x,y]==0, G[x,y]==0}]

A = {{D[F[x,y],x], D[F[x,y],y]},
     {D[G[x,y],x], D[G[x,y],y]}};
A//MatrixForm

A1 = A/.cp[[1]];
A1//MatrixForm

The eigenvalues are 3/2and 3/4, so the critical point {0,0} is an unstable improper node.

A2 = A/.cp[[2]];
A2//MatrixForm

Eigenvalues[A2]//N

The eigenvalues are both negative, so the critical point {2,1/2} is an asymptotically stable improper node.

A3 = A/.cp[[3]];
A3//MatrixForm

The eigenvalues are -3/2 and 3/8 and have opposite signs, so the critical point {3,0} is an unstable saddle point.

A4 = A/.cp[[4]];
A4//MatrixForm

The eigenvalues are 3/4 and -3/4 and have opposite signs, so the critical point {0,3/4} is again an unstable saddle point.

v = PlotVectorField[{F[x,y],G[x,y]},{x,0,3.5},{y,0,1},
    ScaleFunction->(Log[100#+1]&), Frame->True,AspectRatio->1]

We'll choose initial points on lines above, below and parallel to y == -(1/3)x + 1.

solns1 = Table[{x[t],y[t]} /. First[NDSolve[
         {x'[t] == x[t]((3/2) - (1/2)x[t] - y[t]),
          y'[t] == y[t]((3/4) - (1/8)x[t] - y[t]),
          x[0] == i, y[0] == -(1/3)i + 1/3},
          {x[t],y[t]},{t,-10,20}]],
         {i,0.0,1.0,1/6}];

solns2 = Table[{x[t],y[t]} /. First[NDSolve[
         {x'[t] == x[t]((3/2) - (1/2)x[t] - y[t]),
          y'[t] == y[t]((3/4) - (1/8)x[t] - y[t]),
          x[0] == i, y[0] == -(1/3)i + 3/2},
          {x[t],y[t]},{t,-10,20}]],
         {i,1.5,4,1/2}];

p = ParametricPlot[Evaluate[Union[solns1,solns2]],
    {t,-10,20}]

Show[v,p, PlotRange->{{0,3.5},{0,1}}]

If there is only one population, then it will grow or shrink to its 'natural' size. For the x population that size is 3 and for the y population it is 3/4. Both of these points are unstable. If both populations are present they will each grow/shrink until they reach a "peaceful coexistence" at the levels x==2, y==1/2 which is an asymptotically stable node of the system.

Predator-Prey

Consider the predator-prey model:

    x'[t] == x[t]( 2 -  y[t])
    y'[t] == y[t](-1 + (2/3)x[t])

where x[t] represents the population of the prey and y[t] represents the population of the predator.
Find all critical points of this system and identify their type and stability.
Plot the vector field corresponding to this system and a representative sample of trajectories on a single graph. Explain how the populations vary for initial points close to the critical point inside the first quadrant. Superimpose the graphs of one pair x[t], y[t] as a functions of time t to see how the populations increa/decrease with respect to each other.

Solution

Clear[x,y]

F[x_,y_] = x( 2 - y);
G[x_,y_] = y(-1 + (2/3)x);

The critical points are:

cp = Solve[{F[x,y]==0, G[x,y]==0}]

A = {{D[F[x,y],x], D[F[x,y],y]},
     {D[G[x,y],x], D[G[x,y],y]}};
A//MatrixForm

A1 = A/.cp[[1]];
A1//MatrixForm

The eigenvalues are 2 and -1, so the critical point {0,0} is an unstable saddle point.

A2 = A/.cp[[2]];
A2//MatrixForm

Eigenvalues[A2]

The eigenvalues are imaginary, so the critical point {3/2,2} is a center or a spiral point and has indeterminate stability.

v = PlotVectorField[{F[x,y],G[x,y]},{x,0,3},{y,0,4},
    ScaleFunction->(Log[50#+1]&),
    Frame->True,AspectRatio->1]

solns = Table[{x[t],y[t]} /. First[NDSolve[
         {x'[t] == x[t]( 2 - y[t]),
          y'[t] == y[t](-1 + (2/3)x[t]),
          x[0] == i, y[0] == 2}, {x[t],y[t]},{t,0,10}]],
         {i,0,1.25,0.25}];

p = ParametricPlot[Evaluate[solns],{t,0,10},
        PlotRange->{{0,5},{0,5}}]

Show[v,p,PlotRange->{{0,3},{0,4}}]

The trajectories near the critical point {3/2,2} appear to be closed curves which implies that the populations of the predator and prey are periodic functions. Let's plot one pair x[t], y[t] close to the critical point to see how the populations increa/decrease with respect to each other.

Plot[Evaluate[solns[[3]]], {t,0,10}]

Predator-Prey

Consider the predator-prey model:

    x'[t] == x[t](   1   - (1/2)y[t])
    y'[t] == y[t](-(1/4) + (1/2)x[t])

where x[t] represents the population of the prey and y[t] represents the population of the predator. Find all critical points of this system and identify their type and stability.
Plot the vector field corresponding to this system and a representative sample of trajectories on a single graph. Explain how the populations vary for initial points close to the critical point inside the first quadrant. Superimpose the graphs of one pair x[t], y[t] as a functions of time t to see how the populations increa/decrease with respect to each other.

Solution

Clear[x,y]

F[x_,y_] = x(   1   - (1/2)y);
G[x_,y_] = y(-(1/4) + (1/2)x);

The critical points are:

cp = Solve[{F[x,y]==0, G[x,y]==0}]

A = {{D[F[x,y],x], D[F[x,y],y]},
     {D[G[x,y],x], D[G[x,y],y]}};
A//MatrixForm

A1 = A/.cp[[1]];
A1//MatrixForm

The eigenvalues are 1 and -1/4, so the critical point {0,0} is an unstable saddle point.

A2 = A/.cp[[2]];
A2//MatrixForm

Eigenvalues[A2]

The eigenvalues are imaginary, so the critical point {1/2,2} is a center or a spiral point and has indeterminate stability.

v = PlotVectorField[{F[x,y],G[x,y]},{x,0,1.5},{y,0,4},
    ScaleFunction->(Log[50#+1]&), Frame->True,AspectRatio->1]

solns = Table[{x[t],y[t]} /. First[NDSolve[
         {x'[t] == x[t](   1   - (1/2)y[t]),
          y'[t] == y[t](-(1/4) + (1/2)x[t]),
          x[0] == i, y[0] == 2}, {x[t],y[t]},{t,0,15}]],
         {i,0,0.5,0.125}];

p = ParametricPlot[Evaluate[solns],{t,0,15},
        PlotRange->{{0,1.5},{0,4}}]

Show[v,p]

The trajectories near the critical point {1/2,2} appear to be closed curves which implies that the populations of the predator and prey are periodic functions. Let's plot one pair x[t], y[t] close to the critical point to see how the populations increa/decrease with respect to each other.

Plot[Evaluate[solns[[3]]], {t,0,15}]

Predator-Prey

Consider the predator-prey model:

    x'[t] == x[t]( 2 - (1/2)x[t] - y[t])
    y'[t] == y[t](-2 + x[t])

where x[t] represents the population of the prey and y[t] represents the population of the predator. This model includes logistic restrictions on the growth of x[t], i.e., its natural growth rate is 2 but it has a natural 'saturation' level determined by the environment (in this case the saturation level is 2).
Find all critical points of this system and identify their type and stability.
Plot the vector field corresponding to this system and a representative sample of trajectories on a single graph. Explain how the populations vary for initial points close to the critical point inside the first quadrant. Superimpose the graphs of one pair x[t], y[t] as a functions of time t to see how the populations increa/decrease with respect to each other.

Solution

Clear[x,y]

F[x_,y_] = x( 2 - (1/2)x - y);
G[x_,y_] = y(-2 + x);

cp = Solve[{F[x,y]==0, G[x,y]==0}]

A = {{D[F[x,y],x], D[F[x,y],y]},
     {D[G[x,y],x], D[G[x,y],y]}};
A//MatrixForm

A1 = A/.cp[[1]];
A1//MatrixForm

A2 = A/.cp[[2]];
A2//MatrixForm

A3 = A/.cp[[3]];
A3//MatrixForm

Eigenvalues[A3]//N

The eigenvalues are complex with negative real part, so the critical point {2,1} is an asymptotically stable spiral.

v = PlotVectorField[{F[x,y],G[x,y]},{x,0,5},{y,0,3},
    ScaleFunction->(Log[50#+1]&),
    Frame->True,AspectRatio->1]

solns = Table[{x[t],y[t]} /. First[NDSolve[
         {x'[t] == x[t]( 2 - (1/2)x[t] - y[t]),
          y'[t] == y[t](-2 + x[t]),
          x[0] == i, y[0] == 2}, {x[t],y[t]},{t,-10,10}]],
         {i,0.25,1.75,0.25}];

p = ParametricPlot[Evaluate[solns],{t,-10,10},
        PlotRange->{{0,5},{0,3}}]

Show[v,p,PlotRange->{{0,5},{0,3}}]

Here are the graphs of a typical pair {x[t],y[t]}:

Plot[Evaluate[solns[[5]]], {t,0,10}]