Application Center - Maplesoft

# Lesson 13: Direction Fields and Phase Portraits

You can switch back to the summary page by clicking here.

Lesson13.mw

ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Lesson 13 -- Direction Fields and Phase Portraits

Columbia, SC 29208

-------------------------------------------------------------------

 >

Outline of Lesson 13

13.A Direction Fields and Phase Portraits

13.B Equilibrium Points

13.C Phase Portraits for 2d Linear Autonomous Systems

13.C-2 Example 2 (Center)

13.C-3 Example 3 (Node)

13.C-4 Example 4 (Spiral)

13.C-5 Eigenvalues and the Phase Portrait

13.C-6 Representative Examples

13.D Phase Portraits for 2d Nonlinear Autonomous Systems

13.D-1 Example 5

13.D-2 Example 6 (The Nonlinear Pendulum)

13.D-3 Example 7

13.E Nonautonomous Systems

13.E-1 Example 8

13.F Stability of an Equilibrium Point

13.F-1 Example 9

13.F-2 Example 10

13.F-3 Example 11

13.F-4 Summary

13.G Linearization of Nonlinear Systems

13.G-1 Tangent-Line Approximation

13.G-2 Tangent-Plane Approximation

13.G-3 Linearizing Autonomous Systems

13.G-3A Example 12

13.G-4 Interpretation Rules

13.G-5 Example 13

 >

Initialization

 > restart;

 > with( DEtools ): with( plots ): with( plottools ): with( LinearAlgebra ): with( PDEtools ): with( Student[Calculus1] ): with( VectorCalculus ):

Warning, the name changecoords has been redefined

Warning, the assigned names arrow and translate now have a global binding

Warning, these names have been rebound: `&x`, ArcLength, CrossProduct, DotProduct, Tangent

Warning, the assigned names `<,>` and `<|>` now have a global binding

Warning, these protected names have been redefined and unprotected: `*`, `+`, `.`, D, Vector, diff, int, limit, series

 >

13.A Direction Fields and Phase Portraits

A direction field for a two-dimensional system of first-order ODEs, drawn in the phase plane for the system, is similar to the direction field for a single first-order ODE (see Lesson 1, Lesson 3, or Lesson 11).

For the autonomous system, the independent variable does not explicitly appear, and therefore has the representation

The direction field, a figure in the phase plane (which is the -plane), contains arrows whose direction at is determined by the vector .  The slope of the shaft of the arrow is given by

at any point for which .  If , the direction field will have an arrow pointing either straight up or straight down.  Since the direction field shows only the direction in which a trajectory moves, the length of the arrows in unimportant.

A graph of solution curves in the space of the dependent variables is a phase portrait for the system.  (See Lesson 1 and Lesson 11.)

 >

13.B Equilibrium Points

An equilibrium point in the phase plane is a point at which = 0.  Such a point represents an equilibrium solution, a constant solution where the system point sits and does not travel along an orbit.

For the autonomous linear system, we find the only equilibrium point to be the origin, , provided the system matrix is not singular.  In the two-dimensional case, there are four kinds of equilibrium points, and they have the names saddle, center,  node, and spiral.  Table 13.1 summarizes these names and briefly characterizes each type of equilibrium point.

=============================================

SADDLE:  Two paths enter, two leave, all others pass by

CENTER:  Orbits are closed paths around the equilibrium point

NODE:       All paths either enter or leave the equilibrium point

SPIRAL:    Each path spirals around the equilibrium point.

=============================================

Table 13.1

Section 13.C contains examples of the phase portraits for each type of equilibrium point listed in Table 13.1

 >

13.C Phase Portraits for 2d Linear Autonomous Systems

Corresponding to the four types of equilibrium points listed in Table 13.1, there are characteristic phase portraits.  The geometric properties of these distinctive phase portraits are related to the analytic behaviors of the solutions of the corresponding systems.

The following four examples illustrate each of these characteristic phase portraits.

 >

Consider the system of differential equations

 > q1 := diff(x(t),t) = x(t) + 2*y(t); q2 := diff(y(t),t) = 4*x(t) - y(t);

 >

for which solutions through the eight initial points

 > P1 := [1,1]; P2 := [-1,-1]; P3 := [5,-10]; P4 := [-5,10]; P5 := [6,-10]; P6 := [4,-10]; P7 := [-6,10]; P8 := [-4,10];

 >

will be obtained.  Something as repetitive as generating solutions for these eight points demands a for-loop structure.  Each IVP will be solved, and the solution , packaged in a list of the form

 > for k from 1 to 8 do init := x(0)=P||k[1], y(0)=P||k[2]; qq := dsolve({q1,q2, init},{x(t),y(t)}); f||k := eval([x(t),y(t),t=0..5], qq); print(f||k); od:

 >

Next, plot all eight trajectories on a single set of axes, forming Figure 13.1.

 > pp1 := plot([f||(1..8)], x = -10..10, y = -10..10, color=black, tickmarks=[[-10,-5,5,10],[-10,-5,5,10]], labels=[x,`y  `], labelfont=[TIMES,ITALIC,12]): display(pp1, title="Figure 13.1");

 >

