Appendix 1. Modeling With Differential Equations

§A1.1 What is Modeling?

One of the main reasons that differential equations are so important is that they often describe physical phenomena. Certainly, there are many mathematicians who study differential equations from a more theoretical point of view and who solve certain kinds of equations with no known applications to the "real world," but the vast progress made in the last several hundred years of studying differential equations has come from attempts to understand the physical world. Applications of differential equations are present in almost any scientific field of study imaginable. Try to think of an obscure or uncommon realm of scientific inquiry--let's say "underwater acoustics"-- and type "underwater acoustics, differential equations" into any web search engine. You'll probably get at least a few hits. No matter what subject you choose, you can find some connection to differential equations.

Modeling is the process of creating an equation or system of equations that describes or predicts--more or less accurately--some physical situation. There are countless possibilities. We may have gathered some data through observation or experiment and wish to find an equation that matches our observation. We may have already hypothesized a particular differential equation and then use it to predict the future behavior of a physical system. We may wish to vary the parameters of an existing differential equation to save the time and money of having to conduct multiple or lengthy experiments.

§A1.2 The Effectiveness of Modeling

Mathematical models have proven their worth time and time again, and it's no wonder that scientists and mathematicians alike continue to make differential equations one of their primary research tools. Nevertheless, there are a few considerations which need to be observed if the model is to be effective.

It is fair to say that no model can be completely perfect. First of all, certain theoretical assumptions have to be made and these assumptions may have limited accuracy and applicability. In a physics problem, for example, we may have to assume that acceleration due to gravity is 9.8 m/s². This may be enough under certain conditions, but we may obtain incorrect results for events taking place far from the Earth's surface, or on a different planet for that matter. Also, we know that even in the most rigorously controlled environments, there are many subtle factors that can affect the outcome of an experiment. It is very difficult to account for all such possibilities in one model, and even if we could, our equation would likely be so complicated that it would present difficulties when we attempted to solve it.

On the other hand, if we rely too much on simplicity in our model, we may end up with an equation that is powerless to predict future behavior. Or perhaps the data would have to be subject to an extremely poor margin of error to fit an oversimplified model.

Therefore, the creation of a model can be a delicate balancing act of maintaining enough simplicity to be of computational effectiveness, but enough complexity to reflect the subtle dynamics of our physical system. We'll explore such issues as we work through the creation of a basic population model.

§A1.4 Population Model: First Approach

There are a few ways to approach a population model. One common way is to make observations over the course of some period of time, record those observations, and then try to develop an equation that fits. While this is a useful technique in the real world, it won't be of much benefit to us in a computer lab. Instead, we'll make certain assumptions about the type of growth we expect, and proceed from there. However, throughout the lab we will be checking how effective our models are by comparing them to the actual data reflecting the world population from 1950 to the projected population in 2050.

Let us try to come up with a reasonable model for the human population on earth. As with any model, one of the first things to do is to come up with variables. In the case of population, we are interested in population as a function of time. So let t denote time and y(t) the population at time t. Then, of course, the derivative y'(t) will denote the rate of change of population at time t.

Now that we have the notation in place, we need to think about the actual model. In other words, can we think of a reasonable law that will relate the rate of change of population with the actual population at that time? It may seem plausible that the more people we have at a certain time, the faster the population will grow. In other words, the population growth at time t is proportional to the population at time t. Thus we obtain our first model:

y'(t) = r*y(t)

where r is some proportionality constant.

Now let us see if this model can withstand the "test of time". The following data was obtained from the United Nations website and reflects the world population from 1950 to 2050 in billions. The table was last updated in 2005, so that the values after 2005 are projected. (They were probably calculated using more data and more sophisticated mathematical models than the ones we will use in this lab). However, we will treat the entire table as our experimental data.

1950	2.54
1955	2.77
1960	3.03
1965	3.34
1970	3.70
1975	4.08
1980	4.45
1985	4.86
1990	5.29
1995	5.72
2000	6.12
2005	6.51
2010	6.90
2015	7.30
2020	7.67
2025	8.01
2030	8.32
2035	8.59
2040	8.82
2045	9.03
2050	9.20

Exercise A1.1

(a) The equation y'(t) = r*y(t) is separable. (It will be easier to see if you rewrite it as dy/dt = r*y.) Verify that the general solution is y(t) = A*exp(rt).

