2. Visualizing Solutions to ODEs
using DFIELD and PPLANE

§2.1 Introduction

In this lab, we are going to investigate differential equations of the form

(1)

dy/dx = f(x,y),

using direction fields, see Chapter 2 of BdeP for a thourough treatment. The direction field of the above differential equation is a diagram in the (x,y)-plane in which there is a small line segment drawn with slope f(x,y), at the point (x,y). For example, the direction field of the differential equation

(2)

dy/dx = x

looks as follows:

Two possible solutions (of infinitely many) to (2) are show on the graph below. One solution passes through the point (0,0), while the other passes through the point (0,-1). Note that, for each point a solution passes through, the corresponding curve below follows the same direction as the line segments around it.

The above direction fields were drawn using a MATLAB toolbox called DFIELD, with which we will presently become familiar. First, let us think about the above direction field without the aid of a computer.

The differential equation (2) can be solved analytically, with solution

(3)

y(x) = x²/ 2+ C

where C is a constant.

Notice how the direction field above confirms that (3) is the solution to the differential equation (2): if we start at any point on the graph, and follow the lines of the direction field, we get a curve of the form given in (3). The constant C is determined by our starting point (i.e., the initial conditions).

Exercise 2.1
	(a) Sketch (by hand, without using MATLAB) the direction field of the following differential equation: dy/dx = -y/5.
	(b) On your direction field, add a curve that approximates the solution passing through the point x = 0, y = 1.
	(c) Now solve the differential equation given in part (a) (either by hand, or you may use the dsolve command as we saw in lab 1), and compare to your answers to parts (a) and (b).

It is apparent from the above exercises that drawing the direction fields associated to more complicated differential equations will be a tedious task. Let us now take a closer look at the MATLAB toolbox which will help us.

§2.2 DFIELD

DFIELD is a toolbox for MATLAB written by John C. Polking that draws direction fields of a given first order ordinary differential equation. The program has already been installed on the ACS machines where most of you will be working. To run it, start MATLAB and at the prompt type

>> dfield7

If you are using MATLAB on a non-ACS machine, you will need to download the necessary m-file from http://math.rice.edu/~dfield to your working directory. (Remember, if you are using MATLAB at home or outside the designated time and place, you are doing so at your own risk! The TAs cannot answer emails or provide help for technical questions arising from students using MATLAB at home). Alternatively, dfield is available for use without MATLAB from http://math.rice.edu/~dfield/dfpp.html, but be warned that the interface is slightly different from the DFIELD used in MATLAB. The course TAs and Tutors will be unable to provide help for it, so we recommend sticking with MATLAB.

The following window should appear.

Let us now try to use DFIELD. Consider the differential equation

(4)

dy/dx = (e^-x - y)(e^-x + 2 + y)

Drawing the direction field for this differential equation by hand would be very tricky, so let us instead let DFIELD do it for us.

Exercise 2.2
	(a) In this exercise, we are going to use DFIELD to plot the direction field of (4). This can be done by entering the following values into the DFIELD setup window. In the boxes under the words "The differential equation" enter y' = (exp(-x) - y)(exp(-x) + 2 + y). (Note that y goes in the left box, and (exp(-x) - y)(exp(-x) + 2 + y) goes in the right box.) In the independent variable box enter x. For now we are going to leave the Parameters & expressions boxes empty and we are going to let the display window section alone. Click on the Proceed button. Now click in a few places on the plot. What are the lines which DFIELD has drawn? Also, using Options -> Keyboard input from the display window, enter the values x=2 and y=3 and click Compute. To include plots in your write-up, use Edit -> Copy Figure and paste the plot into the Word document which you intend to hand in.
	(b) Considering how complicated differential equation (4) appears to be, what do you think is the utility of plotting direction fields?

Remark 2.1 Notice in the previous exercise, you were instructed to use Options -> Keyboard input. Also, you saw that you can always plot solutions by clicking on the graph itself, so what is the advantage to using Options -> Keyboard input? Well, because of the way the pixels are arranged on the screen, it will never be possible to obtain complete accuracy when you simply click on the screen. In other words, you will never be able to click on the exact location of x = 2, y = 3 no matter how hard you try. In some cases, clicking "close" to the desired location will be sufficient, but you'll see in the exercises that follow examples of equations for which it will matter. Therefore, if a problem asks you to find a solution with given initial conditions, get into the habit now of using Options -> Keyboard input.

Remark 2.2 If your DFIELD window becomes too cluttered, you can clear it by selecting Edit -> Erase All Solutions.

Let us now consider another differential equation and see additional ways DFIELD can help us.

§2.3 Initial Data

