Matlab m files download slopefield

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

using direction fields (also called slope fields). The direction field of this differential equation is a diagram in the plane in which there is a small line segment drawn with slope at the point . For example, the direction field of the differential equation

looks as follows:

Two possible solutions (of infinitely many) to (2) are shown 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 with the help of a MATLAB M-File called , with which we will shortly become familiar. But first, let's think about the above direction field without the aid of a computer.

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

where 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 we follow the lines of the direction field, we get a curve of the form given in (3). The constant is determined by our starting point (i.e., the initial conditions).

1. Sketch (by hand, without using MATLAB) the direction field of the differential equation for and values between -5 and 5. (You do not need to include this sketch in your write-up.)
2. On your direction field, add a curve (by hand) that approximates the solution passing through the point .
3. Now solve the differential equation given in part (a), either working it out by hand or using the command that we saw in Assignment 1. Compare 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 some convenient tools that can help us.

Drawing Direction Fields & Solutions

Before we begin, you will need to make sure you have the four M-Files , , , and installed. If you are working on an ETS machine or using the virtual lab, they should already be there. If not, you can download them all here. Once you've downloaded this folder, move each of the four M-Files into your working directory, which can be located at your current folder. If you can see those files in your current folder, then it's safe to use them. If you choose to use MATLAB online or MATLAB Mobile, you may upload those four M-files to your MATLAB Drive.

SLOPEFIELD

The M-File defines a new function, , which draws direction fields of a given first-order ordinary differential equation. If you like, you can type into MATLAB to read the documentation for the function. We're going to work through an example to see how works.

Let's try to use to draw a direction field for the differential equation

(4)

If we check the documentation, we'll see that takes four inputs: a function, a range of values for the independent variable, a range of values for the dependent variable, and finally a number () that controls how many segments we'll draw on our field.