(b) Now suppose t measures time in intervals of 5 years, and t = 0 corresponds to the year 2000, i.e. when t = 2 we are 10 years from 2000, that is, we are at 2010. And when t = -3, we are 15 years before the year 2000, that is in the year 1985. Thus y(0) = 6.12 and y(1) = 6.51. We can use this data to find the values of A and r. To do that we observe that y(0) = A = 6.12, and y(1) = 6.51 = A*exp(r) = 6.12*exp(r).

Thus r = log (6.51 / 6.12). Use MATLAB to evaluate this expression. Simply type in:

>> r = log(6.51/6.12)

Include the input and output in your Word document.

Before we move on, let us define a function P which will give us the population at time t, according to our model. This entails entering the following commands into MATLAB.

Example A1.1

First define what the constant A is.

>> A = 6.12;

Now define our function P as a string :

Remark A1.1 A string is just the name of a data type that stores collection (string) of characters. The contents of a string type are usually delimited by a set of quotation marks. Sometimes it is important to distinguish a string type from a numeric type because certain functions can only be applied to numbers and cannot be applied to strings. Therefore, it is often necessary to convert from string type to numeric type or vice-versa to be able to take advantage of certain functions. Don't worry too much about these technical details.

>> P = 'A*exp(r*t)'

Next we need to substitute the numerical values of A and r into our string object.

>> P = subs(P)

This command returns a symbolic object, so we need to convert it back to a string and define our inline function P.

>> P = inline(char(P))

At this point our function P is defined and available for us to use.

Let us see how close our model is to the data in the U.N. table. We can find the population in 2015 using the results we obtained in the previous exercise. Try typing in the following:

>> P(3) ans = 7.3661

This gives us the population when t = 3, or in the year 2015, according to our model. It is interesting to see how close it is to the projected population in 2015. Let us compute the error by typing in the following:

>> 100*(P(3) - 7.30)/7.30

ans = 0.90606

The error is less than one percent! That's pretty good for a simple model!

Exercise A1.2
	(a) Just as in the example above, use our model to find the population in years 1990, 2005, 2020. Compute the error relative to the values in the table.
	(b) Now do the same to predict the population in years 2035, 2050. Compute the error relative to the values in the table. What can you conclude about the effectiveness of our model? When does it seem to work well, and when does it seem to break up? Include the input and the output in your Word document.

As we saw, our model is far from perfect. Let us analyze this a bit further in the next exercise.

Exercise A1.3

(a) Let's plot the solution we obtained in the previous exercise for 0 ≤ t ≤ 50. Remember that when t = 50 we are 250 years from 2000. You'll need to enter the following commands:

>> fplot (P, [0, 50]) >> grid on

From the graph, estimate the population for t = 50. Copy the graph into your write up.

(b) The kind of growth exhibited by the solutions to this differential equation is called exponential growth. As you may have seen before, things that grow exponentially do so very rapidly. Explain why our present model may not be the best model for the world's population in the long run.

§A1.5 Population Model: Second Approach

As we have seen our previous model was rather naive. Though it seemed to work well in predicting population for smaller values of t it failed badly as t grew larger. You may have come up with some natural explanations in the previous exercise as to why population could not continue to grow exponentially in the long run. Let us now try to modify our model to reflect that.

The model we will be using now is due to a Belgian mathematician P. F. Verhulst. In 1846 he proposed that the relative growth rate of the population is proportional to the difference between its limiting value and the actual population. This relationship can be expressed as the following differential equation:

where M is the limiting value of the population. (The three expressions for the right-hand side are all the same algebraically, but for the text and exercises that follow, sometimes it will be helpful to look at them in one form, and sometimes in another.)

Note that the first order term in this differential equation is just ky, and the second order term is (k/M)y². We see that the first order term is exactly the same as the one in our first model (it used to be called ry). So, when the values of y are close to zero, the value of y² will be smaller than y, and in fact will be rather negligible as compared to the linear term ky. Thus, for small values of y, this model will be quite close to our initial model.

In the next exercise let's take some time to familiarize ourselves with what solutions to this equation look like and how they behave in the long run.