Adding some arrows and labels by executing the following Maple commands, gives Figure 13.2.

 > pp2 := textplot({[3.5,.6,`A`], [.42,5.9,`B`], [-3.7,-.8,`C`], [ .6,-4.7,`D`]}): for k from 1 to 8 do pp||(k+2) := disk(P||k,.25,color=black): od: pp11 := textplot({[1.5,.7,`1`],[-.8,-2,`2`], [4,-11,`6`], [5,-11,`3`], [6,-11,`5`], [-6,10.6,`7`], [-5,10.6,`4`], [-4,10.6,`8`]}): pp12 := arrow([5,5],vector([.4,.4]),.2,.4,1,color=black): pp13 := arrow([-5,-5],vector([-.4,-.4]),.2,.4,1,color=black): pp14 := arrow([2.5,-5],vector([-.2,.4]),.2,.4,1,color=black): pp15 := arrow([-2.5,5],vector([.2,-.4]),.2,.4,1,color=black): pp16 := arrow([-9.408, -8.265],vector([-.4,-.4]),.2,.4,1,color=black): pp17 := arrow([-9.4, -10.4],vector([-.4,-.4]),.2,.4,1,color=black): pp18 := arrow([9.408, 8.265],vector([.4,.4]),.2,.4,1,color=black): pp19 := arrow([9.4, 10.4],vector([.4,.4]),.2,.4,1,color=black): display([pp||(1..19)],view=[-10..10,-11..11], scaling=constrained, labels=[x,y],labelfont=[TIMES,ITALIC,12], title="Figure 13.2");

 >

The directions on the trajectories correspond to the motion of the system point as the parameter increases along the curve.  These directions are obtained from the differential equations.  Table 13.2 indicates how directions are obtained for the paths through the four points labeled in Figure 13.2.

Table 13.2

The origin is a saddle for the differential equations in this example.  For the saddle, two opposing trajectories enter the equilibrium point, and two opposing trajectories leave it.  One such trajectory is called a separatrix (plural, separatrices).  Trajectories starting on one side of a separatrix remain on that same side.

Aside from the two separatrices that enter the equilibrium point, all other solutions move away from it.  Two trajectories appear to exit the equilibrium point, and the remaining trajectories first move towards it, then ultimately move away.

 >

13.C-2 Example 2 (Center)

The equilibrium point at the origin is a center for the differential equations

 > q5 := diff(x(t),t) = x(t) - 2*y(t); q6 := diff(y(t),t) = x(t) - y(t);

 >

A set of initial points is

 > inits := [seq([0,k,0],k=1..3)];

 >

The initial points have three values, the first, representing the initial time.  Here, that is , the time at which the system point is at that initial point whose and coordinates are the second and third value in the initial condition.

The resulting phase portrait, Figure 13.3, is

 > DEplot([q5,q6],[x(t),y(t)], t=0..2*Pi, x=-4..4, y=-4..4, inits, stepsize=.1, linecolor=black, arrows=medium, title="Figure 13.3");

 >

A center is surrounded by closed trajectories.

 >

13.C-3 Example 3 (Node)

The following linear system of differential equations has an equilibrium point at the origin, and it is a node.

 > q3 := diff(x(t),t) = 7*x(t) - 2*y(t); q4 := diff(y(t),t) = 15*x(t) - 4*y(t);

 >

With the following initial points,

 > inits := [seq([0,5-k,0],k=3..7)];

 >

the behavior of solutions near the node can be determined.

The phase portrait, Figure 13.4, is drawn via

 > DEplot([q3,q4],[x(t),y(t)], t=-10..1/10, x=-4..4, y=-8..8, inits, stepsize=.1, arrows=medium, linecolor=black, title="Figure 13.4");

 >

The equilibrium point at the origin is a node "out," sometimes called a repeller or source. Were the trajectories in this node to point inwards, towards the origin, the node would be a node "in," sometimes called an attractor or sink.

 >

13.C-4 Example 4 (Spiral)

The equilibrium point at the origin is a spiral for the differential equations

 > q7 := diff(x(t),t) = -x(t) - 2*y(t); q8 := diff(y(t),t) = x(t) - 3*y(t);

 >

A set of initial points is

 > inits := [seq([0,k,0],k=7..9), seq([0,k,0],k=-9..-7)];

 >

and the phase portrait, Figure 13.5, is

 > DEplot([q7,q8],[x(t),y(t)], t=-2..2, x=-12..12, y=-30..30, inits, stepsize=.1, linecolor=black, arrows=medium, title="Figure 13.5");

 >

The equilibrium point at the origin is an inward spiral.  Like the inward node, the inward spiral is sometimes called an attractor also.  Clearly, if the trajectories spiral outward, the phase portrait would be that of an outward spiral, and it would likewise be called a repeller.

 >

13.C-5 Eigenvalues and the Phase Portrait

For the linear system x' = Ax, the eigenvalues of the matrix A characterize the nature of the phase portrait at the origin.  These relationships are summarized in the Table 13.3.

Equilibrium Point                    Eigenvalues

============            =============

