Phase-Plane Analysis

phase.mws

Douglas B. Meade (meade@math.sc.edu)

1 Overview

Initial value problems for nonlinear equations — both first-order systems and higher-order equations — are very difficult to solve in closed form. However, significant amounts of information can often be obtained directly from the equations. The phase portrait is a simple, yet powerful, tool for uncovering some of this information. This introduction to phase-plane analysis will involve the examination of Examples 3, 4, and 5 from Section 5.7 of Nagle and Saff. While the worksheet follows the same general outline as the text, the worksheet both includes additional information and omits other information. This worksheet is intended for parallel use with the text.

2 Objectives

• understand what a phase portrait is

• be able to extract qualitative information from a phase portrait

• use the DEplot command to create phase portraits

3 New Maple Commands (and Options)

• dirgrid= optional argument for DEplot, used to specify a grid other than 20 x 20 (which is the default)

• DEplot3d command for creating three-dimensional plots from a system of ODEs

• linecolor= optional argument for DEplot, used to give special coloring schemes for solution curves • piecewise command for creating a piecewise-defined function

• scaling= optional argument for any Maple plot, used to force equal scaling on all axes in a plot

• scene= optional argument for DEplot, often used to create phase portraits • signum the sign function ( 1 for positive numbers, -1 for negative numbers,

and 0 for zero )

• stepsize= optional argument for DEplot, controls the number of points plottted on a solution curve (default is 20)

• view= optional argument for any Maple plot, specifies the viewing window

to be used for the plot

4 Example 3 ( Nagle/Saff, p. 290 )

Determine the trajectories and corresponding equilibrium solutions of

x’ = −y (y − 2), y’ = (x − 2) (y − 2)

4.1 Solution

> restart;
> with( plots ):
> with( DEtools ):

The Maple analysis begins by entering the system:

> SYS := [ diff( x(t), t ) = -y(t) * ( y(t) - 2 ),
> diff( y(t), t ) = ( x(t) - 2 ) * ( y(t) - 2 ) ];

SYS := [\frac{∂}{∂t} x(t) = −y(t) (y(t) − 2), \frac{∂}{∂t} y(t) = (x(t) − 2) (y(t) − 2)]

The critical points of this system were found, in Example 2, to be the single point (2,0) and the line y=2. The trajectories for this system are the general solutions to

   > ODE := subs( SYS, x(t)=x, y(t)=y(x),
   > diff( y(x), x ) = diff( y(t), t )/diff( x(t), t ) );

ODE := \frac{∂}{∂x} y(x) = −\frac{x − 2}{y(x)}

This equation, which is separable, is easy to solve by hand or with Maple:

   > dsolve( ODE, y(x) );

\frac{1}{2} y(x)^{2} + \frac{1}{2} x^{2} − 2 x = C1

4.1.1 Point to Ponder

Do you see how this solution is equivalent to the one given in the text?

The stability of the critical points is determined by the flow along each trajectory. One way in which this can be done is to create a phase portrait for this system. The phase portrait of a system is closely related to the direction field of a single equation – it shows how a particle moves from any point in space. The connection between direction fields and phase portraits is so close that the same Maple command is used for both: DEplot (from the DEtools package).

   > DEplot( SYS, [x,y], t=0..1, x=-2..6, y=-4..4, scene=[x,y],
   > title=`Example 3 (Section 5.7):  Phase Portrait` );

4.1.2 Be Patient!

This plot may take several minutes to appear. Fortunately, the current CPU speeds significantly reduce the delay.

4.1.3 Improving the Appearance of the Phase Portrait – the dirgrid option

Note that the fact that all points on the line y=2 are equilibria is not clear from the default (20 x 20) grid selected by Maple. A better picture can be obtained using the optional argument dirgrid=[21,21]; try it!

4.1.4 Obsolete Feature from Release 3

In Release 3 it was necessary to specify the system of ODEs in a special form before Maple would detect an autonomous system. The advantage of this was a drastic reduction in the time needed to create a phase portrait. It appears as though Release 4 is much better at detecting autonomous systems, and the special form is no longer supported.

For those who might be interested – or who might be needing to translate other Release 3 worksheets for use in Release 4 – here are the commands that would have been used in Release 3:

