Chapter 8 Modeling with ODEs

In this section we investigate a key area of interest for ordinary differential equations – modeling. ODEs are often used to model physical phenomena as they allow us to express a relationship between an unknown quantity, $y$, and the rate of change of the unknown quantity $y'$. In this section we will model with autonomous ODEs. Many physical and biological phenomena can be modeled using autonomous ODE, and the special structure of autonomous ODEs allow us to extract certain types of solutions, called equilibrium solutions, with little effort.

8.1 Exponential Models

8.1.1 Exponential Population Growth

Consider a population of bacteria living a Petri dish. Let $P(t)$ be the population of bacteria in the dish at time $t$. Assume that $t$ is measured in hours. When using differential equations to model a physical phenomenon, we often begin by asking if we can express the rate of change of the quantity in a simple way. That is, we are looking for an expression for $P'(t)$.

One of the most common models for modeling a population is to assume that the rate of change on the population is proportional to the current population. In this context, this means that if there are only two bacteria present at time $t$, the population can’t be growing too quickly, since there are only bacteria present to divide. However, if there are 1000 bacteria present, then there are many cells present to divide, and so will be quickly adding new bacteria to the population.

If the rate of change of the population, $P'(t)$ is proportional to the current population, $P(t)$, then we have $$P'(t)=kP(t)$$ where $k$ is a growth parameter. If $k$ is a large number, the population grows quickly. If $k$ is a small, positive, number, then the population grows slowly. If $k$ is a negative number, then the population decreases, since then we would have $P'(t)<0$. The dimension of $k$ are $1/T$. This means that the dimension $1/k$ is $T$. It turns out that $1/k$ is the average amount of time a cell spends “resting” between divisions.

8.1.1.1 Example

Suppose that a bacterial colony initially has 100 members, that the rate of change of the population is proportional to the current population, and that on average a cell spends 5 hours resting between cell divisions. Find a formula for the population of the colony at time $t$.

Let $P(t)$ be the population at time $t$, where $t$ is in hours. Since the rate of change of the population is proportional to the current population, $$P'=kP.$$ Since the average amount of time between cell divisions is $\frac{1}{k}=5$ hours, we have $k=\frac{1}{5}$ hour$^{-1}$. We have $$P'=\frac{1}{5}P.$$ This equation is separable. Our standard method for solving separable differential equations (assuming $P>0$) gives us an answer. $$\begin{aligned} &\frac{dP}{dt}=\frac15 P\\ &\frac{1}{P}\frac{dP}{dt}=\frac{1}{5}\\ &\int \frac{1}{P}\frac{dP}{dt} dt = \int \frac15 dt\\ &\int \frac{1}{P}dP = \int \frac15 dt\\ &\ln(P) = \frac{1}{5} t + C_1\\ &P=e^{1/5 t+C_1} = Ce^{1/5 t}. \end{aligned}$$

We can use our initial condition to determine $C$. $$100=P(0)=Ce^{1/5\cdot 0} = C.$$ Finally, we have our specific solution to the IVP. $$P=100e^{1/5 t}.$$ You might recall this function from Block 1. This model is called exponential growth.

8.2 Radioactive Decay

Another example of modeling with a 1$^st$ order, autonomous ODE is the decay of a radioactive element. Atoms that are radioactive spontaneously decay into atoms of a different element. Let $y(t)$ be the mass of our radioactive element present at time t.

The most common assumption regarding the rate of change of $y$ is that the rate of change of the amount of element present is proportional to the amount of the element present. Where there are many atoms of our radioactive element present, there are many atoms available to transmute, and thus we can lose a lot of radioactive substance quickly. When there is little of our element present, there are few atoms available to transmute, so we can’t lose material quickly.

As with population growth, if the rate of change of the amount present is proportional to the amount present, then $$\frac{dy}{dt} = -k y.$$

Once again, if $k$ is large then the element decays quickly. If $k$ is positive but small, then $k$ decays slowly. If $k<0$ then the amount of element present actually grows, since we would have $\frac{dy}{dt}>0$.

As with population models, the dimension of $k$ is time$^{-1}$. Therefore $\frac{1}{k}$ has dimension time.