Thus, to use , we first need to define a new function (say, ; the name doesn't matter) that represents the right-hand side of this ODE. Recall that we can do that by typing

Suppose we'd like to see the direction field for values between -4 and 4 and values between -3 and 5, with about 20 marks in each direction. Then we simply need to type

MATLAB produces the following plot:

For the remainder of this lab, we will refer to the independent variable in our differential equations as time, although the name of that variable need not be , and the independent variable for the equations in question may not literally be time.

Now we can use MATLAB to quickly produce direction fields. Let's explore how these direction fields relate to the solutions of the ODE by learning how to plot some solutions.

DRAWODE

MATLAB includes a number of functions designed to numerically approximate solutions to ODEs. However, they are somewhat inconvenient for our purposes because they typically only compute either forwards or backwards in time, but not both. Thus, for instance, you can't call on these functions to plot the solution to for between -4 and 4 with the initial value , because is neither the start time nor the end time.

The M-File rectifies this issue. It defines a new function that can plot solutions to a first-order initial value problem regardless of where the "initial" time is located. In the example we just mentioned, it would be able to go backwards from to and also forwards from to . Let's go ahead and use to do that. Again, it may be useful to check the documentation by entering .

The function takes four inputs: the function that represents (which, in the case of equation (4), we already named ), a range of times, an initial time , and finally an initial value .

Again, we want to be between -4 and 4, and now we'd like to be So we plug all that in to the appropriate places in :

Now MATLAB draws this plot:

Why do you think the graph ends around , even though we asked it to show values down to -4? A look at the direction field from the first example may help to explain. In fact, it would be very nice if we could look at both the direction field and our solution curve at the same time.

Combining Plots

Unfortunately, if we try to run and in succession, the first plot (our direction field) disappears; MATLAB overwrites it with the new graph made by . The keys to solving this are the commands and .

If you use MATLAB to produce a figure and then enter , MATLAB will start drawing any new plots on the same figure, without erasing whatever was there before. After you're done, turns off this behavior so that you can erase your work and draw new plots. Let's see this in action using the work we've already done. (Here, your up-arrow key may prove useful.)

We'd like to draw a direction field for the ODE given in (4) and then plot the solution passing through (1,-1) on that direction field. In MATLAB's command prompt, we'll enter the following:

Finally, we get our combined plot:

We're not limited to graphing just one solution curve, however. Let's try a few more initial values:

You can copy and paste this code into MATLAB to try it for yourself.

Adjusting the Viewing Area

In Assignment 1, we learned that we could adjust the axes on a plot by using the menu option Edit &#; Axes Properties&#;. However, this will not work with , since only draws the direction field over a limited area. To change the boundaries on the axes, we'll have to change the bounds we include in .

We could go through and adjust the range of values every single time we call and , but when you're drawing a lot of solution curves the way we did in Example , that can be a tedious task. Instead, if we had anticipated that we might want to zoom in or out, it would have been better from the beginning to do the following:

Now, if we simply change the value of to , change the value of to 10, and adjust the values in as desired, we can quickly rescale the plot with all six solution curves included.

If you think you'll enjoy the challenge, you might want to devise an even more efficient way to draw a lot of solution curves, perhaps using loops and matrices.

When looking at graphs like these, you might find it useful to try the commands and , which turn on and off the grid lines on the plot, respectively.

Consider the differential equation

(5)

1. Plot a direction field for (5) for and between and Use the commands and to plot at least two solution curves on this direction field, one of which passes through the point (2,3) (that is, and ). Paste your plot into your Word document.

   [Hint: You may manually set the axis limits before calling drawode by where sets scaling for the x- and y-axes on the current plot.]

2. Considering how complicated differential equation (5) appears to be, why do you think we might want to plot a direction field?

Suppose the differential equation

arises as a model in a physics experiment.

Plot a direction field for (6) with and between -5 and 5. On this direction field, plot the solution curve passing through (0,^-1⁄₄). Then draw three more solution curves passing through points very close to (0,^-1⁄₄) on the same figure. Paste your plot into your Word document.

Using this plot, think about what would happen if the initial value in the problem were not exactly (0,^-1⁄₄). Would this greatly affect the solution of the differential equation?

Now suppose that

arises as a model.

Plot a direction field for (7) with and between -6 and 6. Suppose our experiment reveals that the initial value is about (1,1). On your direction field, plot the solution curve passing through (1,1), and also plot several solution curves going through other points near (1,1). Paste your plot into your Word document.

If the initial value were not exactly (1,1), how would this affect the solution?

Modeling: Newton's Law of Cooling

We'll 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:

where is the temperature of the object in degrees Fahrenheit at time in hours, and and are constants.

Plot a direction field for (8) with , , and where the minimum value of is zero (since we are not interested in negative times here). You can choose an appropriate maximum value for and minimum and maximum values for . Include the direction field in your Word document. Now plot some direction fields for other values of , and graph a few solution curves on those direction fields using by choosing some initial values. What property do you think represents in real life? [Hint: Think about the temperature at which the solutions stabilize.]

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 chicken breast has been in the freezer for a while, so its temperature is uniform at -6°F. We'll suppose , based on the properties of the chicken.

1. Recall that an initial value problem consists of a differential equation along with an initial condition. Write out the initial value problem that we must solve here. (We already have the differential equation, so this means you need to find the appropriate initial condition.)
2. To simulate the conditions in the fridge, we must pick the parameter . What do you think the value of should be?
3. Let us consider the chicken breast fully defrosted when the temperature 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 axes on your plots.]
4. How much time would be saved if the chicken breast were thawed on the kitchen counter instead, given that room temperature is around 69°F?

Systems of ODEs

So far we have examined differential equations of the form

For example, if , then this equation says . A system of differential equations involves several equations that tie together one or more variables. For example, a 3×3 system of first-order differential equations

constrains the functions . A specific example is given by:

Here, , , and .

