> map(with,{DEtools,plots,plottools}):

Warning, new definition for translate

For the IVP:

x1' = f1(t, x1, x2, x3), x1(t0) = a1;

x2' = f2(t, x1, x2, x3), x2(t0) = a2;

x3' = f3(t, x1, x2, x3), x3(t0) = a3.

If all of the functions f1, f2, f3, and all (9) partial derivatives diff(f i, xj), for i = 1 to 3 and j = 1 to 3 are continuous on a box B in 4-space, and the point (a1, a2, a3, t0) is inside that box, then:

Existence: There is a solution on some t-interval containing t0 (solutions "live" in 4-space, [x1, x2, x3, t]).

Uniqueness: The IVP has at most one solution on any t-interval containing t0 (again as 4d curves).

Extension: There is a maximally extended (over some time interval) solution which "leaves" the box B at the "endpoints" of the t-interval.

Continuity/Sesitivity: Small changes in the data, [a1, a2, a3] and-or [f1, f2, f3], yield small changes in the solutions.

Here is a "box" in 3 space:

> display3d(cuboid([-1,-2,-3],[1,1,1]),style=patch,axes=framed,
view=[-4..2,-4..2,-4..2],labels=[`x`,`y`,`t`]);

Now we specialize to systems in two state variables, x(t) and y(t), and look at the notions of autonomous (f1 and f2 do not depend on time), equilibrium points (constant "solutions"), cycles (orbits which "loop back" onto themselves) and nullclines (curves where x' or y' are zero). Since autonomous ODEs have no dependence on time (in their rate functions), we see that orbits are either identical, or never meet (in the x-y plane).

>

I'll look at a somewhat complicated ODE in terms of these notions, namely:

x' = -0.6x^2 -x*y + 5x y' = -0.6y^2 -x*y + 5y

> ode:=diff(x(t),t)=-.6*x(t)^2-x(t)*y(t)+5*x(t),
diff(y(t),t)=-.6*y(t)^2-y(t)*x(t)+5*y(t);