CENTER:                            Pure imaginary

NODE:                                Same sign

SPIRAL:                             Complex conjugates

===================================

Table 13.3

 >

13.C-6 Representative Examples

Table 13.4 below gives representative examples of the four types of equilibrium points, and their eigenvalues.

We distinguish between a proper node and an improper node.  For the moment, the improper node is characterized by repeated eigenvalues.  In the next lesson we will see that the improper node poses a minor difficulty when studying nonlinear systems.

Nodes and saddles have exponential solutions that tend either to zero or infinity.  Hence, trajectories for these systems will either tend towards or away from the origin, depending on the sign of the eigenvalues.  The center has trigonometric solutions that are the parametric representations of closed curves.  The spirals have solutions with the trigonometric terms of the center, but have exponential factors that force the solutions either towards or away from the origin.

Equilibrium point         Representative Eigenvalues                       Corresponding Solutions

===================================================================

Proper Node In:                      -1,-2                          =>

Improper Node In:                  -1,-1                          =>

Proper Node Out:                    1,2                            =>

Improper Node Out:                1,1                            =>

Center:                                                          =>

Spiral In:                                         =>

Spiral Out:                                           =>

=====================================================================

Table 13.4

 >

13.D Phase Portraits for 2d Nonlinear Autonomous Systems

The two-dimensional nonlinear autonomous system

can have more than one equilibrium point.  Each real solution of the simultaneous equations = is an equilibrium point, so the phase portrait for such a system is much more complicated than that of the autonomous linear system.  This additional complexity is illustrated in the following three examples.

 >

13.D-1 Example 5

A direction field for the nonlinear autonomous system of differential equations

 > sys1 := { diff( x(t), t ) = x(t) * ( 6 - 3*x(t) - 2*y(t) ),

 > diff( y(t), t ) = y(t) * ( 5 -   x(t) -   y(t) ) };

 >

appears in Figure 13.6.

 > dir_field1 := DEplot( sys1, {x(t),y(t)}, t=0..1,

 > x=-5..5, y=-5..10, scene=[x,y] ):

 > display(dir_field1, title="Figure 13.6");

 >

The DEplot command was used to create this figure.  The key differences in usage for this command are the specification of the scene and the dimensions of the viewing window; the "time" interval is required, but relevant only when solution curves are included in the plot:

This system appears to have several equilibrium solutions, namely,

 > equil_soln1 := solve( map(rhs,sys1), {x(t),y(t)} );

 >

which are easier to read as the list of points

 > equil_pts1 := map( subs, [equil_soln1], [x(t),y(t)] );

 >

Figure 13.7 shows these points as blue circles on the direction field of Figure 13.6.

 > equil_plot1 := pointplot( equil_pts1, symbol=CIRCLE,

 > symbolsize=18, color=BLUE ):

 > display( [dir_field1,equil_plot1], axes=BOXED, title="Figure 13.7" );

 >

Figure 13.7 suggests that the equilibrium solution at the origin is a source, the one at is a sink, and the ones at and are saddle points.

To create a phase portrait for this system it is necessary to specify an initial condition for each trajectory. To create a reasonable set of initial conditions quickly, start with the equilibrium points and a few additional points in the phase plane, thereby developing the list

 > pt_list := [ op(equil_pts1), [5,0], [1,10], [-2,10],

 > [0,.2], [0,.5], [0,-.5] ];

 >

Perturb these points a small amount in each component to obtain the final list

 > h := 0.1:

 > q1 := [seq( seq( X+dX, dX=[ [h,h], [h,-h], [-h,h], [-h,-h] ] ),

 > X=pt_list ) ];

 >

After converting each initial condition into a pair of equations, create the phase portrait, and display it, along with the direction field and equilibrium solutions, as Figure 13.8.

 > IC1 := map( p->[x(0)=p[1],y(0)=p[2]], q1 ):

 > phase_port1 := DEplot( sys1, {x(t),y(t)}, t=0..2, x=-5..5, y=-5..10, IC1, scene=[x,y], arrows=NONE, linecolor=green ):

 > display( [dir_field1, equil_plot1, phase_port1], axes=BOXED, title="Figure 13.8" );

 >

13.D-2 Example 6 (The Nonlinear Pendulum)

The equations

 > sys2 := { diff( x(t), t ) = y(t),

 > diff( y(t), t ) = -sin(x(t)) };

 >

formulate Newton's Second Law of Motion for a simple plane pendulum in which is the angle between the stable hanging position and the pendulum's position at time ;   is the angular velocity of the pendulum bob.  Thus, for every integer , the point corresponds to the stable hanging position and the point corresponds to the unstable position with the pendulum pointing straight up from the pivot.

For this system, the direction field in Figure 13.9 appears better with the dirgrid option used to increase the resolution of the direction arrows, and with the constrained  option used to impose equal scaling on both axes.

 > dir_field2 := DEplot( sys2, {x(t),y(t)}, t=0..1, x=-10..10, y=-5..5,

 > scene=[x,y], dirgrid=[25,25], scaling=constrained ):

 > display(dir_field2, title="Figure 13.9");

 >

