ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Unit 13 -- Direction Fields and Phase Portraits

Outline of Unit 13

Initialization

> restart;

> with( DEtools ):

> with( plots ):

> with( linalg ):

> with( PDEtools ):

Warning, the name changecoords has been redefined

Warning, the name adjoint has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

>

Consider 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) ) };

>

The DEplot command is used to create a direction field for this system. The key differences 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:

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

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

> dir_field1;

>

This system appears to have several equilibrium solutions.

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

>

or, as a list of points,

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

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

> symbolsize=18, color=BLUE ):

> display( [dir_field1,equil_plot1], axes=BOXED );

>

From this diagram, it appears that the equilibrium solution at the origin is a source, the one at (0,5) is a sink, and the ones at (2,0) and (-4,9) are saddle points.

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

> 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

> 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.

> 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 ):

> display( [dir_field1, equil_plot1, phase_port1], axes=BOXED );

>

For the next example, consider the system of ODEs that represents Newton's Second Law of Motion for a pendulum in which is the angle between the stable hanging position and the pendulum's position at time (and is the velocity). That is, for every integer , , corresponds to the stable hanging position and , corresponds to the unstable position with the pendulum pointing straight up from the pivot.

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

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

>

The direction field for this system appears better when the dirgrid option is used to increase the resolution of the direction arrows and the axes are constrained to use the same scaling

> 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 ):

> dir_field2;

>

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:

> _EnvAllSolutions := true:

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

> _EnvAllSolutions := false:

>

The name _Z1 that appears in the result from solve is Maple's name for an arbitrary integer.

> about( _Z1 );

Originally _Z1, renamed _Z1~:

  is assumed to be: integer

>

On subsequent execution, the name might be _Z2 , _Z3 , .... In these cases, it will be necessary to modify all occurrences of this Maple-generated name.

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

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

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

> symbolsize=18, color=BLUE ):

> display( [dir_field2, equil_plot2] );

>

The direction field suggests different behaviors for the equilibria with even and odd integers. Solution curves will further elaborate on these differences. With initial conditions selected from the "inflow" portions of the window of the direction field, i.e. =10, -5 < < 0 and =-10, 0 < < 5, and from points on the -axis

> 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 is

> 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] );

>

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

>

For the next example, consider the previous example with an additional "damping" term :

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

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

>

To facilitate comparisons with the previous example, create the direction field with the same window, but with a slightly higher resolution:

> 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 ):

> dir_field3;

>

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

> _EnvAllSolutions := true:

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

> _EnvAllSolutions := false:

>

where _Z2 is a Maple-generated integer parameter.

> about( _Z2 );

Originally _Z2, renamed _Z2~:

  is assumed to be: integer

>

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

>

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

>

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

> symbolsize=18, color=BLUE ):

> display( [dir_field3,equil_plot3] );

>

For the phase portrait, consider the initial conditions from the top and bottom edges of the viewing window

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

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

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

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

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

> display( [dir_field3,equil_plot3,sol_traj3] );

>

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

>

Any non-autonomous 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 .

For this final example, consider the non-autonomous system

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

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

>

The conversion to an autonomous system is achieved with the help of the dchange command from the PDEtools package:

> 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] ]:

>

is obtained with the DEplot3d command:

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

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

> sol_traj4a;

>

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 non-autonomous 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 );

>