Many of these notions depend on the rate functions, which I'll denote fx (for what x' is equal to) and fy (the rhs of the second differential equation in the system, or what y' is suppose to be equal to). We'll want to manipulate these as if x(t) and y(t) were numbers, so I'll use x and y, as opposed to x(t) and y(t) here.

> fx:= -.6*x^2-x*y+5*x;
fy:= -.6*y^2-y*x+5*y;

Nullclines are the points [x,y] where x' = 0 (a so called x-nullcline) and the set of points where y' = 0 (y-nullcline). Normally, you would plot fx = 0 and fy = 0, using implicitplot if needed and look for the zeros. Our rate functions are simple enough to see or solve for where they are zero, directly. We can factor an x out of fx, giving fx = x*(-0.6x - y + 5), so the vertical line x = 0 is an x-nullcline, as is the result of setting the second factor equal to zero and solving for y:

> solve(fx = 0,{y});

Similarly, we can get both y-nullclines at once with:

> solve(fy = 0, {y});

Since we normally plot curves in the form y = f(x), I chose to solve for y above. Plotting x = 0 is problematic, and can be done using a parametric plot [0, t, t = a..b], which means x = 0, y = t and t runs from a to b.

> xnull:=plot({[0,t,t=-5..10], -.6*x+5.},x=-5..10,
color=red,title=`x-nullclines`,thickness=3):

xnull;

As the plot below shows, Maple has trouble drawing a vertical line on its own, and gets a bit shaky where nullclines cross one another:

> implicitplot(fx=0,x=-5..10,y=-5..10,thickness=3);

Nullclines are NOT (normally) orbits (a plot of the points [x(t), y(t)] where x and y solve the ODE), and thus can be crossed by orbits. The red lines above have x' = 0, so orbits which do cross these lines, must be moving straight up or straight down (no left-right movement allowed) as they cross. Before explaining how one decides this direction, let's add the y-nullclines. Plotting y = 0 is much easier.

> ynull:=plot([0,-1.666666667*x+8.333333333],x=-5..10,
color=blue,title=`y-nullclines`,thickness=3):

ynull;

> display([xnull,ynull],
title=`nullclines`,
view=[-5..10, -5..10]);

>

The nullclines divide the xy-plane into regions, where x' and y' have fixed signs (always positive or always negative). If we factor fx and fy, one can see the signs quickly: fx = (x - 0)( -y - 0.6x + 5) = (x - 0)(-1)(y - (-0.6x + 5) ). I factored out a negative one, so that the equation of the x-nullcline could be forced to appear. The third factor is positive, whenever y is above the diagonal red line, and negative whenever y is below the line. Likewise, fy = y(-0.6y - x + 5) = -0.6*(y - 0)(y - (-1.67x + 8.33) ). Consider the 4 sided region in which [2, 2] lies. For fx, this region is to the right of x = 0, so (x - 0) is positive, and it is below the diagonal red line, so the third factor is negative, and fx = pos*(-1)*neg = positive. All orbits in this region have x' positive, so they will be moving to the right. For fy, this region is above y = 0 and below the blue diagonal line, so y' = fy is positive (neg*pos*neg). Thus, orbits in this region move up and right. If we switch to the triangular region above this (where [1,6] lives, for example), the only change is that we now have y values above the diagonal x-nullcline, so fy is still positive, but fx is negative. Orbits will move left and up in the triangular region. This suggests that orbits which cross the red diagonal x-nullcline, do so in "straight upward" fashion, as if you're going counter-clockwise around the unit circle, from just below [1,0] to just above [1,0], or as a clock would run backwards through 3'oclock. Check out a DEplot of the orbit through [2,4]:

> crossorbit:=
DEplot({ode},[x(t),y(t)],
t=-.5 ..0.5,[[x(0)=2,y(0)=4]],
x=0..4,y=0..8,scene=[x,y],stepsize=0.05,
linecolor=aquamarine,arrows=none):

display([crossorbit,xnull,ynull],title=`crossing`);

Because we know that orbits move left and up in the triangle and right and up in the quadrilateral region, we can tell that this orbit moves "counter-clockwise" as time passes.

To check the signs of x' = fx and y' = fy, you can also just evaluate fx and fy at any point [x,y] in a region, and see which directions orbits will travel. All points in a region must have fx and fy of constant sign (positive or negative), fx positive means right, fx negative means left, fy positive means up and fy negative means down. Here is what would happen in the triangular region containing [6, 1]:

> subs({x=6,y=1},fx);
subs({x=6,y=1},fy);

That says right and down. As you cross the blue diagonal y-nullcline (between the 4 sided region and triangular region spoke of above), orbits go from right and up to right and down. They must cross the blue line with y' = 0, i.e. a flat spot, so this is like a clock passing through 12 noon (in the clockwise direction).

Check out the orbit through [4,2]:

> display([DEplot({ode},[x(t),y(t)],
t=-.5 ..0.5,[[x(0)=4,y(0)=2]],
x=0..4,y=0..8,scene=[x,y],stepsize=0.05,
linecolor=aquamarine,arrows=none),xnull,ynull],
title=`crossing`);

You can "see it all" with a DEplot, which provides arrows, so long as the system is autonomous:

> all:=DEplot({ode},[x(t),y(t)],
t=-5 ..10,
x=-5..10,y=-5..10,scene=[x,y],stepsize=0.05,
arrows=MEDIUM,color=black):

all;

Here it is with the nullclines

> display([xnull,ynull,all]);

When an x-nullcline and a y-nullcline cross, we have both x' = 0 and y' = 0. If this happens at [a, b], then the point [a, b] is called an equilibrium point. The reason is that the pair of constant functions, x(t) = a and y(t) = b, satisfy the ODE (because x' and y' are zero, and fx and fy are also zero). To find equilibrium points, just solve the system of equations fx = 0 and fy = 0 (for numbers x and y):

> equilpts:=solve({fx = 0, fy = 0},{x,y});

equilpts :=

Look at the DEplot above. It is clear that [0,0] pushes orbits away, as all arrows "run outwards" from [0, 0]. For the other equilibrium points, it isn't so clear what happens. Many things can occur, like
 

• "getting drawn near" (but an orbit can NOT make it to an equilibrium point, because an equilibrium point is an orbit, and orbits don't cross);

• approaching and then veering off; or

• "circling".

Numeric techniques are often times too crude to accurately visualize this behavior. Many times, near an equilibrium point, x' and y' are so close to zero, that the orbit seems to stand still. The maximally extended solution will then leave "the box B" via the "t directions". In the plot below for example, the orbit starts at [4, 4.2], heads toward the [3.125, 3.125] equilibrium point, but then veers off toward the one at [0, 8.33333], where it hardly moves at all. This means that if we plotted x(t), it would level off near 0 (but never get there). A plot of y(t) would level off at 8.3333 and never quite make it either. In 3 space, [x, y, t], this solution would be heading nearly straight-up, with x close to 0, y close to 8.3333 and t marching up.

> display([DEplot({ode},[x(t),y(t)],
t=0 ..5,[[x(0)=4,y(0)=4.2]],
linecolor=aquamarine,
x=-1..5,y=0 .. 10,scene=[x,y],stepsize=0.05,
arrows=MEDIUM,color=black),
xnull,ynull],title=`nearby`,
view=[-1..5,0..10]);

Here is the 3d view:

> DEplot3d({ode},[x(t),y(t)],
t=0 ..10,[[x(0)=4,y(0)=4.2]],
linecolor=aquamarine,
x=-1..5,y=0 .. 10,scene=[x,y,t],stepsize=0.05);

If you drag the image around and look straight down on it, arranging the x and y axes in their usual 2-d positions first, you'll see the same plot as the direction field above.

>

Here, I add a bunch of starting points, most of which get drawn toward the equilibrium at [0, 8.3333]

> DEplot({ode},[x(t),y(t)],
t=0 ..5,
[seq([x(0)=-1+k*.5, y(0)=.8*(k*.5)+2],k=0..6),
seq([x(0)=k+.3,y(0)=6.3],k=-1..4)],
linecolor=aquamarine,
x=-2..4,y=0 .. 10,scene=[x,y],stepsize=0.05,
arrows=MEDIUM,color=black);

Here it is in 3d:

> DEplot3d({ode},[x(t),y(t)],
t=0 ..5,
[seq([x(0)=-1+k*.5, y(0)=.8*(k*.5)+2],k=0..6),
seq([x(0)=k+.3,y(0)=6.3],k=-1..4)],
linecolor=aquamarine,
x=-2..4,y=0 .. 10,scene=[x,y,t],stepsize=0.05);

>

When orbits come back to themselves, they are called cycles. Cycles occur often with periodic solutions. Here is a direction field which shows some cycles, and some orbits which appear not to be cycles:

> DEplot({diff(x(t),t)=-x(t)+y(t)*sin(y(t)),
diff(y(t),t) = y(t)-x(t)^2*cos(x(t))},
[x(t),y(t)],t=-2..2,
[seq([x(0)=1+.5*k, y(0)=-1.6*.5*k-4],k=0..10)],
x=1..6,y=-12..-4,linecolor=aquamarine,
stepsize=0.01);

Try problems 2a,b and 3a from section 5.2