SYSa := [ -y*(y-2), (x-2)*(y-2) ]; DEplot( SYSa, [x,y], t=0..5, x=-2..6, y=- 4..4, scene=[x,y], title=‘Example 3 (Section 5.7): Phase Portrait‘ );

Note that the above DEplot command generates an error in Release 4.

Understanding the information presented in a phase portrait can be a little confusing — until you know what to read the graph. The discussion in the text should help with this – READ IT!. It also helps to be able to explicitly see the dependent variables plotted against time. The scene= argument of the DEplot command can be used to produce the plot of any combination of two or three variables for any initial condition. To illustrate, we will create the phase portrait and plot each component of the solution to this system with initial conditions x(0) = 4.5, y(0) = 0.

The initial condition can be specified as before, i.e., IC := [ 0, 4.5, 0 ], but the following form is somewhat more natural – particularly if only one initial condition is being specified. (The extra square brackets are necessary as the initial conditions must be a list of points or list of lists of conditions.)

   > IC := [ [ x(0)=4.5, y(0)=0 ] ];

IC := [[x(0) = 4.5, y(0) = 0]]

The phase portrait is essentially the same as the one created before. Note that the list of initial conditions replaces the ranges on x and y, the stepsize option is added to produce a smoother approximation to the solution curve, and the dirgrid option has been used as suggested previously. The linecolor option is used to give an indication of time along the curve (early points will be red, then moving through orange, yellow, green, cyan, blue, ...); see the on-line help for DEplot for a full description of these, and other, optional arguments.

   > Pxy := DEplot( SYS, [x,y], t=0..5, IC, scene=[x,y], stepsize=0.1,
   > dirgrid=[21,21], linecolor=t,
   > title=`Example 3 (Section 5.7):  Phase Portrait` ):

The corresponding plots of x vs. t and y vs. t are obtained by changing the right-hand side of the scene argument.

4.1.5 Bug in Release 4

The need for the arrows=NONE argument is a result of the fact that, in Release 4, DEplot uses arrows=SMALL as its default. This is a bug, since Maple should automatically use this as its default for any DEplot which includes the independent variable in the scene. This has been reported to Maple’s Technical Support; it has been fixed in the development version of Maple; the patch should be available in early 1997.

   > Ptx := DEplot( SYS, [x,y], t=0..5, IC, scene=[t,x], stepsize=0.1,
   > arrows=NONE, linecolor=t,
   > title=`Example 3 (Section 5.7):  x vs.  t` ):
   > Pty := DEplot( SYS, [x,y], t=0..5, IC, scene=[t,y], stepsize=0.1,
   > arrows=NONE, linecolor=t,
   > title=`Example 3 (Section 5.7):  y vs.  t` ):

The nicest way to view these plots is side-by-side. This can be achieved by displaying the plots as a 1 x 3 array; the scaling=constrained option ensures that the same scale is used for both axes.

   > display( array(1..3, [ Pxy, Ptx, Pty ] ), scaling=constrained
   );

4.1.6 Improving the Appearance of the Plot

It is usually necessary to expand the plot window within the Maple worksheet. To do this, click in the graphic region, then drag the anchors horizontally until the plot fill most of the horizontal width of the screen.

Notice how, during the first stage, the solutions for both x and y are decreasing with x>0 and y<0 – this corresponds to the quadrant of the phase plane with x>2 and y<0. The next stage of the solution takes place for x<2 and y<0; here x is decreasing and y is increasing. When y>0, x increases. The last stage of the solution is its approach to one of the critical points on the line y=2. Notice how the coloring of the curve indicates that the trajectory “slows down” as it approaches the line of equilibria at y=2.

4.1.7 Point to Ponder

How would the phase portrait and trajectories look for the solution starting with initial conditions

• x(0) = 3.5, y(0) = 0?

• x(0) = 3.5, y(0) = 3?

• x(0) = 4.5, y(0) = 3?

5 Example 4 ( Nagle/Saff, p. 293 )

Use a software package to construct a phase plane diagram about the critical point (0,0) for the linear system

x’ = −2 x + y, y’ = −5 x − 4 y.

5.1 Solution

> restart;
> with( plots ):
> with( DEtools ):

The Maple version of this system of ODEs is

> SYS := [ diff( x(t), t ) = -2*x(t) + y(t),
> diff( y(t), t ) = -5*x(t) - 4*y(t) ];