Mathematical models like Newton's laws produce systems like these, and the initial conditions are determined by whatever data you measure at . Systems with hundreds of equations or more are commonplace in engineering, biochemistry, and most types of technology. One major conceptual point is that such a system produces a direction field (in this case in three dimensions). Then the solutions can be seen as curves, with very similar properties to what you saw with for 2×2 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 understand 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 similar to the smaller systems, and even though we can no longer view our system graphically, examining solutions visually in lower dimensional systems provides us with an intuition that applies to the larger system. The next few sections will help to familiarize you with the basic definitions 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

with an initial condition . Under conditions on and (namely, that each is continuous and differentiable), specifying determines the solutions uniquely—that is, there can be no other functions that satisfy both differential equations and also satisfy the initial conditions .

An excellent way to think of a system of ODEs is in terms of the left-hand side giving a tangent vector and the right side giving a vector field. To be more precise, to each point , associate a 2-vector

This captures all the information in the right side of the ODE. To understand the left side of the ODE, recall from Math 20C that a curve parametrized by has derivative pointing tangent to . Putting these two facts together, we can deduce that a solution curve to the ODE must at each point have the vector tangent to it. Different initial conditions 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 and look at curves tangent to it, which indeed are trajectories of solutions to our system of ODEs. This is much like what we did with , although ignored the magnitude of vectors. Our next objective is to put these ideas into action.

Drawing Phase Portraits

Double-check that your working directory contains the M-Files and . If you are working on an ETS machine, they should already be installed; otherwise, you can download them using this link.

For a given system of first-order differential equations and a given point on the plane, there is at most one solution to that differential equation that passes through that point. When we draw a phase plane (or phase portrait), 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 that goes through that point. (This is analogous to the one-dimensional phase diagrams that you may have already seen in lecture; in those diagrams, we only drew our arrows up and down.) For our two-dimensional case, the resulting plot of vectors helps us to visualize how different solution curves behave. We'll use the command to see some examples.

PHASEPLANE

The function plots the phase plane at time zero for a two-equation system of first-order ODEs , where is the independent variable and is a 2×1 vector containing our dependent variables. In MATLAB, we write such a vector as . Let's see how to use with an example.

We'll use to draw a phase plane for the system

(10)

We can think of a system like this in terms of vectors. Let be the 2×1 vector with entries and (our two dependent variables). Then we can write this system as

where is a function that spits out the 2×1 vector with entries and (from the right-hand side of our system).

The syntax for is very much like that for . The command takes four inputs: the vector-valued function that defines the right-hand side; a range of values for the first dependent variable ; a range of values for the second dependent variable ; and finally the number of vectors we want to draw in each direction on our field.

We need to define a vector-valued function, which we haven't yet learned how to do. The procedure is much like the one for real-valued functions that we've used previously. In this case, we want to enter

Here, represents the second element of the vector variable ; for this system, that's our second dependent variable, . Of course, is the first element of and stands in for . You can do this kind of thing in general in MATLAB: given a vector , then returns the i^th element.

Suppose we want to draw the phase plane for our system with between -5 and 5 and between -6 and 4, with a 15×15 grid of arrows. Then we must enter the following:

MATLAB will draw the plot shown below.

We'd now like to see how individual solutions to our system of ODEs appear on our phase plot. We can't graph the solutions themselves, but we can graph their projections onto the plane (their phase paths). Let's see how to do that now.

DRAWPHASE

The command defined by plots the projections of solutions onto our vector field, helping us to see exactly how the solutions behave. Let's use it on our system (10) of ODEs.

The syntax for is as follows:

In this command, is the same function we defined earlier. We always start at time , but specifies the time at which we'd like to stop. Lastly, and give initial values for and , respectively, at time .

Suppose we want to plot the solution curve that starts at , from time zero to time But remember: we want to draw this on the phase plane we've already made, without erasing the previous plot. Therefore, we need to use . We should enter that like so:

Or, if we were starting from scratch, we would enter

