Mathematics 304 – Ordinary Differential Equations Lab 1 – Direction Fields and Approximate Solutions for ODE September 5, 2000

Background

Last time, in Discussion 1, we introduced the direction field for the 1st order ODE \frac{dy}{dx} = f(x, y) and we saw one method for sketching the direction field by hand, using the isoclines (level curves of f(x, y)). Today, we want to use Maple to plot direction fields for several more involved examples and continue to study the qualitative information about solutions we can derive from the direction field.

Correction: In the handout for Discussion 1, the “comment” in part 2 of question B that solutions are defined for all x ∈ R is not true without extra hypotheses on f. I forgot to delete this after I modified an earlier version of the problem (dang!). Anyway, just cross that out, and be aware that there are autonomous equations where solutions are not defined for all values of the independent variable. We’ll see several examples in class.

Maple Commands for Plotting Direction Fields, Isoclines, and Solutions

The Maple commands that we will use today are called:

DEplot – used to plot direction fields and solutions

contourplot – used to plot isoclines

display – used to combine separate plots into a single plot

These commands are parts of two separate packages; to use them you will need to start by entering the commands:

                               with(plots):
                              with(DEtools):

to load the plots, DEtools routines. Note: the colon at the end of a Maple command suppresses the output – the command is executed but the results are not displayed. In this case the output would be (long) lists of the commands that are parts of the two packages – maybe interesting, but wasteful of paper when you go to print out your worksheet!

A Worked Example

For instance, suppose we wanted to study the direction field, isoclines, and solutions of various initial value problems for the 1st order ODE \frac{dy}{dx} = xy − 1. For instance, say we wanted to plot the solution of the initial value problem with y(0) = 1 for x ∈ [−1, 1] (this is chosen arbitrarily here). We could use the following commands. First we define the slope function f(x, y) = xy − 1:

                          f := (x,y) -> x*y - 1;

Then we set up the differential equation in Maple’s format:

                  eq := diff(y(x),x) = f(x,y(x));

Then we use DEplot to generate the plot we want:

         DEplot(eq,y(x),-1..1,[[y(0)=1]],linecolor=black);

When you execute the DEplot command the output will be a plot showing the direction field (with line segments drawn as arrows pointing in the direction of increasing x), together with a plot of the solution.

Comments:

1. The 3-step process above could actually be combined into a single command if you put the definition of the differential equation directly into the DEplot command (instead of building it up using the slope function). I recommend that you do things this way though, at least at first, to keep everything straight in your mind and minimize the chance for typing errors. We will also use the slope function in a different way in what follows, so having a separate definition for it will be useful.

2. The general format of the DEplot command is

      DEplot(equation,depvar(indepvar),range,inits,options);

where equation is the differential equation, in the format given above. After the equation comes the name of the dependent variable (i.e. the unknown function) with the independent variable in parentheses, then the range of values of the independent variable you want to see plotted, then the initial conditions, and options last. The order of these is fixed – you cannot mix them up freely. You can plot several different solutions together by putting more than one initial condition inside the outer pair of square brackets, separated by commas, e.g. [[y(0)=1],[y(0)=2],[y(-1)=2]]. Options can be used to control the way solutions are plotted and the appearance of the plot. The linecolor=black in the above is an option, for instance. Without that, Maple would use a default yellow color that does not show up when you print your worksheet(!)

2. To draw the solution, Maple is using an approximate numerical method similar to, but more powerful than, the Euler’s Method you may have studied in calculus (described in section 1.5 of our text). Methods of this type produce a table of values for an approximate solution function by “following the direction field”. Then Maple plots the approximate solution by connecting the points with straight line segments. This method is far from fool-proof! In particular, if the slope function is changing rapidly, discontinuous at some points, etc., it can lead to inaccurate results. Fortunately, there are ways to try to fine-tune its behavior. If you notice that your approximate solution looks very jagged (adjacent straight line segments have very different slopes) or if it jumps around wildly, you can try including an option

                           stepsize=.01

or some number even smaller than .01. The step size is the spacing between successive values of the independent variable when Maple computes the approximate solution. Reducing the stepsize from its default value can improve plots. It also slows down the computation, though, because more points must be computed! If the slope function is discontinous at some point on your solution, this may not help; in some cases, this type of method will always fail no matter how small the step size is!

Isoclines, and Combining Plots