SYS := [\frac{∂}{∂t} x(t) = −2 x(t) + y(t), \frac{∂}{∂t} y(t) = −5 x(t) − 4 y(t)]

The discussion in the text begins with a view of the phase portrait for a variety of initial conditions on the rectangle -4 <= x <= 4, -5 <= y <= 5.

> IC := [ [0, 4, -5], [0, 4, 0], [0, 4, 5], [0, -4, 5],
> [0, -4, -5], [0, -4, 0], [0, 0, -5], [0, 0, 5] ];

IC := [[0, 4, −5], [0, 4, 0], [0, 4, 5], [0, −4, 5], [0, −4, −5], [0, −4, 0], [0, 0, −5], [0, 0, 5]]

> DEplot( SYS, [ x, y ], t=0..1, IC, scene=[x,y], stepsize=0.1,
> title=`Example 4 (Section 5.7):  Phase Portrait` );

5.1.1 Points to Ponder

The trajectories, in phase space, appear to spiral towards the origin.

• How would the trajectories look when plotted against time?

• Use the techniques from the previous example to check your conjecture.

To zoom in on the behavior close to the origin, all that is needed is to restrict the view to a rectangle, say, -0.1 <= x <= 0.1, -0.1 <= y <= 0.1. (Note that it is also necessary to increase the time interval.)

   > DEplot( SYS, [ x, y ], t=0..10, IC, scene=[x,y], view=[-0.1..0.1,-0.1..0.1],
   > stepsize=0.1, title=`Example 4 (Section 5.7):  Phase Portrait
   -- ZOOM` );

Note that the view option must change the default arrows setting to arrows=NONE (specify arrows=SMALL to include the direction field in the plot). This process can be continued (almost) indefinitely.

6 Example 5 ( Nagle/Saff, p. 294 )

The motion of a spring–mass system with “sticky friction” (see Nagle and Saff, p. 294) can be modelled by

m x”( t) = −k x( t) + F_{friction}

where the frictional force depends on the instantaneous velocity, v(t) = \frac{∂}{∂t} x(t), and the spring force, F_{spring} = −k x(t), relative to the maximum static friction, μ, that can be exerted.

6.1 Note: Specific from of F[friction]

[ Since there is no simple means of getting Maple to produce a typeset form of F[friction], the specific formula for F[friction] will be given only in the solution. ]

Describe the motion for the case where m= μ= k= 1.

6.2 Solution

> restart;
> with( plots ):
> with( DEtools ):

The Maple representation of this differential equation is

> ODE := m*diff( x(t), t$2 ) = -k * x(t) + F[friction];

ODE := m (\frac{∂^{2}}{∂t^{2}} x(t)) = −k x(t) + F_{friction}

where

> F[friction] := piecewise( diff( x(t), t ) <> 0, -mu*signum( diff(
x(t),t ) ),
> abs(k*x(t))<mu, k*x(t),
> abs(k*x(t))>=mu, mu*signum(x(t)) );