Differential equations often arise as models for real-life behavior. In this section we are going to see, using direction fields, a difficulty that can result if experimental findings are not accurate enough. First let us suppose that the following equation

(5)

dy/dx = x + 2y

arises as a model in a physics experiment.

Exercise 2.3 Plot the direction field of (5). Suppose that the experiment also reveals that the initial value is about (0,1). As in Exercise 2.2, use Options -> Keyboard input to plot the solution passing through the point (0,1). (That is, use x = 0, y = 1 and click Compute.) Now click on points near (0,1) on the graph itself. Using these plots, think about what would happen if the initial value in the problem were not exactly (0,1)? Would this greatly affect what the solution looks like for this differential equation?

Now suppose that

(6)

dy/dx = y - 1

arises as a model.

Exercise 2.4 Plot the direction field of (6). Suppose that, according to our model, the initial value is about (1,1). Use Options -> Keyboard input to plot the solution passing though (1,1) and then click on a few points near it. If the initial value were not exactly (1,1), how would does this affect the solution?

§2.4 Zooming In

Consider the initial value problem

dy/dx = x² - y³, y(0) = 0

Enter this into DFIELD and use Options -> Keyboard input for the initial condition. This produces the following direction field:

From the above picture, it doesn't appear that the curve follows the direction field. This is because we are not looking closely enough. We can zoom in by going to EDIT -> Zoom in and then either clicking somewhere to zoom in on the direction field or holding the left mouse button and then dragging the cursor to create a rectangle; DFIELD will then zoom in on this region.

In the above example, if we use the cursor to create a box with one corner at (4,2) and the other at (6,4) we obtain:

Upon this closer inspection we can see that the curve does indeed follow the direction field. The lesson here is that, while the DFIELD tool works well, direction fields often have subtle behavior which can require you to be careful when using them.

§2.5 Newton's Law of Cooling

We now look at a real-world application in which direction field plots will give us valuable information about the solutions to our differential equation. Newton's Law of Cooling models the temperature change of an object at a certain temperature when placed in a surrounding environment of a different temperature. The law can be stated as follows:

(7)

dy/dt = k ( A - y )

where y(t) is the temperature of the object in degrees Fahrenheit at time t in hours; and A and k are constants.

Exercise 2.5 Enter the above equation into DFIELD. In the Parameters & expressions section enter A = 1 and k = 2. In the section below that, change the minimum value of t to 0 (since we are not interested in negative values of t here). Plot the direction field and include it in your Word document. Plot the direction fields for different values of A and click on the direction field to plot some solutions for each of these values. What property do you think A represents in real life? (Think about the temperature at which the solutions stabilize.) Note: the independent variable is now t so be sure to make the corresponding change in the DFIELD setup window.

Let us try to figure out how long it will take to defrost a frozen chicken breast in the fridge, which keeps a constant temperature of 41°F. The succulent skinless boneless chicken breast has been in the freezer so its temperature is uniform at -6°F. We'll suppose k = 0.4, based on the properties of the chicken.

Exercise 2.6
	(a) Recall that an initial value problem consists of a differential equation along with an initial condition. Write out the initial value problem which we must solve here. (We already have the differential equation, so this means you need to find the appropriate initial condition.)
	(b) To simulate the conditions in the fridge we must pick the parameter A. What do you think the value of A should be?
	(c) Let us consider the chicken breast fully defrosted when the temperature at its center reaches 39°F. How long does it take to defrost a chicken breast under the above conditions? A rough estimate from a direction field plot is sufficient. (Hint: You may need to adjust the size of the display window, using the minimum value of y and maximum value of y boxes in the setup window. Also you may find using Options -> Keyboard input useful when entering initial data.)
	(d) How much time would be saved if the delicious chicken breast were thawed on the kitchen counter instead, given that room-temperature is around 69°F?

§2.6 Systems of ODEs

So far we have examined differential equations of the form

dx(t)

= f(x(t))

For example, if f(x) = 2x + 3 equation x'(t) = 2x(t) + 3. A system of differential equations involves several equations which tie together one or more variables. For example, a 3x3 system of first order differential equations

dx(t)

= f₁(x(t), y(t), z(t))

dy(t)

= f₂(x(t), y(t), z(t))

dz(t)

= f₃(x(t), y(t), z(t))

constrains the functions (x(t), y(t), z(t)). A more specific example is given by:

(8)

x'(t) = 2x(t) + 3y(t) - 4z(t)

y'(t) = -3x(t) - y(t) +2z(t)

z'(t) = x(t) + y(t) - z(t)

Here f₁(x, y, z) = 2x + 3y - 4z, f₂(x, y, z) = -3x - y + 2z, and f₃(x, y, z) = z + y - z.