Exercise A1.4
	(a) What happens to the rate of change if we use values of y that are very close to M? What if we use very large values of y, i.e., values of y much larger than M? Will the population be increasing, decreasing, or staying about the same in each of these cases?
	(b) Enter the Verhulst equation into DFIELD. (Be careful about placing parentheses and , the multiplication sign.) Leave the independent variable as t. Assign the values k = 1 and M = 10 as parameters and adjust the display window to the ranges 0 < t <* 10 and 0 < y < 15. Assume that we start with 2 billion people and graph the corresponding solution curve. Do the same with 8, 10 and 12 billion people. (For these cases, you will need to enter the initial conditions manually. Remember you will need to go to Options -> Keyboard Input, and then enter your initial conditions.) Copy the display window figure to your Word document.
	(c) What happens to the population in the long run if we start with fewer than 10 billion people? More than 10 billion? Exactly 10 billion? The constant M is called the carrying capacity of the population. Explain why you think it has such a name.
	(d) Recall that in Exercise A1.3(b) you explained why our previous model was insufficient to model population. Does the Verhulst model fix the problem you identified with our previous model? Why or why not?

Now that we have an idea about the qualitative behavior of the solutions to the Verhulst equation, let us see if it is in sync with the data we have from the U.N. table. As before, t will measure time in intervals of 5 years with t = 0 corresponding to the year 2000, and y will measure population in billions. Remember that we need to determine the constants k and M, and the data from the table will help us do just that.

Exercise A1.5
	(a) Solve the equation y' = (k/M)y(M - y). What is the limit of y as the time goes to infinity? Hint: This is a separable equation, which you may already know how to solve (the solution involves partial fractions). However, to save time you may use the dsolve command. If you choose to use the dsolve command, don't forget to include your input and output in your final write-up.
	(b) It can be shown that: where y₀= y(0), y₁= y(1) and y₂= y(2). (You do not have to work it out, since the algebra may get a little long). Let us just compute the numerical value of M. We will need to enter the following commands in MATLAB: `>> y0 = 6.12` `>> y1 = 6.51` `>> y2 = 6.90` `>> M = (2y0y1y2-y1^2y0-y1^2y2)/(y2y0-y1^2)` Include your MATLAB commands along with the output in your write up. According to this model, what is the carrying capacity of the world population?
	(c) From part (a), note that 6.12 = y(0) = M / (1 + B), where B is a constant. (Note that if you used dsolve, we have B = C1M). Find the value of B. Also show that and evaluate it numerically using MATLAB. (Hint:* Plug in t = 1 into the solution in part (a) and solve for k) Include your input and output in the final write up. For future use, make sure these constants are stored in MATLAB.

At this point we can test our model against the actual data from the UN table. This will be the goal of the next exercise.

Exercise A1.6
	(a) Just as we did before for our first model, let us now define the function V which will tell us the population at time t according to the Verhulst model. First, make sure your constants B, M, and k are already defined in MATLAB from the above exercise! Now enter the following commands: `>> V = 'M/(1 + Bexp(-kt))'` `>> V = subs(V) >> V = inline(char(V))` Don't be alarmed if your output looks complicated. It's supposed to. Please don't forget to include your input and output in the write up.
	(b) Using the inline function V you defined, compute the population in the years 1985, 2005, 2025. Compute the error relative to the values in the table. Copy and paste your commands along with the output in your Word document. Which model seems to give better results for these years?
	(c) Repeat the previous question for years 2040, 2050.

In the last exercise, we will try to get a visual picture of how the two models differ.

Exercise A1.6 To compare our two models let us plot the solutions we found in section A1.4 and A1.5 on the same set of axes, with -10 ≤ t ≤ 20.

>> fplot (P, [-10, 20]) >> hold on >> fplot (V, [-10, 20], 'r') >> grid on

Include your input and the resulting figure in your lab write up. You should notice that the results from using the two models are very close when the population is small, but then seem to diverge from each other as the population grows bigger. Can you explain, mathematically, what accounts for such closeness of the two solutions for small values of the population? (Hint: Compare the differential equations for the two models. How do their first order terms compare?)

§A1.6 Conclusion

Of course, population models only represent a very small subset of all the types of models that are out there. Hopefully, though, this assignment gave you an idea of some of the principles that go into the creation of a model. We've also explored some of the difficulties of producing a good model that accurately reflects some physical process. As you encounter differential equations in whatever field of study you pursue, you should always assess what assumptions have to be made for the model in question to be valid.

Last Modified:

07/02/2008