Although solve finds only one equilibrium solution for this system

 > solve( map(rhs,sys2), {x(t),y(t)} );

 >

it is readily apparent that there is an infinite collection of equilibrium solutions. Maple can find the entire family of equilibria when the _EnvAllSolutions environment variable is set to true.  The result is then

 > _EnvAllSolutions := true:

 > temp_soln := solve( map(rhs,sys2), {x(t),y(t)} ):

 > _EnvAllSolutions := false: unassign('k'); Q := remove(has,indets(temp_soln),{t})[1]: equil_soln2 := eval(temp_soln, Q=k);

 >

(Note on syntax:  Instead of the parameter k, Maple returns a symbol such as _Z1 to denote an arbitrary integer.  Unfortunately, on subsequent executions of the solve command, this symbol will change sequentially to _Z2, _Z3, .  To avoid such complications, the operation of solving for the equilibrium solutions has been coded to return the symbol for the arbitrary integer.)

To plot the equilibrium solutions in the direction field, convert the general expression to a point

 > q1 := op(map( subs, [equil_soln2], [x(t),y(t)] ));

 >

and substitute appropriate integers to get the list of points

 > equil_pts2 := [ q1 \$ k=-3..3 ];

 >

The desired plot then appears in Figure 13.10.

 > equil_plot2 := pointplot( evalf(equil_pts2), symbol=CIRCLE,

 > symbolsize=18, color=BLUE ):

 > display( [dir_field2, equil_plot2], title="Figure 13.10" );

 >

The direction field suggests different behaviors for the equilibria at even and odd integer multiples of . Solution curves will help to further elaborate on these differences. With initial conditions selected along the -axis, and from the "inflow" portions of the window for the direction field, i.e., from and , , the totality of initial points is

 > IC2 := [ [x(0)=10,y(0)=i] \$ i=-4..-1,

 > [x(0)=-10,y(0)=i] \$ i=1..4,

 > [x(0)=2*i,y(0)=0] \$ i=-4..4 ];

 >

The combined direction field and phase portrait so generated is seen in Figure 13.11.

 > phase_port2 := DEplot( sys2, {x(t),y(t)}, t=0..15,

 > x=-10..10, y=-5..5, IC2, scene=[x,y],

 > arrows=NONE, linecolor=CYAN ):

 > display( [dir_field2, equil_plot2, phase_port2], title="Figure 13.11" );

 >

Thus, for any integer , the point is a center and is a saddle point.

 >

13.D-3 Example 7

If damping is added to the simple plane pendulum of Example 6, the system equations become

 > sys3 := { diff( x(t), t ) = y(t),

 > diff( y(t), t ) = -sin(x(t)) - y(t) };

 >

To facilitate comparisons with Example 6, create the direction field with the same window, but with a slightly higher resolution.  The result appears in Figure 13.12.

 > dir_field3 := DEplot( sys3, {x(t),y(t)}, t=0..1,

 > x=-10..10, y=-5..5, scene=[x,y],

 > dirgrid=[30,30], scaling=constrained ):

 > display(dir_field3, title="Figure 13.12");

 >

The differences in the direction fields are obvious. Nonetheless, the system has the same equilibrium solutions, namely,

 > _EnvAllSolutions := true:

 > temp_soln := solve( map(rhs,sys3), {x(t),y(t)} ):

 > _EnvAllSolutions := false: unassign('k'); Q := remove(has,indets(temp_soln),{t})[1]: equil_soln3 := eval(temp_soln, Q=k);

 >

where the integer has the same explanation as appeared in Example 6.

Again writing this solution as the point

 > q3 := op(map( subs, [equil_soln3], [x(t),y(t)] ));

 >

and generating the following list of equilibrium points,

 > equil_pts3 := [ q3 \$ k=-3..3 ];

 >

Figure 13.13 showing both the equilibrium points and the direction field, is generated.

 > equil_plot3 := pointplot( evalf(equil_pts3), symbol=CIRCLE,

 > symbolsize=18, color=BLUE ):

 > display( [dir_field3,equil_plot3], title="Figure 13.13" );

 >

For the phase portrait, select initial conditions from the top and bottom edges of the viewing window.  With the initial points

 > IC3 := [ [x(0)=i,y(0)=5] \$ i=-10..10,

 > [x(0)=i,y(0)=-5] \$ i=-10..10 ];

 >

the phase portrait in Figure 13.14 is generated.

 > sol_traj3 := DEplot( sys3, {x(t),y(t)}, t=0..5,

 > x=-10..10, y=-5..5, IC3,

 > scene=[x,y], arrows=NONE, linecolor=green ):

 > display( [dir_field3,equil_plot3,sol_traj3], title="Figure 13.14" );

 >

From the direction field and equilibrium solution, it appears as though the equilibrium solutions at that were centers in the undamped system are inward nodes (attractors, or sinks) in the damped system while the equilibria at are saddle points for both systems.

 >

13.E Nonautonomous Systems

Any nonautonomous system can be reformulated as an autonomous system via the introduction of the independent variable, say , and augmenting the system with the (trivial) differential equation