Different mathematical models, such as Newton's Laws, produce a system like this and the initial conditions (x(0), y(0), z(0)) = (x_0,y_0,z₀) are determined by whatever data you measure at t = 0. Systems with hundreds of equations or more are commonplace in engineering, biochemistry, and most types of technology. A major conceptual point is that such a system produces a direction field (in this case in three dimensions). Then, solutions can be seen as curves, with very similiar properties to what you saw with DFIELD for 2x2 systems. The course has not yet covered systems of equations, but don't worry! Our goal in these labs is to give you a taste of what a system of differential equations is, a few examples and an idea where the course is going. These labs will discuss everything you need in order to grasp this.

A system with hundreds of equations can be difficult to grasp and impossible to visualize. The way to understand large systems is to get a good understanding of systems of two or three equations, both graphically and algebraically. In systems of hundreds of variables, the algebra is very similiar to the smaller systems, and even though we can no longer view our system graphically, examining solutions visually in lower dimensional systems provide us with an intuition which applies to the larger system. The next few sections will help to familiarize you with the basic definition and ideas of two dimensional systems, which will prepare you for when systems are covered more thoroughly in Assignment 4.

A two dimensional first order system is a pair of differential equations of the form

dx(t)

= f₁(x(t), y(t))

dy(t)

= f₂(x(t), y(t))

with an initial condition (x(0), y(0)) = (x_0,y₀). Under conditions on f₁ and f₂ (i.e. each is continuous and differentiable), specifiying (x_0,y₀) determines the solutions (x(t), y(t)) uniquely - that is, there can be no other functions (x(t), y(t)) which satistify both differential equations and satisfy the initial conditions x(0) = x₀, y(0) = y₀.

An excellent way to think of a system of ODE is in terms of the left hand side giving a tangent vector and the right side giving a vector field. To be more detailed, to each point (x, y) associate a 2-vector

( f₁(x, y), f₂(x, y) ).

This captures all of the information in the right side of the ODE. To understand the left side of the ODE, recall from M20C that a curve C parameterized by (x(t), y(t)) has derivative (x'(t), y'(t)) pointing tangent to C. Put these two facts together to see that a solution curve to the ODEs must at each point (x,y) have the vector ( f₁(x, y), f₂(x, y) ) tangent to it. Different initial conditions (x(0), y(0)) produce different solution curves.

Since a vector has both a magnitude and a direction, it can be represented visually in two dimensions by an arrow. Thus we can "draw" the vector field ( f₁(x, y), f₂(x, y) ) and look at curves tangent to it, which indeed are trajectories of solutions to our system of ODEs. This is much like we did with DFIELD. (Note that DFIELD ignores the magnitude of vectors.) Rather than go into great detail our goal now is modest and all you need to do is look at PPLANE, the companion program to DFIELD. That is the next topic.

§2.7 PPLANE

For a given system of first order differential equations and a given point on the plane, there is one solution to that differential equation which passes through that point. When we draw a vector field, what we are doing is picking a sample of points on the plane and then, for each point, plotting the vector tangent to the solution curve which goes through that point. The resulting plot of vectors then helps us to visualize how different solution curves behave. Now we use the PPLANE program to see some examples.

PPLANE is an add-on for MATLAB written by John C. Polking. It will, given a 2x2 system of ODEs, plot the vector

(f₁(x, y), f₂(x, y)) at the point (x, y).

This vector is tangent to any solution curve through the point (x(t), y(t)), since

(x'(t), y'(t)) = (f₁(x(t), y(t)) , f₂(x(t), y(t))) .

What this means is, we get vector fields, like the direction fields we drew using DFIELD, but there is a big difference - DFIELD gave us the direction field for, say, the t-x plane, where x is the dependent variable and t is the independent variable. With PPLANE our plane is the x-y plane, where x and y are both dependent variables (dependent on t). Like DFIELD, PPLANE will also plot the trajectory C swept out by a solution for each particular initial condition. However, with PPLANE at any particular time t all we know is that we are on C, but we do not know exactly where.

If you are working on an ACS-machine then PPLANE is already installed. To run it, start MATLAB and at the prompt type

>> pplane7

If you are using MATLAB on a non-ACS machine, you will need to download the necessary m-file from the author's website to your working directory. (Remember, if you are using MATLAB at home or outside the designated time and place, you are doing so at your own risk! The TAs cannot answer emails or provide help for technical questions arising from students using MATLAB at home). Like DFIELD, the program PPLANE can be accessed online by visiting http://math.rice.edu/~dfield/dfpp.html.

But again the user interface is a bit different than the MATLAB version.

The following window should now appear:

If you click proceed, the following will appear:

Let us try a simple example. Consider the 2x2 system

(9)

