Lesson13.mw
ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL
Lesson 13  Direction Fields and Phase Portraits
Prof. Douglas B. Meade
Industrial Mathematics Institute
Department of Mathematics
University of South Carolina
Columbia, SC 29208
URL: http://www.math.sc.edu/~meade/
Email: meade@math.sc.edu
Copyright 2001 by Douglas B. Meade
All rights reserved

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.C1 Example 1 (Saddle)
13.C2 Example 2 (Center)
13.C3 Example 3 (Node)
13.C4 Example 4 (Spiral)
13.C5 Eigenvalues and the Phase Portrait
13.C6 Representative Examples
13.D Phase Portraits for 2d Nonlinear Autonomous Systems
13.D1 Example 5
13.D2 Example 6 (The Nonlinear Pendulum)
13.D3 Example 7
13.E Nonautonomous Systems
13.E1 Example 8
13.F Stability of an Equilibrium Point
13.F1 Example 9
13.F2 Example 10
13.F3 Example 11
13.F4 Summary
13.G Linearization of Nonlinear Systems
13.G1 TangentLine Approximation
13.G2 TangentPlane Approximation
13.G3 Linearizing Autonomous Systems
13.G3A Example 12
13.G4 Interpretation Rules
13.G5 Example 13
Initialization
> 
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 twodimensional system of firstorder ODEs, drawn in the phase plane for the system, is similar to the direction field for a single firstorder 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 twodimensional 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.
13.C1 Example 1 (Saddle)
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 forloop structure. Each IVP will be solved, and the solution
, packaged in a list of the form
ready for parametric plotting.
> 
for k from 1 to 8 do
init := x(0)=Pk[1], y(0)=Pk[2];
qq := dsolve({q1,q2, init},{x(t),y(t)});
fk := eval([x(t),y(t),t=0..5], qq);
print(fk);
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(Pk,.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.C2 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.C3 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,5k,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.C4 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.C5 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
============ =============
SADDLE: Opposite signs
CENTER: Pure imaginary
NODE: Same sign
SPIRAL: Complex conjugates
===================================
Table 13.3
13.C6 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 =>
Saddle: 1,2 =>
Center:
=>
Spiral In:
=>
Spiral Out:
=>
=====================================================================
Table 13.4
13.D Phase Portraits for 2d Nonlinear Autonomous Systems
The twodimensional 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.D1 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
> 
q1 := [seq( seq( X+dX, dX=[ [h,h], [h,h], [h,h], [h,h] ] ), 
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.D2 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.D3 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.E1 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 threedimensional system
> 
sys4a := dchange( t=s, sys4, s ) union { diff( t(s), s ) = 1 }; 
Because this is a threedimensional 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 xy 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 xyt space (the phase space for the augmented system).
The

view is obtained with
> 
display( sol_traj4a, orientation=[90,0] ); 
and the ty 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.F1 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]): 
and compute its eigenvalues.
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.F2 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]): 
and compute the eigenvalues of 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.F3 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]): 
and compute the eigenvalues of 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 yaxis 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 halfplane
tend to
, whereas solutions starting in the halfplane
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.F4 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.G1 TangentLine Approximation
Near the point of contact, the curve
is approximated by its tangent line. For example, the function
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 pointslope formula
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 tangentline approximation is known as linearizing in one dimension.
13.G2 TangentPlane Approximation
The tangent plane approximates a function of two variables. This is the twodimensional 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 firstdegree polynomial approximation of
is just the multivariable Taylor polynomial of degree one.
To construct this tangent plane, first find the corresponding zcoordinate 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
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
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 firstdegree Taylor polynomial gives the linear approximation known geometrically as the tangent plane.
13.G3 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 higherorder remainder terms into the functions
and
, the firstdegree Taylor expansions
become
because
is an equilibrium point at which
= 0. The nonlinear autonomous system now reads
Define the change of variables
so that
and the system now reads
=
+
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 xyplane has been translated to the equilibrium point
in the uvplane. 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 xysystem at (
) can often be deduced by analyzing the behavior of the linear uvsystem at
, as will be shown shortly.
13.G3A 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
Pk := eval([x,y], q[k]);
od; 
The Jacobian matrix is
which is calculate using Maple's builtin 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
Ak := 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(Ak);
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.G4 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.G5 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
Pk := 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
Ak := eval(A, q[k]);
od; 
Obtain the eigenvalues of each of the matrices
.
> 
for k from 1 to 2 do
Eigenvalues(Ak);
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.
[Back to ODE Powertool Table of Contents]