Calculus tells us that the derivative of a function measures how the function changes. An equation relating a function to one or more of its derivatives is called a differential equation. The subject of differential equations is one of the most interesting and useful areas of mathematics. We can describe many interesting natural phenomena that involve change using differential equations. In addition, the theory of the subject has broad and important implications.

We begin our study of ordinary differential equations by modeling some real world phenomena. For a particular situation that we might wish to investigate, our first task is to write an equation (or equations) that best describes the phenomenon. Suppose that we wish to study how a population \(P(t)\) grows with time \(t\text{.}\) We might make the assumption that a constant fraction of population is having offspring at any particular time. If we also assume that the population has a constant death rate, the change in the population \(\Delta P\) during a small time interval \(\Delta t\) will be

\begin{equation*}
\Delta P \approx k_{\text{birth}} P(t) \Delta t - k_{\text{death}} P(t) \Delta t,
\end{equation*}

where \(k_{\text{birth}}\) is the fraction of the population having offspring during the interval and \(k_{\text{death}}\) is the fraction of the population that dies during the interval. Equivalently, we can write

\begin{equation*}
\frac{\Delta P}{ \Delta t} \approx k P(t),
\end{equation*}

where \(k = k_{\text{birth}} - k_{\text{death}}\text{.}\) Since the derivative of \(P\) is

is one of the simplest differential equations that we will consider. The equation tells us that the population grows in proportion to its current size. It is not too difficult to see that \(P(t) = Ce^{kt}\) is a solution to this equation, where \(C\) is an arbitrary constant. Indeed, if we differentiate \(P(t)\text{,}\) we obtain

In addition, if we know the value of \(P(t)\text{,}\) say when \(t = 0\text{,}\) we can also determine the value of \(C\text{.}\) For example, if the population at the time \(t = 0\) is \(P(0) = P_0\text{,}\) then

is an example of an initial value problem, and we say that \(P(0) = P_0\) is an initial condition. Since the solution to equation (1.1.1) is \(P(t) = Ce^{kt}\text{,}\) and we say that the population grows exponentially.

As an example, suppose that \(P(t)\) is a population of a colony of bacteria at time \(t\text{,}\) whose initial population is 1000 at \(t = 0\text{,}\) where time is measured in hours. Then

\begin{equation*}
1000 = P(0) = C e^0 = C,
\end{equation*}

and our solution becomes \(P(t) = 1000e^{kt}\text{.}\) If the population grows at three percent per hour, then

Of course, it is important to realize that this is only a model. If \(t\) is small, our model might be reasonably accurate. However, if we let \(t\) be very large, our colony of bacteria could very well exceed the mass of the earth.

The growth rate of a population need not be positive. For example, Japan has experienced negative growth in recent years.^{ 1 }See http://data.worldbank.org/indicator/SP.POP.GROW. The equation \(dP/dt = kP\) can also be used to model phenomena such as radioactive decay and compound interest—topics which we will explore later.