The starting point for our curve is marked with an asterisk (), and the ending point (at whichever value we put in for ) is marked with a circle.

If we want, we can plot several more phase paths, adjust the viewing area, and define a variable so that it's easy to adjust our final time, just as we did before with :

Again, it may be helpful to copy and paste this code, rather than re-typing all of it yourself.

It's time for you to draw some phase portraits on your own. Consider the 2×2 system

Notice that in this problem, our dependent variables are named and rather than and . That's okay! The names of the variables don't actually matter. Like the function in Example , we'll have to define a new function in MATLAB to represent the system. This time, will stand for , since that's our first dependent variable; will be .

Before we try to use MATLAB to solve this system, notice that we can solve this system easily by hand: we get and . These are parametric equations of a straight line.

Use the techniques we've learned involving to plot a phase portrait of (11), where the and values are between -5 and 5. Then, on the same plot, use to draw at least three different solution curves. Include the resulting plot in your Word document.

Now, try changing the values of and from 2 and -3 to other constant values. Describe in your Word document how this changes the phase portrait.

Next, let's consider the 2×2 system

This problem is again fairly easy to solve by hand. Whenever we first try out a program, it's a good idea to try some problems for which you already know the answer; that way you can check that the program behaves the way you expect. In this case, by dividing the two equations, we get . Now we have a separable ODE:

When we integrate both sides, we find that the solution curves should look like , where is a constant. What kind of curve is this?

Use to plot a phase portrait of (12), where the and values are between and [Remember that when you're defining your function , you need to use to represent your first dependent variable and for the second one.] Then, on the same plot, use to draw at least three different solution curves. Include the resulting plot in your Word document.

There is an important connection between direction fields and phase portraits.

Use to draw a direction field for the differential equation

and then draw some solution curves on the resulting direction field using . Include your figure in your Word document.

Now consider the system

(13)

+ y3

where .

Tell MATLAB to create a new figure using the command. Then use to draw a phase portrait for the system (13), and plot a phase path on your diagram using ; the value represents and should therefore be zero, while the value can be any initial value for of your choosing. Finally, try adding a few more phase paths using other values for . Paste the resulting figure into your Word document.

What's the relationship between your phase portrait figures and your direction field figure? [Hint: When you drew the solutions on your direction field, you should have picked several different starting values for and . Try drawing some phase paths on your phase portrait with the same initial values. Do the resulting curves match up?]

We'll explore this idea more in Assignment 4.

Modeling: Predator–Prey Systems

We'll 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 we'll say is available in unlimited amounts. The following system of ODEs, which is called the Lotka–Volterra model, can be used to model this situation.

Here, is the population of rabbits, is the population of foxes, and , , , and are positive constants. We will now use and to investigate the solutions to this system.

1. Use MATLAB to produce a phase plane of the system (14) above, with the parameter values . Set the minimum values of and to -1 and the maximum values to 5. Where in the plane are the physically possible solutions? (Remember, and represent populations.)
2. On your phase plane, use to plot three solution curves in the first quadrant from time zero to time Note that because of rounding error, might fail to close a loop exactly. Include the graph in your Word document.
3. The predator–prey system has the following states:
   1. The populations of foxes and rabbits are both relatively small.
   2. The small number of foxes allows the rabbit population to increase.
   3. The increased number of rabbits allows the number of foxes to increase.
   4. The increase in the fox population causes the rabbit population to decrease.
   5. The decreased supply of rabbits causes the fox population to decrease, returning to state A.
   On your plot in your Word document, mark one of the solutions (that is, one loop) with the states A, B, C, D, and E at the locations where they occur on that solution.

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 some relatively simple MATLAB M-Files 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. 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? Behind the scenes, the M-Files we used in this lab work by (approximately) solving ODEs via numerical approximation. You can actually open the M-Files and read them yourself, if you want to see what they're doing. These numerical methods will be the topic of Assignment 3. They typically work well for large systems.