F_{friction} := \{ \begin{matrix}−μ signum(\frac{∂}{∂t} x(t)) & \frac{∂}{∂t} x(t) ≠ 0 \\ k x(t) & |k x(t)| < μ \\ μ signum(x(t)) & μ ≤ |k x(t)| \end{matrix}

6.2.1 Simplification – moving mass

Note, as indicated in the text, that the friction reduces to −μ signum(\frac{∂}{∂t} x(t)) assuming the mass is always moving. To implement this simplification, replace the previous command with (the much simpler)

F[friction] := -mu * signum( diff( x(t), t ) );

DEplot will accept second- (and higher-)order equations in the first argument, but will only create plots of solution curves – not phase portraits. To obtain a phase portrait, the equation must be converted to a first-order system of ODEs. In this example this is accomplished by introducing the velocity v= x’ as the second dependent variable. Thus, the problem can be re-written as

> SYS := [ diff( x(t), t ) = v(t),
> m * diff( v(t), t ) = subs( diff(x(t),t)=v(t), rhs(ODE) ) ];

SYS := [\frac{∂}{∂t} x(t) = v(t), m (\frac{∂}{∂t} v(t)) = −k x(t)+(\{ \begin{matrix}−μ signum(v(t)) & v(t) ≠ 0 \\ k x(t) & |k x(t)| < μ \\ μ signum(x(t)) & μ ≤ |k x(t)| \end{matrix})]

The specific parameters of interest at this time are

> PARAM := { m=1, mu=1, k=1 };

PARAM := {m = 1, μ = 1, k = 1}

so the corresponding system is

> SYS2 := subs( PARAM, SYS );

SYS2 := [\frac{∂}{∂t} x(t) = v(t), \frac{∂}{∂t} v(t) = −x(t) + (\{ \begin{matrix}−signum(v(t)) & v(t) ≠ 0 \\ x(t) & |x(t)| < 1 \\ signum(x(t)) & 1 ≤ |x(t)| \end{matrix})]

The direction field for this problem is now easily obtained

> DEplot( SYS2, [ x, v ], t=0..10, x=-5..5, v=-5..5, scene=[x,v],
> title=`Example 5 (Section 5.7) -- Direction Field` );

> IC := [ [ 0, 7.5, 0 ], [ 0, -7.5, 0 ], [ 0, 2.5, 0 ], [ 0, -2.5,
0 ] ];

IC := [[0, 7.5, 0], [0, −7.5, 0], [0, 2.5, 0], [0, −2.5, 0]]

> DEplot( SYS2, [ x, v ], t=0..10, IC, scene=[x,v],
> arrows=NONE, linecolor=[GREEN,BLUE,MAGENTA,CYAN],
> title=`Example 5 (Section 5.7) -- Direction Field` );

6.2.2 Points to Ponder

• The discussion in the text shows that the point ( -1/2, 0 ) is an equilibrium point.

• Does the phase portrait exhibit this feature? (Explain.)

• Do the trajectories ever reach the equilibrium, or are they only asymptotic?

• How would the direction field and trajectories change for different sets of parameters?

The phase portrait, together with the solution curves, can be produced exactly as in the first example of this worksheet.

> Pxv := DEplot( SYS2, [x,v], t=0..10, IC, scene=[x,v], arrows=NONE,
> linecolor=[GREEN,BLUE,MAGENTA,CYAN], stepsize=0.25,
> title=`Example 5 (Section 5.7):  Phase Portrait` ):
> Ptx := DEplot( SYS2, [x,v], t=0..10, IC, scene=[t,x], arrows=NONE,

> linecolor=[GREEN,BLUE,MAGENTA,CYAN], stepsize=0.25,
> title=`Example 5 (Section 5.7):  x vs.  t` ):
> Ptv := DEplot( SYS2, [x,v], t=0..10, IC, scene=[t,v], arrows=NONE,
> linecolor=[GREEN,BLUE,MAGENTA,CYAN], stepsize=0.25,
> title=`Example 5 (Section 5.7):  v=x' vs.  t` ):
> display( array( 1..3, [ Pxv, Ptx, Ptv ] ), scaling=constrained
);

For additional discussion of these solutions, please consult the text (p. 296 in Nagle and Saff).

6.2.3 Three-Dimensional Plots for ODEs (optional) To conclude, DEplot3d can be used to produce a three-dimensional plot that, when viewed with different orientations, can display either the phase plane or either of the two solution spaces. The syntax for DEplot3d is essentially the same as the one used for DEplot, with the exception that scene now has three elements; see ?DEplot3d for a full description of this command..

> DEplot3d( SYS2, [x,v], t=0..10, IC, scene=[t,x,v], stepsize=0.25,
> linecolor=[GREEN,BLUE,MAGENTA,CYAN],
> title=`Example 5 (Section 5.7):  3D View of Trajectories` );

6.2.4 Point to Ponder

Can you rotate this plot obtain each of the three views displayed in the previous plot.

7 Copyright

phase.mws – a Maple worksheet for use with phase-plane analysis (for use with section 5.7, p. 287, of Nagle and Saff) Prepared by: Douglas B. Meade (7 March 1996) as part of the Maple Supplement for Nagle/Saff (Addison-Wesley) Updated to Maple V, Release 4, by Douglas B. Meade (27 December, 1996) For more information, please contact the author at one of following addresses: Department of Mathematics, Univ. of South Carolina, Columbia, SC 29208 (803) 777-6183 meade@math.sc.edu (copyright, Douglas B. Meade, 1996)