The strange appearance of this differential equation that results from using the same name for the independent variable and one of the unknown functions can be avoided by introducing a new name for the independent variable, say . Then all occurrences of in the original system are replaced with and the auxiliary ODE becomes

13.E-1 Example 8

The nonautonomous system

 > sys4 := { diff( x(t), t ) = y(t),

 > diff( y(t), t ) = -sin(x(t)) - sin(t) };

 >

can be converted to an autonomous system with the help of the dchange command from the PDEtools package.  This results in the three-dimensional system

 > sys4a := dchange( t=s, sys4, s ) union { diff( t(s), s ) = 1 };

 >

Because this is a three-dimensional system, it is not realistic to draw a direction field. Instead, a sample of solution curves with initial conditions

 > IC4a := [ [x(0)=1,y(0)=1,t(0)=0], [x(0)=1,y(0)=-1,t(0)=0],

 > [x(0)=-1,y(0)=1,t(0)=0], [x(0)=-1,y(0)=-1,t(0)=0] ]:

 >

obtained with the DEplot3d command, appears in Figure 13.15.

 > sol_traj4a := DEplot3d( sys4a, {x(s),y(s),t(s)}, s=0..10,

 > IC4a, scene=[t,x,y], stepsize=0.2 ):

 > display(sol_traj4a, title="Figure 13.15");

 >

This plot can be rotated in real time by using the mouse to "grab" and move a corner of the box.

Specific views of the solution can be obtained by explicitly listing the viewing coordinates. For example, the x-y view is obtained with

 > display( sol_traj4a, orientation=[0,90] );

 >

Note that the crossing of solutions in this view is not reason for concern -- this is not an autonomous system in x and y -- and the solution curves do not intersect in the x-y-t space (the phase space for the augmented system).

The - view is obtained with

 > display( sol_traj4a, orientation=[-90,0] );

 >

and the t-y view is obtained with

 > display( sol_traj4a, orientation=[-90,90] );

 >

As a concluding note, it should be mentioned that the DEplot3d command is capable of handling nonautonomous systems. A simpler means of plotting the solution trajectories in this example would be to use

 > IC4 := [ [x(0)=1,y(0)=1], [x(0)=1,y(0)=-1],

 > [x(0)=-1,y(0)=1], [x(0)=-1,y(0)=-1] ]:

 > DEplot3d( sys4, {x(t),y(t)}, t=0..10, IC4,

 > scene=[t,x,y], stepsize=0.2 );

 >

13.F Stability of an Equilibrium Point

An isolated equilibrium point for an autonomous system of differential equations is said to be

Orbitally (Neutrally) Stable if all solutions sufficiently close to the equilibrium point remain close to that equilibrium point as ,

Asymptotically Stable if it is orbitally stable, and in addition, all solutions sufficiently close to the equilibrium point actually tend towards that equilibrium point as ,

Unstable if at least one solution starting close to the equilibrium point does not remain close to that equilibrium point as .

Special note

Some texts drop the modifier in conjunction with the orbitally (neutrally) stable equilibrium point, choosing to call such a point a stable equilibrium point.  However, without insisting on some modifier in this case, it can be difficult to determine if the word "stable" is being used in a precise sense.  Students who intend the term "asymptotic stability" but who carelessly drop the modifier then end up using the term "stable" ambiguously.  They could have intended to use "stable" in the sense of orbital stability, or in the sense of asymptotic stability, improperly expressed.  Hence, to avoid this possible ambiguity, we will always modify the term "stable" with either asymptotic or with orbital.

Table 13.5 listing Examples 9 - 11, which will be used to illustrate the terms given above.

Example                             System Matrix A             Eigenvalues                 Phase Portrait

=================             ============          =========               ==========

(13.9)

+                         Figure 13.16

___________________________________________________________________________

(13.10)

Figure 13.17

___________________________________________________________________________

(13.11)

+ 1                          Figure 13.18

____________________________________________________________________________

Table 13.5

 >

13.F-1 Example 9

The system in Example 9, Table 13.5, with equations

 > q1 := diff(x(t),t) = y(t); q2 := diff(y(t),t) = -x(t);

 >

has an orbitally stable equilibrium point at the origin.  To see what this means, form the coefficient matrix, A,

 > (Aneg,b) := GenerateMatrix(eval([q1,q2],[x(t)=_x,y(t)=_y]),[_x,_y]):

 > A := -Aneg;

 >

and compute its eigenvalues.

 > Eigenvalues(A);

 >

The eigenvalues are a pure imaginary conjugate pair.  Hence, the phase portrait is that of the center.  The origin is the center of closed orbits that neither spread apart nor approach each other, the essential meaning of orbital stability.  Any solution starting near the origin remains near the origin.  It does not necessarily get arbitratily close to the origin, but it does not get far from the origin either.  The following phase portrait, Figure 13.16, is useful in this regard since is shows that solutions starting near the origin remain near the origin, tending neither towards, nor away from the origin.

To draw the figure, pick the initial points

 > inits := [seq([0,k,0],k=1..3)];

 >