As we said in class, the isoclines are the level curves or contours of the slope function. So to plot them in Maple, we could use the contourplot command you used in multivariable calculus. Continuing in the example above, enter the command

         contourplot(f(x,y),x=-1..1,y=-3..3,contours=[-3,-5/2,-2,
                 -3/2,-1,-1/2,0,1/2,1,3/2,2],color=blue);

to plot the contours of f(x, y) for the values c = −3, −5/2, … and color all the contours blue. Do you see why these have the shape they do here?

Now we want to combine the direction field, solutions, and isoclines into a single plot. To do this, we need to assign the results of the plot commands to variables in Maple:

         DSPlot:=DEplot(eq,y(x),-1..1,[[y(0)=1]],linecolor=black):

Note the colon at the end of this command. If you make this a semicolon, you will see what Maple’s internal representation of plot data structures looks like(!)

      CPlot:=contourplot(f(x,y),x=-1..1,y=-3..3,contours=[-3,-5/2,-2,
                 -3/2,-1,-1/2,0,1/2,1,3/2,2],color=blue):

The names DSPlot, CPlot could be anything, of course. Then to show the plots together:

                         display([DSPlot,CPlot]);

Note: You can combine any number of 2D plots together this way. Similarly for 3D plots. Maple usually does a reasonable job of combining the plots in an intelligent way. Here, since the plots we are combining cover different regions in the (x, y)-plane, it only shows the part of the contour plot in the region given by the plot of the direction field and solutions.

Lab Questions

Suggestions: To avoid confusion as you work on the worksheet for this assignment, use separate names for the slope functions, equations, etc. in the different problems. Also, to save some paper (and grading time), please do not show the worked example above or the intermediate plots produced on the way to the requested combined plots in the worksheet printout you submit.

A) Produce a combined plot as above showing the direction field, solutions for the three initial value problems

\{ \begin{matrix}\frac{dy}{dx} = (1 + x) sin(y) + x \\ y(−3) = 5; y(−3) = 2; y(−3) = 0; \end{matrix}

and isoclines for this ODE. Use -3..3 for the range of values of x in the DEplot command, and plot the contours on the region x=-3..3,y=-5..5. Choose your contour c values to make sure you are seeing a reasonable selection of isoclines in the plot.

B) Produce a combined plot as above showing the direction field, solutions for the initial value problems

\{ \begin{matrix}\frac{dy}{dx} = xy^{2} − x^{2}y \\ y(0) = −0.1; y(0) = 1; y(0) = 1.05; y(0) = 1.1; y(0) = 2; \end{matrix}

and isoclines for this ODE. Use -4..4 for the range of values of x in the DEplot command, and plot the contours on the region x=-4..4,y=-4..4. Choose your contour c values to make sure you are seeing a reasonable selection of isoclines in the plot. Your first approximate solution plots here will probably develop a bad case of the “jaggies” for x ≥ 3. That’s an artifact of the default choices Maple makes in its method for producing approximate solutions, not a property of the actual solutions. Adjust the options in the DEplot command to see that this behavior goes away if the step size is reduced. Explain why this might be expected. (What does the slope function look like at points near the x-axis, with x ≥ 3?)

C) The population P (t) (t in years) of a certain species of insect satisfies

\frac{dP}{dt} = (0.4)P (1 − P/30)

for 0 ≤ t ≤ 5. (This is called a logistic growth equation. It models the fact that the environment where the insects live might only be able to support a limited population.) At t = 5, a parasite enters the insects’ habitat and begins to kill a fixed proportion of the total population in each year. The growth law for the population changes to

\frac{dP}{dt} = (0.4)P (1 − P/30) − (0.25)P

This means that we effectively have the following ODE with a piecewise-defined slope function:

\frac{dP}{dt} = \{ \begin{matrix}(0.4)P (1 − P/30) & if t ≤ 5 \\ (0.4)P (1 − P/30) − (0.25)P & if t > 5 \end{matrix}

Use Maple to sketch the direction field and solutions for the initial value problems for this equation with P (0) = 1, 3, 25, 40 (this is four different initial conditions). The slope function can be defined like this in Maple:

fnC:=(t,P)->piecewise(t<=5,0.4*P*(1-P/30),t>5,0.4*P*(1-P/30)-0.25*P);

You do not need to include the isoclines in your plot on this problem. Use a range of t values t = 0..20. What seems to be happening to all the solutions here? We will study this equation in more detail on the first Problem Set.

Assignment

One Maple printout with solutions to the lab questions from each discussion group, due in class on Tuesday, September 12. Everyone in each group should work individually in the lab to make sure you are familiar with the Maple commands, then compare results and work up the version of the worksheet you want to submit to me together.