Commonly we refer to the half-life of an isotope which can be found directly. We want to solve for $t$ when the number of atoms is half of the initial value. Solving for an unknown exponent is a common application of logarithms. Here we find the half-life, $t_{h}$, by applying the definition. At the half-life, $t_{h}$, the number of atoms remaining equals half of the original: $N\left( t_{h} \right) = \frac{N_{0}}{2}$ .

Half-life is found as

$${N\left( t_{h} \right) = \frac{N_{0}}{2} = N_{0}e^{- \lambda t_{h}} }{\frac{1}{2} = e^{- \lambda t_{h}} }{\ln(1/2) = - \lambda t_{h} }{t_{h} = - \frac{\ln(1/2)}{\lambda}}$$

We see that the half-life of the isotope depends directly on the decay rate.

Exponential decay is often a good approximation to physical phenomena. Radioactive isotopes decaying, medication metabolizing, and resale value of automobiles are examples where exponential decay provides a respectable model. However, exponential growth is different. Physical systems exhibit exponential growth for only a limited time. Consider the example with rabbits. If the rabbit population actually grew exponentially, then our planet, our solar system, our galaxy, … would be buried in an ever-increasing pile of rabbits. To improve the modeling approach for these types of dynamics, we use the logistic equation.

8.3 Logistic Equation