and use Maple's DEplot command to obtain a phase portrait.

 > DEplot([q1,q2], [x(t),y(t)], t=0..2*Pi, inits, linecolor=black, stepsize=.1, scaling=constrained, arrows=medium, title="Figure 13.16");

 >

13.F-2 Example 10

The system in Example 10, Table 13.5, whose equations are

 > q3 := diff(x(t),t) = -8*x(t) + 3*y(t); q4 := diff(y(t),t) = -10*x(t) + 3*y(t);

 >

has an asymptotically stable equilibrium point at the origin.  To see what this means, form the coefficient matrix, A

 > (Aneg,b) := GenerateMatrix(eval([q3,q4],[x(t)=_x,y(t)=_y]),[_x,_y]):

 > A := -Aneg;

 >

and compute the eigenvalues of A.

 > Eigenvalues(A);

 >

The eigenvalues of the system matrix are of the same sign, and negative.  Hence, the phase portrait is that of the node in.  That makes the origin an attractor, and hence, all solutions will be drawn into the origin, making the system asymptotically stable.  Any solution that starts near the origin will tend towards the origin as .  For example, consider the solution starting at ( ), obtained via

 > q5 := dsolve({q3,q4, x(0)=a,y(0)=c},{x(t),y(t)});

 >

Extract the component solutions

 > xt := eval(x(t), q5); yt := eval(y(t), q5);

 >

and compute the limits

 > Limit(x,t=infinity) = limit(xt,t=infinity); Limit(y,t=infinity) = limit(yt,t=infinity);

 >

As promised, an arbitrary solution starting near the origin, not only remains near the origin, but tends to the origin as .  Hence, all solutions, not only the ones near the origin, tend towards the origin with increasing , as seen in the phase portrait, Figure 13.17.  To obtain this figure, pick the initial points

 > inits := [seq([0,-3,k],k=-3..3), seq([0,3,k],k=-3..3)];

 >

and use Maple's DEplot command to obtain

 > DEplot([q3,q4], [x(t),y(t)], t=0..1, inits, linecolor=black, stepsize=.1, scaling=constrained, dirgrid=[10,20], arrows=medium, title="Figure 13.17");

 >

 >

13.F-3 Example 11

The system in Example 11, Table 13.5, whose equations are

 > q6  := diff(x(t),t) = x(t); q7 := diff(y(t),t) = 2*x(t) - y(t);

 >

has an unstable equilibrium point at the origin.  To see what this means, form the coefficient matrix, A,

 > (Aneg,b) := GenerateMatrix(eval([q6,q7],[x(t)=_x,y(t)=_y]),[_x,_y]):

 > A := -Aneg;

 >

and compute the eigenvalues of A.

 > Eigenvalues(A);

 >

The eigenvalues are of opposite signs.  Hence, the phase portrait is that of the saddle.  There are two trajectories that will enter the origin, but all other solutions, no matter how close to the origin they begin, will tend towards infinity as increases.  Since at least one solution starting near the origin does not remain near the origin, the origin is an unstable equilibrium point.

Solutions starting on the y-axis will tend towards the origin.  In fact, such solutions are precisely

 > q8 := dsolve({q6,q7, x(0)=0,y(0)=a},{x(t),y(t)}): x1 := eval(x(t), q8); y1 := eval(y(t), q8);

 >

for which

 > Limit(x,t=infinity) = limit(x1,t=infinity); Limit(y,t=infinity) = limit(y1,t=infinity);

 >

Any other solution with , given by

 > q9 := dsolve({q6,q7, x(0)=a,y(0)=c},{x(t),y(t)}): x2 := eval(x(t), q9); y2 := eval(y(t), q9);

 >

will not remain near the origin, as shown by the following limits.  Indeed, Maple computes the limits

 > Limit(x,t=infinity) = limit(x2,t=infinity) assuming a>0; Limit(y,t=infinity) = limit(y2,t=infinity) assuming a>0;

 >

for , whereas if , Maple computes

 > Limit(x,t=infinity) = limit(x2,t=infinity) assuming a<0; Limit(y,t=infinity) = limit(y2,t=infinity) assuming a<0;

 >

As promised, all other solutions, no matter how near the origin they may have started, tend towards infinity.

As seen in the following figure, Figure 13.18, the phase portrait, that of a saddle, makes the instability of this equilibrium clear.

To obtain this figure, pick the following initial points

 > inits := [seq([0,k,2],k=-2..2), seq([0,k,-2],k=-2..2)];

 >

and use the DEplot command.

 > DEplot([q6,q7], [x(t),y(t)], t=-2..2, x=-5..5, y=-5..5, inits, linecolor=black, stepsize=.1, arrows=medium, scaling=constrained, title="Figure 13.18");

 >

Solutions starting in the half-plane tend to , whereas solutions starting in the half-plane tend to .

An equilibrium point that is not stable, is called unstable.  Thus, if at least one solution starting near the equilibrium point does not remain near it, the equilibrium point is unstable.

 >

13.F-4 Summary

Table 13.6 summarizes the relationships between phase portraits, eigenvalues, and stability.

Equilibrium Point           Eigenvalues                                 Conditions                                  Stability