x' = 2

y' = -3

We can solve this system easily by hand:

x = 2t+C₁, y = -3t+C₂, where C₁ and C₂ are arbitrary constants. These are the parametric equations of a straight line.

Exercise 2.7 Use PPLANE to plot the direction field of (9). (To do this you will need to set the values in the x' and y' boxes in the setup window, leaving everything else the same). Click at some points on the diagram to see the solution curve through those points. Now try changing the values of x' and y' from 2 and -3 to other constant values. How does this change the direction field?

Next consider the following system:

(10)

x' = y

y' = -x

This was chosen to be a problem easy enough to be solved by hand, since doing so when you first try a program is extremely valuable in checking that the program behaves a you expect. Let us now compute the trajectory C of the solution: dividing the two equations, we have dy/dx = (dy/dt) / (dx/dt) = -x / y. This gives us a separable ODE,

ydy = -xdx

and integrating we find that solution curves should look like y²+ x² = R, where R is a constant.

Before we plot the direction field, let us modify the viewing window. This is done by changing the x and y minimum and maximum values.

Exercise 2.8 Set the minimum values for both x and y to -10, and the maximum values to 10. Run PPLANE: (10). Click on some points on the diagram to see the solution curve through those points.

There is an important connection between DFIELD and PPLANE:

Exercise 2.9 Plot the solution to the differential equation y' = y^3 + x^2 using DFIELD. Now apply PPLANE to the system

x' =1 x(0)=0

y' = y^3 + x^2

Take a few different initial values for y(0). What do you notice?

We'll explore this idea more in Assignment 4.

§2.8 Predator-Prey Systems

We now consider a mathematical model of population dynamics: a predator and prey system. Suppose that two species of animals, say, foxes and rabbits, coexist on an island. The foxes eat rabbits, and the rabbits eat vegetation, which is available in unlimited amounts. The following system ODE's, known as the Lotka-Volterra model, can be used to model this situation.

(11)

x' =x*(a - b*y)

y' = c*y*(x - d)

where x is the population of rabbits, y is the population of foxes, and a, b, c and d are positive constants (See section. 9.5 in the textbook for more about this). We will now use PPLANE to investigate the solutions to this system.

	Exercise 2.10
	(a) Enter the system (11) above into PPLANE, setting the parameter values a = b = c = d = 1. Plot the vector field for the system. (Remember to use '*' whenever you need to multiply). Set the minimum values of x and y to -1, and the maximum values to 5. Where in the x-y plane are the physically possible solutions (remember - x and y represent populations!)?
	(b) Click on a few different places in the graph to view some solutions. Make sure you click in the first quadrant so that your populations are positive. Note that occasionally due to rounding error, PPLANE will fail to close a loop and will keep drawing loops for an extra minute or so. Don't worry about this occurence. Include the graph in your word document.
	(c) The predator-prey system has the following behaviors: A. The populations of foxes and rabbits are both relatively small. B. The small number of foxes allows the rabbit population to increase C. The increased number of rabbits allows the number of foxes to increase D. The increase in the fox population causes the rabbit population to decrease. E. The decreased supply of rabbits causes the fox population to decrease, returning to behavior A. On your plot in your Word document mark one of the solutions (that is, one loop) with the behaviors A, B, C, D and E where they occur on that solution.

§2.9 Feedback

Introduction to MATLAB

We're almost done, but before we finish we'd like to get some (optional) feedback from you so that we can improve this lab. If you have any comments, write them (by hand if you like) on a separate page at the end of your assignment, so that we can tear it off before we give back your assignment.

Feedback (Optional)
	Do you have any suggestions to help improve this lab? This includes errors or potential errors, typos, unclear statements or questions, or anything else that comes to mind. The more specific you can be, the easier it will be for us to improve the lab.

§2.10 Conclusion

In this lab, we have seen that plots of direction fields and vector fields can be useful when trying to understand the behavior of solutions of ordinary differential equations. This is particularly important when we cannot solve the equation analytically. We have also seen how useful MATLAB’s DFIELD and PPLANE tools can be in generating direction fields and plots of solutions. It is important to understand that these techniques do not in any obvious way scale to systems of hundreds of equations. Chapter 7 of BdeP and Assignment 4 will present techniques which are effective for large linear systems and as such are a mainstay of modern technology and science.

Hopefully, many of you are curious. How are those solution trajectories we just saw produced? The answer is by (approximately) solving ODEs via numerical approximation, that is the topic of Assignmnent 3. Numerical methods typically work well for large systems.

Last Modified:

06/30/2009

2. Visualizing Solutions to ODEs using DFIELD and PPLANE

2. Visualizing Solutions to ODEs
using DFIELD and PPLANE