The logistic equation expands our model for the rate of change. Previously we introduced a simple model where the rate of change is modeled by a line, e.g. $P^{'} = kP$. Here the rate of change has the form $kP + c$ where $k$ is the slope of the line and the constant $c = 0$.

The logistic equation uses a quadratic model for the rate of change. For example, we can represent a population rate as

$$P^{'} = kP - qP^{2}$$

where $k > 0$ and $q > 0$ are parameters of our model. As we saw previously, this differential equation is autonomous because the rate of change of the population, $P^{'}$, only depends explicitly on the population, $P$. Moreover, for populations that are “large enough” the term $qP^{2}$ will be greater than $kP$. In this situation, the rate of change, $P^{'}$, will be negative, and the population will be decreasing. Accurately modeling a population focuses on selecting the parameters $k$ and $q$.

Consider the logistic equation written in a factored form:

$$P^{'} = kP\left( 1 - \frac{q}{k}P \right)$$

Now we recall our previous discussion about equilibrium solutions. At what values of the population will the rate of change equal 0? To answer this question, we solve $P^{'} = 0$:

$$0 = kP\left( 1 - \frac{q}{k}P \right)$$

$$kP = 0\ \ \ \ or\ \ \ \ 1 - \frac{q}{k}P = 0$$

Solving, we find two different equilibrium solutions. Let us consider both.

1)	$$P(t) = 0$$	This makes sense. If there are no organisms, then there will be no rate of change.
2)	$$P(t) = \frac{k}{q}$$	This population is more interesting. At the $k/q$ population, births and deaths are occurring at the same rate. The $k/q$ population is called the carrying capacity.

With a seemingly minor adjustment to our model, the concept of a maximum sustainable population is introduced through carrying capacity. The logistic equation is a powerful step forward in modeling.

Example: Consider a population of squirrels, $S(t)$, with a growth rate of $k = 0.3\ week^{- 1}$. The squirrels live in a park that can support a maximum of 225 squirrels. What is the IVP that models this situation?

Start by considering that the squirrels increase at a rate proportional to their population, and the proportionality constant is the growth rate:

$$S^{'} = 0.3S$$

This initial model is updated to a logistic equation by including the carrying capacity of the park:

$$S^{'} = 0.3S\left( 1 - \frac{S}{225} \right)$$

The initial value problem combines the model with some initial number of squirrels, $S(0) = S_{0}$.

Using the graphical analysis discussed earlier, we consider the direction field and phase line for this model. These graphs provide rate information as well as the values of the equilibrium populations.

We can see that an initial population of 40 squirrels, i.e. $S_{0} = 40$, leads to the population approaching 225 in approximately 14 weeks.

Continuing with this example, consider that a hawk has moved into the park where the squirrels are living. The hawk is an efficient hunter, and he eats 2 squirrels each day. How can we adjust our model to account for this change? We must keep in mind that $S^{'}$ is the rate of change of squirrels, and the units on the growth rate, $k = 0.3\ week^{- 1}$, give us the time scale for our model.

The hawk provides a weekly change to the squirrel population of $- 14$. This change is independent of the squirrel population; he kills and eats 14 squirrels per week regardless of the population. Thus, the hawk’s impact on the population is: $S^{'} = - 14$ … the rate of change of squirrels is a constant negative 14 each week. We can incorporate this new development into our model:
$$S^{'} = 0.3S\left( 1 - \frac{S}{225} \right) - 14$$

Again, we turn to our graphical analysis to understand the impact of the hawk. We will not develop analytic solutions, but we will understand the dynamics of the squirrel population. This first graph is scaled to match the original, pre-hawk, squirrel model. We notice that there are still two equilibrium solutions, but those solutions have moved toward the middle population (~112 squirrels). We also notice that the solutions away from the equilibria flow in the same direction as they did previously.

The second graph provides a closer look at the population dynamics near the equilibrium solutions. To start our analysis, we solve for those equilibria:

$${0.3S\left( 1 - \frac{S}{225} \right) - 14 = 0 }{- \frac{S^{2}}{750} + \frac{3S}{10} - 14 = 0}$$

This quadratic can be solved directly:

$$S = \frac{225 \pm \sqrt{345}}{2}$$

$$S \approx 66.1\ \ \text{ and } \ \ S \approx 158.9$$

The equilibrium solutions are:

$$S(t) \approx 66.1\ \ \text{ and } \ \ S(t) \approx 158.9$$

Now we interpret the dynamics for various initial conditions. Consider starting with 40 squirrels, $S_{0} = 40$. With the hawk present, the squirrel population dies off. The reproduction rate of $0.3S$ squirrels per week is not sufficient to offset the impact of the hawk. With an initial population above 66, the squirrels will survive, and their population will approach the equilibrium of 159 squirrels. Consider the original equilibrium at 225 squirrels, i.e. $S_{0} = 225$. When the hawk moves into the park, the squirrel population decreases toward the new equilibrium of 159.

8.4 Falling Bodies

Consider the governing equation for objects falling under the constant force of gravity. If we define the height of an object as $h(t)$, then we have

$$h^{''}(t) = - g$$

where $t$ is time and $g$ is the acceleration due to gravity. If we let $v_{0}$ be the initial, vertical speed and $h_{0}$ be the initial height of an object, then we can solve the resulting initial value problem (IVP): $$h^{''}(t) = - g\ \ \text{ with } \ h(0) = h_{0}\ \text{ and } \ h^{'}(0) = v_{0}$$

Integrating both sides of the ODE, $h^{''}(t) = - g$, with respect to $t$ yields

$$h^{'}(t) = - gt + C_{1}$$

Applying the initial speed gives the value for the constant, $C_{1}$:

$$h^{'}(0) = - g*0 + C_{1} = v_{0}$$

$$C_{1} = v_{0}$$

At this point we have

$$h^{'}(t) = - gt + v_{0}$$

Now integrating both sides of the ODE, $h^{'}(t) = - gt + v_{0}$, with respect to $t$ yields

$$h(t) = - \frac{g}{2}t^{2} + v_{0}t + C_{2}$$

Applying the initial height gives the value for the constant, $C_{2}$:

$$h(0) = - \frac{g}{2}*0^{2} + v_{0}*0 + C_{2} = h_{0}$$

$$C_{2} = h_{0}$$

Combining all the results together generates a formula that is familiar to most Physics students:

$$h(t) = - \frac{g}{2}t^{2} + v_{0}t + h_{0}$$

This straightforward example shows how differential equations are directly related to Physics.

For the rest of this section we focus on modeling with the autonomous ODE, $y^{'} = f(y)$. While the use of generic variable names, e.g. $y$ and $t$, is typical in mathematics, we always want to use descriptive variables when modeling and working with application problems.

Following are examples of modeling using growth, decay, and the logistic equation.

Example #1

During a summer picnic, potato salad is left outside on a table. The bacteria in the salad reproduce rapidly as the potato salad warms. If the bacteria population doubles every 20 minutes, how long will it take 1,000 bacteria to become 1,000,000,000,000 bacteria?

We start by modeling the rate of change of the bacteria population. Let $B(t)$ represent the bacteria population and let $t$ be the independent variable measuring time in hours. The rate of change of bacteria is proportional to the number of bacteria, so the model is

$$B^{'} = kB$$ We can easily calculate the general solution to this autonomous ODE as we did in Section 5. The solution is
$$B(t) = Ce^{kt}$$ The initial condition is given as $B(0) = 1000$, and we can calculate the value for $C$ as
$${B(0) = 1000 = Ce^{k*0} }{C = 1000}$$

We calculate the growth rate, $k$, by applying the fact that the population doubles every 20 minutes (1/3 hours). Given any population, $B^{\blacksquare}$, we know that in 1/3 of an hour the population will double and become $2B^{\blacksquare}$. We use this fact to solve for $k$:

$${2B^{\blacksquare} = B^{\blacksquare}e^{k\left( \frac{1}{3} \right)} }{2 = e^{k\left( \frac{1}{3} \right)} }$$ To solve for the exponent, we apply the natural logarithm to both sides of the equation:

$${\ln(2) = \ln\left( e^{k\left( \frac{1}{3} \right)} \right) = k\left( \frac{1}{3} \right) }{k = 3\ln(2)}$$

The resulting model for the bacterial population is
$$B(t) = 1000e^{3\ln(2)\ t}$$ Finally, we want to solve for the value of $t$ when $B(t) = 1,000,000,000,000$.
$${1000000000000 = 1000e^{3\ln(2)\ t} }{1000000000 = e^{3\ln(2)\ t} }{\ln(1000000000) = 3\ln(2)t }{t = \frac{\ln(1000000000)}{3\ln(2)} }{t \approx 9.966\ hours }$$ The initial population of 1000 bacteria grows to 1 trillion bacteria in just under 10 hours.

Example #2

Create a model for the temperature of a hot cup of tea sitting on a counter at room temperature. Note that the temperature of the tea changes at a rate proportional to the difference between the temperature of the tea and the temperature of the air surrounding the tea.

We start by recognizing that the temperature of the tea will decay from hot to colder. Next, we use the fact that the rate of change is proportional to a temperature difference. Let the temperature of the tea by given by $T(t)$ where $t$ is time measured in minutes. Also, let the temperature of the air be given by $T_{A}$. We assume $T_{A}$ is constant. Using the ideas for modeling proportional decay, we have

$$T^{'} = - \lambda\left( T - T_{A} \right)$$

This model says, “The rate-of-chance of $T$ is negative and is proportional to the temperature difference between $T$ and $T_{A}$ .” This is exactly what we want.

Example #3

For some subjects the rate at which a person learns new material is proportional to the proportion heshe has left to learn. What is a model for this type of learning?

Let the proportion already learned by given by $L(t)$ where $t$ is time in hours. Note that this proportion is bounded between 0 (knowing nothing about the subject) and 1 (knowing everything about the subject). Thus, $0 \leq L(t) \leq 1$ for any value of $t$. An appropriate model is

$$L^{'} = k(1 - L)$$

where $k$ is the learning rate. This model says, “The rate of change of learning is proportional to the amount remaining to be learned.” Amounts can be scaled to proportions.

Example #4

Consider a student with the flu returns to an isolated college campus with 4000 students. Assume all the students interact equally and on a constant basis, and all students will end up catching the flu once. What is a model/IVP for the number of people who have been infected by the flu?

In this case, the flu will spread exponentially as the college students interact. Each infected person will infect some proportion of those with whom heshe interacts. However, the number of students who have been infected is limited to 4000. This situation can be modeled with the logistic equation.

Let $F(t)$ be the number of people who have caught the flu at some point where $t$ is time measured in days. The logistic equation gives the model

$$F^{'} = kF\left( 1 - \frac{F}{4000} \right)$$

where $k$ represents the rate at which the flu spreads. This model says, “The rate of change of students who have contracted the flu is proportional to the product of the number of students who have already caught the flu and the number who have not yet caught the flu.” The product represents the interaction between those with the flu and those without the flu.

Example #5

Newton’s Second Law, $F = ma$, states that force equals the product of mass and acceleration. Moreover, acceleration is the time rate of change of velocity, i.e. $a = v^{'}$. An object falling through the atmosphere has acceleration due to gravity, $g = - 32.2\frac{ft}{s^{2}}$ , but the object faces a friction force due to air resistance. If the air resistance force is proportional to the objects’ velocity, what equation models the velocity of an object falling through the atmosphere?

Let $v(t)$ be the object’s velocity where $t$ is time in seconds. We model the forces:

$${Force\ on\ the\ falling\ object = Force\ due\ to\ gravity + Force\ from\ air\ resistance }{ma = mg - kv }{mv^{'} = mg - kv }$$ The air resistance force is opposite the force due to gravity. We simplify the model to reach our standard form for ODEs

$$v^{'} = g - \frac{k}{m}v$$

The model says, “The rate of change of velocity (acceleration) equals the acceleration due to gravity less a value proportional to velocity.”

Example #6

Initially 1g of a radioactive substance is present. After 12 hours, the mass decreased by 2%. If the rate of decay is proportional to the amount of the substance present at time $t$, then how much of the radioactive substance is remaining after 1 week?

Let $S(t)$ be the amount of substance in milligrams where $t$ is time measured in days. Since the rate of change is proportional to the amount of substance and the substance is decaying, we can use
$$S^{'} = - \lambda S$$

where $\lambda$ is the rate of decay measured in $days^{- 1}$. We can solve this ODE:
$${\frac{dS}{dt} = - \lambda S }{\int_{}^{}{\frac{1}{S}dS} = \int_{}^{}{- \lambda dt} }{\ln(S) + C_{1} = - \lambda t + C_{2} }{\ln(S) = - \lambda t + C_{3} }{S = e^{- \lambda t + C_{3}} }{S = e^{C_{3}}e^{- \lambda t} }{S(t) = Ce^{- \lambda t} }$$ Applying the initial condition (recalling $S$ is measured in mg) gives
$${S(0) = 1000 = Ce^{- \lambda*0} }{C = 1000}$$

After 1/2 day the substance has been reduced by 2%. Therefore, we can find the decay rate as
$${0.98*1000 = 1000e^{- \lambda\left( \frac{1}{2} \right)} }{0.98 = e^{- \frac{\lambda}{2}} }{\ln(0.98) = - \frac{\lambda}{2} }{\lambda = - 2\ln(0.98) }$$ The complete model for the decaying substance is
$$S(t) = 1000e^{2\ln(0.98)\ t}$$ Now we evaluate the function at $t = 7$ days.
$${S(7) = 1000e^{2\ln(0.98)*7} }{S(7) \approx 753.642mg}$$

Of the original 1000mg, there are 753.64mg remaining after 1 week.

Example #7

Two chemicals, P and Q, react with one another to form the compound Z. The chemical formula is $1P + 2Q \rightarrow 1Z$. The rate of the reaction is proportional to the product of the amounts of P and Q that have not yet reacted. Initially, there are 30g of P and 50g of Q. What model can be used to determine the rate at which compound Z is formed?

Let $Z(t)$ be the amount of compound Z in grams where $t$ is time measured in minutes. We know the rate of the reaction depends on the product of the remaining portions of P and Q. The chemical formula, $1P + 2Q \rightarrow 1Z$, tells us that 1g of P along with 2g of Q makes 3g or Z. Thus, at any time the remaining amount of P is: $30 - \frac{1}{3}Z(t)$ grams. Similarly, at any time the amount remaining for Q is: $50\ –\frac{2}{3}Z(t)$ grams. Using the information about rate the rate of the reaction, we have the model

$$Z^{'} = k\left( 30 - \frac{1}{3}Z \right)\left( 50 - \frac{2}{3}Z \right)$$

This model says, “The rate of producing compound Z (the rate of the reaction) is proportional to the product of the remaining amount of P and the remaining amount of Q.” Our model is autonomous, and we can solve for $Z(t)$ as we did in Section 5:

$${\frac{dZ}{dt} = \frac{k}{9}(90 - Z)(150 - 2Z) }{\int_{}^{}{\frac{1}{(90 - Z)(150 - 2Z)}dZ} = \int_{}^{}{\frac{k}{9}dt}}$$

The integral on the left requires partial fraction decomposition, and we state the final solution here:

$$Z(t) = \frac{75Ce^{\frac{10}{3}kt} - 90}{Ce^{\frac{10}{3}kt} - 1}$$

$C$ is found by applying the initial condition, $Z(0) = 0$, since no compound Z exists at the beginning of the reaction. The rate parameter $k$ is found using experimental information about the reaction, such as 10g of Z is formed during the first 5 minutes.

8.5 Exercises

1. A cadet can invest part of his/her Firstie car loan. Assume the money can be placed in a secure investment that earns 6% compounded continuously. Note: interest rates are always given as annual rates. Complete these 3 steps: a) Develop a differential equation that models the value of the investment, $I$, b) Given the initial value, $I_{0} = \$ 12,000$, find the investment’s value as a function of time, c) State the investment’s value after a 25-year AF career, i.e. after 26 years.

2. Consider Example #3 in this section regarding learning rate. Given a person knows everything that can be known about a subject, develop a differential equation model that could be used to represent forgetting. Describe your model and your thought process for developing the model.

3. Two chemicals, P and Q, react with one another to form the compound Z. The chemical formula is $2P + 3Q \rightarrow 1Z$. The rate of the reaction is proportional to the product of the amounts of P and Q that have not yet reacted. Initially, there are 300g of P and 600g of Q. What model can be used to determine the rate at which compound Z is formed?

4. Consider that a modern SUV’s value, $V$, changes at a rate that is directly proportional to its current value. The constant of proportionality for the SUV’s decreasing value is

$\lambda = 0.023\ months^{- 1}$. Complete these 3 steps: a) Provide the differential equation model for the SUV’s value, $V$, b) Given the initial value, $V_{0} = \$ 56,000$, provide the SUV’s value as a function of time, and c) Provide the approximate value of the SUV after 60 months.

5. A Ferrari 428 has a top speed of 195 mph, and its acceleration is proportional to the difference between its current speed and 195 mph. Complete these 3 steps: a) Provide a differential equation model for the speed, $v$, of the Ferrari, b) Given the Ferrari starts at rest and can reach 100 mph in 7 seconds, find the function that represents the Ferrari’s speed, $v(t)$, c) Determine the amount of time required for the Ferrari to reach 180 mph.

6. Consider Example #5 in this section regarding air resistance. The force due to air resistance is actually proportional to the square of speed and the exposed area of the object. Develop a differential equation model that accurately models the velocity of a falling object based upon the square of velocity and area of the object. Explain your model.

7. The population of deer on the Air Force Academy does not face any natural predators. The population of deer, $D$, changes at a rate proportional to the population. The growth rate of deer is $k = 0.4\ year^{- 1}$, and we assume there is no limit to the number of deer USAFA can hold. Each year, hunters harvest 35 deer. Complete these 3 steps: a) Provide a differential equation model for the number of deer, $D$, b) Given an initial deer population of 100, state the deer population as a function of time, c) State the deer population after 5 years.

8. Metabolism of acetaminophen (Tylenol) is an important consideration for patients. The rate at which acetaminophen is metabolized is directly proportional to the amount of acetaminophen, $A$, in a person’s body. The half-life of this drug is 2.5 hours. Complete these 3 steps: a) Provide the differential equation model for the amount of acetaminophen, $A$, in a person’s body, b) Given the initial amount, $A_{0} = 500mg,\ $provide the amount of acetaminophen as a function of time, c) Provide the length of time required for the amount to fall below 70mg.

9. Consider the population of Cutthroat Trout in the Colorado River between Grand Junction and Fruita (~10 miles). The population of Cutthroat Trout, $T$, changes at a rate proportional to the population. The growth rate of cutthroat is $k = 2\ year^{- 1}$, and the 10 miles of river can hold a maximum of 12000 fish. Complete these 3 steps: a) Provide the differential equation model for the Cutthroat population, $T$, along this 10-mile stretch of river, b) Plot the phase line (and direction field?) for the model and state the equilibrium solutions. c) Given an initial population of 500 Cutthroat, use the direction field to calculate how long it will take for the population of fish to reach 11000?

10. Continue using the background information from problem 9). Complete these 3 steps:
    

    1. Modify your model from part 9a) by including the fact that fishermen and Osprey catch and remove approximately 4000 Cutthroat each year, b) Determine the equilibrium solutions for this updated model and plot the phase line of the model. What is the minimum sustainable population of Cutthroat? c) (more difficult) … What is the maximum number of Cutthroat that can be taken by fishermen and Osprey each year?