===========           ========                             ============                  ===============

Node (in)                      real and negative                                                                 Asymptotically stable

Spiral (in)                       +                                                                        Asymptotically stable

Center                            +                                         pure imaginary                   Orbitally stable

Node (out)                    real and positive                                                                  Unstable

Spiral (out)                     +                                                                        Unstable

Saddle                           real, with opposite signs                                                       Unstable

==============================================================================

Table 13.6

 >

13.G Linearization of Nonlinear Systems

13.G-1 Tangent-Line Approximation

Near the point of contact, the curve is approximated by its tangent line.  For example, the function

 > f := (x-1)^2+2;

 >

and its tangent at the point appear in Figure 13.19.

 > Student[Calculus1][Tangent]( f, x=2, output=plot, title="Figure 13.19", view=[0..3,0..5], labels=[x,y], scaling=constrained );

 >

Near the point the function very much resembles its tangent line .  Analytically, the tangent line is the Taylor polynomial of degree 1, .

The tangent line can be constructed by the following calculations.  First, the slope at ,

 > m := eval(diff(f,x), x=2);

 >

then, the equation of the tangent line via the point-slope formula

 > Y := m*(x - 2) + 3;

 >

Analytically, the tangent line is a Taylor polynomial of degree 1, obtained in Maple via

 > y = convert(taylor(f, x = 2, 2), polynom);

 >

Replacing a curve with its tangent-line approximation is known as linearizing in one dimension.

 >

13.G-2 Tangent-Plane Approximation

The tangent plane approximates a function of two variables.  This is the two-dimensional version of linearization.

For example, at the point  where , the tangent plane approximating the function

 > f := 3*x^2 + 2*y^2 - 5*x -7*y -9;

 >

is .  This first-degree polynomial approximation of is just the multivariable Taylor polynomial of degree one.

To construct this tangent plane, first find the corresponding z-coordinate at , the function value .  Thus,

 > z1 := eval(f, {x = 2, y = 2});

 >

To construct a plane tangent at obtain a vector normal to the surface .  The general normal vector is

so compute

 > N := < -diff(f,x), -diff(f,y), 1 >;

 >

and evaluate it at to obtain

 > N1 := eval(N, {x = 2, y = 2});

 >

Once the normal for the tangent plane is known, write the equation of the plane as

 > q := N1 . = d;

 >

To determine the value of d, force the plane to go through the point .  Hence,

 > eval(q, {x = 2, y = 2, z = -13});

 >

from which the implicit form of the tangent plane is found to be

 > q1 := eval(q, d = -29);

 >

Explicitly solving this equation for yields

 > Z := solve(q1,z);

 >

That is one way to find the equation of the tangent plane.  Alternatively, at the point where , generate the multivariable Taylor polynomial of degree one.  In Maple, use the mtaylor command to create a linear approximation to the function , thereby getting the explicit form of the equation of the tangent plane, namely, .  This results in

 > z = mtaylor(f,[x=2, y=2], 2);

 >

The two computations of the equation of the tangent plane yield the same result.  Hence, the first-degree Taylor polynomial gives the linear approximation known geometrically as the tangent plane.

 >

13.G-3 Linearizing Autonomous Systems

Near an equilibrium point of the autonomous nonlinear system

the behavior of solutions can be studied by linearizing the functions and at each equilibrium point.

Let ( ) be an equilibrium point determined by the equations = 0.  Then, lumping the higher-order remainder terms into the functions and , the first-degree Taylor expansions

become

because is an equilibrium point at which = 0.  The nonlinear autonomous system now reads

Define the change of variables

so that

= +

or

u' = Au + P

where u = , the matrix A is the Jacobian Matrix

and the functions and are just the functions and under the change of variables .

The equilibrium point ( ) in the xy-plane has been translated to the equilibrium point in the uv-plane.  If the perturbation term P is dropped, the resulting system u' = Au is a linear system with an equilibrium point at .

The behavior of the nonlinear
xy-system at ( ) can often be deduced by analyzing the behavior of the linear uv-system at , as will be shown shortly.

 >

13.G-3A Example 12

The differential equations

model the dynamics of two species competing for the same resources.  Such models will be studied at length in Lesson 27.  For this model, the functions and are

 > F := x*(60 - 3*x - 4*y); G := y*(42 - 3*y - 2*x);

 >

The equilibrium points are obtained by solving the equations = 0, done in Maple via

 > q := solve({F=0, G=0}, {x,y});

 >

Extract and label these four equilibrium points as

 > for k from 1 to 4 do P||k := eval([x,y], q[k]); od;

 >

The Jacobian matrix is

which is calculate using Maple's built-in jacobian command.

Unfortunately, the Jacobian is defined mathematically as the determinant of the Jacobian matrix, but Maple returns just the matrix itself with the command

 > A := Jacobian([F,G],[x,y]);

 >

At each of the equilibrium points, the Jacobian matrix evaluates to

 > for k from 1 to 4 do A||k := eval(A, q[k]); od;

 >