To summarize, we say that the expression \(x(t) = Ce^{kt}\) is a general solution of the equation \(x' = kx\text{,}\) and \(x(t) = x_0 e^{kt}\) is a particular solution to the differential equation. The general solution to our equation \(x(t) = Ce^{kt}\) graphs as an infinite family of curves, which we will call integral curves or solution curves (Figure 1.1.1).

Not all populations grow exponentially; otherwise, a bacteria culture in a petri dish would grow unbounded and soon be much larger than the size of the laboratory. To see what happens if there are limiting factors to population growth, let us consider the population of fish in a children's trout pond. The number of trout will be limited by the available resources such as food supply as well as by spawning habitat. A small population of fish might grow exponentially if the pond is large and food is abundant, but the growth rate will decline as the population increases and the availability of resources declines. We can use the logistic equation to model population growth in a resource limited environment.^{ 2 }The logistic model was first used by the Belgian mathematician and physician Pierre François Verhulst in 1836 to predict the populations of Belgium and France.

To see how the logistic model works, let us try to adjust our model of exponential growth to account for the limited resources of the pond. We will make the following assumptions.

If the population of trout is small and the pond is large with abundant resources, the rate of growth will be approximately exponential,

for \(P \gt N\text{.}\) We say that \(N\) is the carrying capacity for the population.

Our assumptions suggest that we might try an equation of the form

\begin{equation*}
\frac{dP}{dt} = k f(P) P,
\end{equation*}

where \(f(P)\) is a function of \(P\) that is close to 1 if the population is small, but negative if the population is greater than \(N\text{.}\) The simplest function satisfying these conditions is

Thus, the logistic population model is given by the differential equation

\begin{equation*}
\frac{dP}{dt} = k \left( 1 - \frac{P}{N} \right) P.
\end{equation*}

Example1.1.2

Suppose we have a pond that will support 1000 fish, and the initial population is 100 fish. In order to determine the number of fish in the lake at any time \(t\text{,}\) we must find a solution to the initial value problem

It is easy to verify that \(P(t) = 1000/(9e^{-kt} + 1)\) is the solution to our initial value problem.^{ 3 }We will learn how to solve initial-value problems such as the one described here in Section 1.2. For the time being, we will be satisfied with being able to verify the fact that we have a solution. Certainly \(P(0) = 100\text{,}\) and if we differentiate \(P\text{,}\) we will obtain the righthand side of the differential equation,

Sometimes it is necessary to consider the second derivative when modeling a phenomenon. Suppose that we have a mass lying on a flat, frictionless surface and that this mass is attached to one end of a spring with the other end of the spring attached to a wall. We will denote displacement of the spring by \(x\text{.}\) If \(x \gt 0\text{,}\) then the spring is stretched. If \(x \lt 0\text{,}\) the spring is compressed. If \(x = 0\text{,}\) then the spring is in a state of equilibrium (Figure 1.1.4). If we pull or push on the mass and release it, then the mass will oscillate back and forth across the table.

We can construct a differential equation that models our oscillating mass. First, we must consider the restorative force on the spring. We will make the assumption that this force depends on the displacement of the spring, \(F(x)\text{.}\) Using Taylor's Theorem from calculus, we can expand \(F\) to obtain

where \(F'(0) = -k\) and \(F(0) = 0\text{.}\) If the displacement is not too large, then \(x^n\) will be small for \(n \geq 2\text{,}\) and we can ignore higher ordered terms. Thus, we can consider the restorative force on the spring to be proportional to displacement of the spring from its equilibrium length,

\begin{equation*}
F = -kx.
\end{equation*}

This equation is known as Hooke's Law. We can test this law experimentally, and it is reasonably accurate if the displacement of the spring is not too large.

By Newton's second law of motion, the force on the mass must be

\begin{equation*}
F = ma = m \frac{d^2 x}{dt^2} = m x''.
\end{equation*}

Setting the two forces equal, we have a second order differential equation,

\begin{equation*}
mx'' = -kx,
\end{equation*}

which describes our oscillating mass. The spring-mass system is an example of a harmonic oscillator.

Example1.1.5

Suppose that we have a spring-mass system where \(m =1\) and \(k = 1\text{.}\) If the initial velocity of the spring is one unit per second and the initial position is at the equilibrium point, then we have the following initial value problem,

Since \(x''(t) = - x(t)\) for both the sine and cosine functions, we might guess that a general solution of our differential equation has the form

\begin{equation*}
x(t) = A \cos t + B \sin t.
\end{equation*}

Noting that

\begin{equation*}
x'(t) = -A \sin t + B \cos t,
\end{equation*}

and using our initial conditions, we can determine that \(A = 0\) and \(B = 1\) or

\begin{equation*}
x(t) = \sin t.
\end{equation*}

The graph of the position of the mass as a function of time is given in Figure 1.1.6.

Now let us add a damping force to our system. For example, we might add a dashpot, a mechanical device that resists motion, to our system. Think of a dashpot as that small cylinder that keeps your screen door from slamming shut. The simplest assumption would be to take the damping force of the dashpot to be proportional to the velocity of the mass, \(x'(t)\text{.}\) In other words, the harder you try to slam the screen door, the more resistance you will feel. Thus, we have will have an additional force,

\begin{equation*}
F = -b x'
\end{equation*}

acting on our mass, where \(b \gt 0\text{.}\) Our new equation for the spring-mass system is

Here we let \(m = 1\text{,}\) \(b = 3\text{,}\) and \(k = 2\text{,}\) We will learn how to solve equations of the form \(mx'' + bx' + kx = 0\) in Chapter 4, but let us assume that the solution is of the form \(x(t) = e^{rt}\) for now. In this case,

Since \(e^{rt}\) is never zero, it must be the case that \(r = -2\) or \(r = -1\text{,}\) if \(x(t) = e^{rt}\) is to be a solution to our equation. Thus, we might guess that

\begin{equation*}
x(t) = A e^{-t} + B e^{-2t}
\end{equation*}

is a general solution to our equation. If the initial velocity of our mass is one unit per second and the initial position is zero, then we have the initial value problem

Notice that the additional damping negates any oscillation in the system. In this case, we say that the system is over-damped (Figure 1.1.8).

Of course, if we have a very strong spring and only add a small amount of damping to our spring-mass system, the mass would continue to oscillate, but the oscillations would become progressively smaller. In other words, our system would be under-damped. For example, our spring-mass system might be described by the initial value problem

Some situations require more than one differential equation to model a particular phenomenon. We might use a system of differential equations to model two interacting species, say where one species preys on the other.^{ 4 }The predator-prey model was discovered independently by Lotka (1925) and Volterra (1926). For example, we can model how the population of Canadian lynx (lynx canadenis) interacts with a the population of snowshoe hare (lepus americanis) (see https://www.youtube.com/watch?v=ZWucOrSOdCs).

We have good historical data for the populations of the lynx and snowshoe hare from the Hudson Bay Company. This large Canadian retail company, which owns and operates a large number of retail stores in North America and Europe, including Galeria Kaufhof, Hudson's Bay, Home Outfitters, Lord & Taylor, and Saks Fifth Avenue, was originally founded in 1670 as a fur trading company. The Hudson Bay Company kept accurate records on the number of lynx pelts that were bought from trappers from 1821 to 1940. The company noticed that the number of pelts varied from year to year and that the number of lynx pelts reached a peak about every ten years [11]. The ten year cycle for lynx can be best understood using a system of differential equations.

The primary prey for the Canadian lynx is the snowshoe hare. We will denote the population of hares by \(H(t)\) and the population of lynx by \(L(t)\text{,}\) where \(t\) is the time measured in years. We will make the following assumptions for our predator-prey model.

If no lynx are present, we will assume that the hares reproduce at a rate proportional to their population and are not affected by overcrowding. That is, the hare population will grow exponentially,

Since the lynx prey on the hares, we can argue that the rate at which the hares are consumed by the lynx is proportional to the rate at which the hares and lynx interact. Thus, the equation that predicts the rate of change of the hare population becomes

The lynx receive benefit from the hare population. The rate at which lynx are born is proportional to the number of hares that are eaten, and this is proportional to the rate at which the hares and lynx interact. Consequently, the growth rate of the lynx population can be described by

We will learn how to analyze and find solutions of systems of differential equations in subsequent chapters; however, we will give a graphical solution in Figure 1.1.10 to the system

Our graphical solution is obtained using a numerical algorithm (see Section 1.4 and Section 2.3). Notice that the predator population, \(L\text{,}\) begins to grow and reaches a peak after the prey population, \(H\) reaches its peak. As the prey population declines, the predator population also declines. Once the predator population is smaller, the prey population has a chance to recover, and the cycle begins again.^{ 5 }An excellent account of the actual lynx and snowshoe hare data and model can be found in [5].

The interaction of the HIV-1 virus with the body's immune system can be modeled by a system of differential equations similar to a predator-prey system. After an individual is infected with the HIV-1 virus, the amount of the virus in the bloodstream rises dramatically and the person will often suffer from flu-like symptoms. However, these symptoms will disappear after a period of weeks or months as the body begins to manufacture antibodies against the virus. Tests have been developed to determine the presence of HIV-1 antibodies. If an individual has such antibodies, then they are said to be HIV-1 positive. Once infected with the HIV-1 virus, it can be years before an HIV-positive patient exhibits the full symptoms of AIDS. The body's immune system fights the HIV-1 virus with white blood cells. The CD4-positive T-helper cell, a specific type of white blood cell, is especially important since it helps other cells fight the virus. However, the HIV-1 virus use the CD4-positive T-helper cells to create more virions, destroying the CD4-positive T-helper cells in the process.

We can develop a system of differential equations to better understand the dynamics of the HIV-1 virus [20]. Let \(V = V(t)\) be the population of the HIV-1 virus at time \(t\text{.}\) We will assume that the virus concentration is governed by the following differential equation,

\begin{equation*}
\frac{dV}{dt} = P - cV.
\end{equation*}

The first term, \(P\) is some function of \(t\) that determines the rate at which new viral particles are created. The term \(-cV\) is the death rate for the virions. If someone discovers a drug that blocks the creation of new HIV-1 virions, then \(P\) would be zero and the virions would clear the body at the following rate,

and \(V(t) = V_0 e^{-ct}\text{,}\) where \(V_0\) is the initial viral population.

Now let us consider a model for the concentration \(T = T(t)\) of (uninfected) CD4-positive T-helper cells,

\begin{equation*}
\frac{dT}{dt} = s + pT\left(1 - \frac{T}{T_{\text{max}}} \right) - d_T T.
\end{equation*}

The constant \(s\) represents the rate at which T-cells are created from sources in the body, such as the thymus. New CD4-positive T-helper cells can also be created from the proliferation of existing CD4-positive T-helper cells, and the second term in the equation represents the logistic growth of the T-cells, where \(p\) is the maximum proliferation rate and \(T_{\text{max}}\) is the T-cell population density where proliferation ceases. Finally, \(d_T\) is the death rate of the T cells.

Like the influenza virus, the HIV-1 virus is an RNA virus. An RNA virus cannot reproduce on its own and must use the DNA from a host cell. To do this, the virus attaches itself to a CD4-positive T-helper cell and injects its RNA into the cell. This way the virus can use the T-cell's DNA to replicate itself using a process called reverse transcription, where a DNA copy of the virus's RNA is made. New virus particles are created, and the T-cell eventually bursts releasing the virions into the body. If we let \(T^*\) be the concentration of infected T-cells, we can model this process with the following system of equations,

where \(\delta\) is the rate of loss of the virus producing T-cells and \(N\) is the number of virions produced per infected T-cell during its lifetime. The term \(kTV\) tells us the rate at which the HIV-1 virus infects T-cells. This is the same idea as modeling how predators interact with prey in a predator-prey model. Thus, our complete model becomes

One class of drugs that HIV infected patients receive are reverse transcriptase (RT) inhibitors. RT inhibitors block the action of reverse transcription and prevent the virus from duplicating. If one could find the perfect RT inhibitor, then \(k =0\) and our system becomes

Unfortunately, no one has discovered a perfect RT inhibitor, so we will need to modify the system to account for the effectiveness of the RT inhibitor. We can accomplish this by adding an effectiveness factor, \(1 - \eta\text{,}\) to the \(kVT\) term. Our system now becomes

If \(\eta = 1\text{,}\) then the RT inhibitor is completely effective. On the other hand, if \(\eta = 0\text{,}\) then the RT inhibitor is completely ineffective. We now have a model for how the HIV-1 virus interacts with the immune system. Researchers can use data to estimate the parameters and see exactly what types of solutions are possible.

In this section we have provided a general notion of what a differential equation is as well as several modeling situations where differential equations are useful; however, we have left many questions unanswered. First, we need a more rigorous definition of a differential equation. What is the proper way to define a system of differential equations? Does a differential equation or a system of differential equations always have a solution? Are solutions to differential equations unique? If a unique solution to a differential equation exists, can we find it? If it is not possible to find a precise solution algebraically, can we estimate the solution numerically? If neither is possible, can we still say anything useful about the solution? Of course, other questions will come to mind as we continue our study of differential equations.

where \(N\) is the carrying capacity of the population.

Some phenomenon, such as the relationship between a population of predators and a population of prey, are best modeled by systems of differential equations.

The three principle steps in modeling any phenomenon with differential equations are:

Discovering the differential equation or equations that best describe a specified physical situation.

Finding—either exactly or approximately—the appropriate solution of the equation or equations.

Interpreting the solution in terms of the phenomenon.

Many differential equations have solutions of the form \(y(t) = e^{at}\text{,}\) where \(a\) is some constant. For example, \(y(t) = e^{3t}\) is a solution to the equation \(y' = 3y\text{.}\) Find all values of \(a\) such that \(y(t) = e^{at}\) is a solution to the given equation in Exercise Group 9–14.

9

\(y' = -2y\)

10

\(y' + 7 y = 0\)

11

\(y'' - 3y' + 2y = 0\)

12

\(y'' + 4y' - 12y = 0\)

13

\(y''' - y'' - 4y' + 4y = 0\)

14

\(y^{(5)} - 5y''' + 4y' = 0\)

Initial Value Problems

Use direct substitution to verify that \(y(t)\) is a solution of the given differential equation in Exercise Group 15–20. Then use the initial conditions to determine the constants \(C\) or \(c_1\) and \(c_2\text{.}\)

Use direct substitution to verify that \(y(t)\) is a solution of the given differential equation in Exercise Group 21–24. Then use the boundary conditions to determine the constants \(c_1\) and \(c_2\) (if possible).

Consider the differential equation \(y'' + 9y = 0\text{.}\)

Verify that \(y(t) = c_1 \cos 3t + c_2 \sin 3t\) is a solution to this equation.

Sketch solution curves for \(c_1 = 1\) and \(c_2 = 1, 2, \ldots, 5\text{.}\)

27

The growth of a population of rabbits with unlimited resources and space can be modeled by the exponential growth equation, \(dP/dt = kP\text{.}\)

Write a differential equation to model a population of rabbits with unlimited resources, where hunting is allowed at a constant rate \(\alpha\text{.}\)

Write a differential equation to model a population of rabbits with unlimited resources, where hunting is allowed at a rate proportional to the population of rabbits.

Write a differential equation to model a population of rabbits with limited resources, where hunting is allowed at a constant rate \(\alpha\text{.}\)

Write a differential equation to model a population of rabbits with limited resources, where hunting is allowed at a rate proportional to the population of rabbits.

28

We can modify the logistic growth model to account for a threshold population. That is, if the population of a species drops below a certain level, the species is endanger of extinction. Consider the case of the black rhinoceros, once the most numerous of all rhinoceros species and now critically endangered. The black rhino, native to eastern and southern Africa, was estimated to have a population of about 100,000 around 1900. Because of hunting, habitat changes, competing species, and most of all illegal poaching, the number of black rhinos today is estimated to be below 3000. If the wild population becomes too low, the animals may not be able to find suitible mates and the black rhino will become extinct. There must be a minimum population for the species to continue. Suppose the this minimum or threshold population for the black rhino is 1000 animals and that remaining habitant in Africa will support no more that 20,000 rhinos.

For what values of \(P\) is the rhino population increasing? What can be said about the value of \(dP/dt\) for these values of \(P\text{?}\)

For what values of \(P\) is the rhino population decreasing? What can be said about the value of \(dP/dt\) for these values of \(P\text{?}\)

For what values of \(P\) is the rhino population in equilibrium? What can be said about the value of \(dP/dt\) for these values of \(P\text{?}\)

Find a differential equation that models the population of rhinos at time \(t\text{.}\)

(a) The population is increasing if \(dP/dt \gt 0\) and \(1000 \lt P \lt 20000\text{.}\)

29

Given the equation \(x' + p x = q(t)\text{,}\) where \(p\) is a constant and \(q(t)\) is a continuous function defined on an interval \(I\text{,}\) show that

Carbon 14 is a radioactive isotope of carbon, the most common isotope of carbon being carbon 12. Carbon 14 is created when cosmic ray bombardment changes nitrogen 14 to carbon 14 in the upper atmosphere. The resulting carbon 14 combines with atmospheric oxygen to form radioactive carbon dioxide, which is incorporated into plants by photosynthesis. Animals acquire carbon 14 by eating plants. When an animal or plant dies, it ceases to take on carbon 14, and the amount of isotope in the organism begins to decay into the more common carbon 12. Carbon 14 has a very long half-life, about 5730 years. That is, given a sample of carbon 14, it will take 5730 years for half of the sample to decay to carbon 12. The long half-life is what makes carbon 14 dating very useful in dating objects from antiquity.

Consider a sample of material that contains \(A(t)\) atoms of carbon 14 at time \(t\text{.}\) During each unit of time a constant fraction of the radioactive atoms will spontaneously decay into another element or a different isotope of the same element. Thus, the sample behaves like a population with a constant death rate and a zero birth rate. Make use of the model of exponential growth to construct a differential equation that models radioactive decay for carbon 14.

Solve the equation that you proposed in (a) to find an explicit formula for \(A(t)\text{.}\)

The Chauvet-Pont-d'Arc Cave in the Ardèche department of southern France contains some of the best preserved cave paintings in the world. Carbon samples from torch marks and from the paintings themselves, as well as from animal bones and charcoal found on the cave floor, have been used to estimate the age of the cave paintings. If a particular sample taken from the Cauvet Cave contains 2% of the expected cabon 14, what is the approximate age of the sample?

31

Consider the following predator-prey systems of differential equations

Technology can prove very useful when studying differential equations. Software packages such Maple, Mathematica, and Matlab each have their advantages and disadvantages. We will use Sage, a readily available open source computer algebra system, as our choice of software. Sage can be run on an individual computer or over the Internet on a server. You can even access Sage from your smart phone. For our purposes, Sage cells are embedded into the textbook, so there is nothing to install on your computer. Simply, evaluate the cell. You can even change the preloaded commands in the cell if you wish. For example, let us evaluate the derivative of \(f(x) = x^2 \cos x\text{.}\)

Note that anything following a pound sign # is a comment.

We can use Sage to plot functions. For example, we can plot the function \(f(x) = x^2 \cos x\) as well as its derivative on the same graph.

We can use Sage to solve differential equations. Suppose that we wish to solve the initial value problem

We will provide abundant examples of how to use Sage to solve and analyze differential equations throughout the book, and we encourage the reader to experiment by altering the Sage commands inside the individual Sage cells. If you make a mistake, you can simply reload the webpage and start again.

The reader will find plenty of resources to learn how to use Sage. A good place to start is http://www.sagemath.org/help.html or [1]. Although we will be using Sage as the technology of choice, much of this book can be read independently of Sage. Finally, we would like to emphasize once again that the reader who chooses not to use some sort of technology will be at a disadvantage.

SubsubsectionExercises

1

In the Sage cell below enter 2 + 2 and then evaluate the cell. Your answer should be 4 of course. Next try entering 2^1000. You should see a very large number.