For each of the four matrices obtained by evaluating A at an equilibrium point, obtain the eigenvalues.  Thus, compute the eigenvalues of each of the four matrices .

 > for k from 1 to 4 do Eigenvalues(A||k); od;

 >

To determine the signs of the eigenvalues from the matrix , obtain floating point equivalents of its exact eigenvalues.

 > evalf(Eigenvalues(A4));

 >

For the linearized systems we have

P1 = (0, 0) - node out

P2 = (20, 0) - saddle

P3 = (0, 14) - saddle

P4 = (12, 6) - node in

Rules are needed for deducing the nature of the equilibrium points of the nonlinear system from the information generated for the linear system.  These rules appear in the next section.

 >

13.G-4 Interpretation Rules

The following five rules summarize how to convert information about the linearized system into information about the nonlinear system.  They form a hierarchy, and are applied from top to bottom.

1.  Equilibrium points that are unstable or asymptotically stable in the linearized system will be the same in the nonlinear system.

2.  Saddles and spirals in the linearized system remain the same in the nonlinear system.

3.  Proper nodes ( 's real and distinct) in the linearized system remain the same in the nonlinear system.

4.  Improper nodes ( 's real and equal) in the linearized system can be nodes or spirals in the nonlinear system.

Such nodes are called improper because of the slight ambiguity in passing from the linearized to the nonlinear system.  However, instability and asymptotic stability will persist by Rule 1.  Thus, an improper outward node will be either an outward node or an outward spiral.  An improper inward node will be either an inward node or an inward spiral.  So even though there is a slight ambiguity, it is still possible to determine whether the nonlinear system is unstable or asymptotically stable.

5.  Centers in the linearized system can be centers or spirals in the nonlinear system

This is the difficult case.  Since centers are not asymptotically stable to begin with, the center that becomes a spiral can become an inward or an outward spiral.  Centers are neither unstable nor asymptotically stable, so Rule 1 does not apply.  Thus, the very nature of the behavior of the nonlinear system is in complete doubt if the linearized system has a center.  No decision can be made on the basis of the linearization and far more powerful tools are then necessary to determine the behavior of the nonlinear system at such an equilibrium point.

It is often exceedingly difficult to distinguish between a spiral and a stable center in a nonlinear system.  It is rarely clear if the spiral behavior is an artifact of numerical (roundoff) error in an approximating numerical scheme.  Analyzing the stability of such a point requires a great deal more mathematical skill than will be developed here.

 >

13.G-5 Example 13

To analyze the stability of the nonlinear autonomous system

obtain the equilibrium points and , and the Jacobian matrix

Evaluate A at each equilibrium point, and compute the corresponding eigenvalues, obtaining the matrices

and

and the respective eigenvalues , and + .

The point is a saddle, and unstable, for the linearized system; therefore, it is a saddle, and unstable, for the nonlinear system.

The point is a center for the linearized system; therefore, the behavior for the nonlinear system cannot be determined from the linearization.  The best that can be done with the present suite of tools is try to infer the behavior from the phase portrait in Figure 13.20, obtained below.  The phase portrait suggests that may well be a center in the nonlinear case.

To implement these calculations in Maple, enter the funtions and .

 > F := (x,y) -> 5*x*(1 - y): G := (x,y) -> y*(3*x - 12): 'F'(x,y) = F(x,y); 'G'(x,y) = G(x,y);

 >

Then, solve the equations x' = y' = 0 for equilibrium points, obtaining

 > q := solve({F(x,y) = 0, G(x,y) = 0}, {x,y});

 >

There are two equilibrium points that can now be carefully identified and labeled.

 > for k from 1 to 2 do P||k := eval([x,y], q[k]); od;

 >

Compute the generic Jacobian Matrix of first partial derivatives.

 > A := Jacobian([F(x,y),G(x,y)],[x,y]);

 >

Evaluate the Jacobian Matrix at each of the two equilibrium points.

 > for k from 1 to 2 do A||k := eval(A, q[k]); od;

 >

Obtain the eigenvalues of each of the matrices .

 > for k from 1 to 2 do Eigenvalues(A||k); od;

 >

Point P1 = (0, 0) is a saddle for the linearized system; therefore it is a saddle for the nonlinear system.

Point P2 = (4, 1) is a center for the linearized system; therefore the behavior for the nonlinear system cannot be determined from the linearization.  The best that can be done with the present suite of tools is to draw a phase portrait, hoping to infer the behavior from that plot.

For this, write the two differential equations.

 > q1 := diff(x(t),t) = F(x(t),y(t)); q2 := diff(y(t),t) = G(x(t),y(t));

 >

and select the initial conditions

 > inits := [seq([0,10,2*k],k=-3..3), seq([0,-3,2*k],k=-3..3)];

 >

Then, use the DEplot command to construct a phase portrait.

 > display(DEplot([q1,q2],[x(t),y(t)],t=-1..1, x = -3..12, y = -6..8, inits, stepsize=.01, linecolor = black, arrows=medium), labels=[`

 > x`,y], labelfont=[TIMES,ITALIC,12], title="Figure 13.20");

 >

The phase portrait suggests that (4, 1) may indeed be a center in the nonlinear case.

 >