2.9 Modeling with Systems of First-Order DEs

INTRODUCTION

In this section we are going to discuss mathematical models based on some of the topics that we have already discussed in the preceding two sections. This section will be similar to Section 1.3 in that we are only going to discuss systems of first-order differential equations as mathematical models and we are not going to solve any of these models. There are two good reasons for not solving systems at this point: First, we do not as yet possess the necessary mathematical tools for solving systems, and second, some of the systems that we discuss cannot be solved analytically. We shall examine solution methods for systems of linear first-order DEs in Chapter 10 and for systems of linear higher-order DEs in Chapters 3 and 4.

Systems

Up to now all the mathematical models that we have considered were single differential equations. A single differential equation could describe a single population in an environment, but if there are, say, two interacting and perhaps competing species living in the same environment (for example, rabbits and foxes), then a model for their populations x(t) and y(t) might be a system of two first-order differential equations such as

(1)

When g1 and g2 are linear in the variables x and y, that is,

g1(t, x, y) = c1x + c2y + f1(t) and g2(t, x, y) = c3x + c4y + f2(t),

then (1) is said to be a linear system. A system of differential equations that is not linear is said to be nonlinear.

Radioactive Series

In the discussion of radioactive decay in Sections 1.3 and 2.7, we assumed that the rate of decay was proportional to the number A(t) of nuclei of the substance present at time t. When a substance decays by radioactivity, it usually doesn’t just transmute into one stable substance and then the process stops. Rather, the first substance decays into another radioactive substance, this substance in turn decays into yet a third substance, and so on. This process, called a radioactive decay series, continues until a stable element is reached. For example, the uranium decay series is U-238 → Th-234 → … → Pb-206, where Pb-206 is a stable isotope of lead. The half-lives of the various elements in a radioactive series can range from billions of years (4.5 × 109 years for U-238) to a fraction of a second. Suppose a radioactive series is described schematically by Z, where k1 = −λ1 < 0 and k2 = −λ2 < 0 are the decay constants for substances X and Y, respectively, and Z is a stable element. Suppose too, that x(t), y(t), and z(t) denote amounts of substances X, Y, and Z, respectively, remaining at time t. The decay of element X is described by

whereas the rate at which the second element Y decays is the net rate,

since it is gaining atoms from the decay of X and at the same time losing atoms due to its own decay. Since Z is a stable element, it is simply gaining atoms from the decay of element Y:

In other words, a model of the radioactive decay series for three elements is the linear system of three first-order differential equations

(2)

Mixtures

Consider the two tanks shown in FIGURE 2.9.1. Let us suppose for the sake of discussion that tank A contains 50 gallons of water in which 25 pounds of salt is dissolved. Suppose tank B contains 50 gallons of pure water. Liquid is pumped in and out of the tanks as indicated in the figure; the mixture exchanged between the two tanks and the liquid pumped out of tank B is assumed to be well stirred. We wish to construct a mathematical model that describes the number of pounds x1(t) and x2(t) of salt in tanks A and B, respectively, at time t.

Two tanks are labeled A and B. 4 pipes are connected to the tanks and the rate of flow is as follows. Into A: pure water at three gallons per minute. From B to A: mixture at four gallons per minute. From A to B: mixture at one gallon per minute. Out of B: mixture at three gallons per minute.

FIGURE 2.9.1 Connected mixing tanks

By an analysis similar to that on page 22 in Section 1.3 and Example 5 of Section 2.7, we see for tank A that the net rate of change of x1(t) is

Similarly, for tank B, the net rate of change of x2(t) is

Thus we obtain the linear system

(3)

Observe that the foregoing system is accompanied by the initial conditions x1(0) = 25, x2(0) = 0.

A Predator–Prey Model

Suppose that two different species of animals interact within the same environment or ecosystem, and suppose further that the first species eats only vegetation and the second eats only the first species. In other words, one species is a predator and the other is a prey. For example, wolves hunt grass-eating caribou, sharks devour little fish, and the snowy owl pursues an arctic rodent called the lemming. For the sake of discussion, let us imagine that the predators are foxes and the prey are rabbits.

Let x(t) and y(t) denote, respectively, the fox and rabbit populations at time t. If there were no rabbits, then one might expect that the foxes, lacking an adequate food supply, would decline in number according to

(4)

When rabbits are present in the environment, however, it seems reasonable that the number of encounters or interactions between these two species per unit time is jointly proportional to their populations x and y; that is, proportional to the product xy. Thus when rabbits are present there is a supply of food, and so foxes are added to the system at a rate bxy, b > 0. Adding this last rate to (4) gives a model for the fox population:

(5)

On the other hand, were there no foxes, then the rabbits would, with an added assumption of unlimited food supply, grow at a rate that is proportional to the number of rabbits present at time t:

(6)

But when foxes are present, a model for the rabbit population is (6) decreased by cxy, c > 0; that is, decreased by the rate at which the rabbits are eaten during their encounters with the foxes:

(7)

Equations (5) and (7) constitute a system of nonlinear differential equations

(8)

where a, b, c, and d are positive constants. This famous system of equations is known as the Lotka–Volterra predator–prey model.

Except for two constant solutions, x(t) = 0, y(t) = 0 and x(t) = d/c, y(t) = a/b, the nonlinear system (8) cannot be solved in terms of elementary functions. However, we can analyze such systems quantitatively and qualitatively. See Chapters 6 and 11.

EXAMPLE 1 Predator–Prey Model

Suppose

represents a predator–prey model. Since we are dealing with populations, we have x(t) ≥ 0, y(t) ≥ 0. FIGURE 2.9.2, obtained with the aid of a numerical solver, shows typical population curves of the predators and prey for this model superimposed on the same coordinate axes. The initial conditions used were x(0) = 4, y(0) = 4. The curve in red represents the population x(t) of the predator (foxes), and the blue curve is the population y(t) of the prey (rabbits). Observe that the model seems to predict that both populations x(t) and y(t) are periodic in time. This makes intuitive sense since, as the number of prey decreases, the predator population eventually decreases because of a diminished food supply; but attendant to a decrease in the number of predators is an increase in the number of prey; this in turn gives rise to an increased number of predators, which ultimately brings about another decrease in the number of prey. ≡

A graph consists of two periodic curves. The horizontal axis is labeled time t, and the vertical axis is labeled population P. The first curve represents prey. One cycle begins at approximately (0, 4), reaches a high point at approximately (0.7, 6.1) and then goes down to the right, reaching a low point at approximately (4.5, 0). The second curve represents predators. One cycle begins at approximately (0, 4), reaches a high point at approximately (2.1, 6.4) and then goes down to the right, reaching a low point at approximately (7.5, 3.8).

FIGURE 2.9.2 Population of predators (red) and prey (blue) appear to be periodic in Example 1

Competition Models

Now suppose two different species of animals occupy the same ecosystem, not as predator and prey but rather as competitors for the same resources (such as food and living space) in the system. In the absence of the other, let us assume that the rate at which each population grows is given by

    and     (9)

respectively.

Since the two species compete, another assumption might be that each of these rates is diminished simply by the influence, or existence, of the other population. Thus a model for the two populations is given by the linear system

(10)

where a, b, c, and d are positive constants.

On the other hand, we might assume, as we did in (5), that each growth rate in (9) should be reduced by a rate proportional to the number of interactions between the two species:

(11)

Inspection shows that this nonlinear system is similar to the Lotka–Volterra predator–prey model. Last, it might be more realistic to replace the rates in (9), which indicate that the population of each species in isolation grows exponentially, with rates indicating that each population grows logistically (that is, over a long time the population is bounded):

    and     (12)

When these new rates are decreased by rates proportional to the number of interactions, we obtain another nonlinear model

(13)

where all coefficients are positive. The linear system (10) and the nonlinear systems (11) and (13) are, of course, called competition models.

Networks

An electrical network having more than one loop also gives rise to simultaneous differential equations. As shown in FIGURE 2.9.3, the current i1(t) splits in the directions shown at point B1, called a branch point of the network. By the first of the two circuit laws formulated by Gustav Robert Kirchhoff, we can write

(14)

A circuit diagram. Resistor R subscript 1 is in series with 2 parallel paths. The first one consists of inductor L subscript 1 and resistor R subscript 2 in series between nodes B subscript 1 and B subscript 2. The second path consists of L subscript 2 between nodes C subscript 1 and C subscript 2. Current i subscript 1 from source E flows through node A subscript 1, resistor R subscript 1 and to node B subscript 1. From there, it branches into i subscript 2 and i subscript 3. i subscript 2 flows through the path with L subscript 1 and r subscript 2 to node B subscript 2. i subscript 3 flows through node C subscript 1, through the path with L subscript 2, to node C subscript 2.

FIGURE 2.9.3 Network whose model is given in (17)

In addition, we can also apply Kirchhoff’s second law to each loop. For loop A1B1B2A2A1, summing the voltage drops across each part of the loop gives

(15)

Similarly, for loop A1B1C1C2B2A2A1 we find

(16)

Using (14) to eliminate i1 in (15) and (16) yields a system of two linear first-order equations for the currents i2(t) and i3(t):

(17)

We leave it as an exercise (see Problem 16 in Exercises 2.9) to show that the system of differential equations describing the currents i1(t) and i2(t) in the network containing a resistor, an inductor, and a capacitor shown in FIGURE 2.9.4 is

(18)

A circuit diagram. Inductor L is in series with 2 parallel paths. The first one consists of resistor R. The second path consists of capacitor C. Current i subscript 1 from source E flows through inductor L. It then branches into i subscript 2 and i subscript 3. i subscript 2 flows through the path with resistor R. i subscript 3 flows through the path with capacitor C.

FIGURE 2.9.4 Network whose model is given in (18)

2.9 Exercises Answers to selected odd-numbered problems begin on page ANS-4.

Radioactive Series

  1. We have not discussed methods by which systems of first-order differential equations can be solved. Nevertheless, systems such as (2) can be solved with no knowledge other than how to solve a single linear first-order equation. Find a solution of (2) subject to the initial conditions x(0) = x0, y(0) = 0, z(0) = 0.
  2. In Problem 1, suppose that time is measured in days, that the decay constants are k1 = −0.138629 and k2 = −0.004951, and that x0 = 20. Use a graphing utility to obtain the graphs of the solutions x(t), y(t), and z(t) on the same set of coordinate axes. Use the graphs to approximate the half-lives of substances X and Y.
  3. Use the graphs in Problem 2 to approximate the times when the amounts x(t) and y(t) are the same, the times when the amounts x(t) and z(t) are the same, and the times when the amounts y(t) and z(t) are the same. Why does the time determined when the amounts y(t) and z(t) are the same make intuitive sense?
  4. Construct a mathematical model for a radioactive series of four elements W, X, Y, and Z, where Z is a stable element.

Contributed Problem

Jeff Dodd, Professor
Department of Mathematical Sciences, Jacksonville State University

  1. Potassium-40 Decay     The mineral potassium, whose chemical symbol is K, is the eighth most abundant element in the Earth’s crust, making up about 2% of it by weight, and one of its naturally occurring isotopes, K-40, is radioactive. The radioactive decay of K-40 is more complex than that of carbon-14 because each of its atoms decays through one of two different nuclear decay reactions into one of two different substances: the mineral calcium-40 (Ca-40) or the gas argon-40 (Ar-40). Dating methods have been developed using both of these decay products. In each case, the age of a sample is calculated using the ratio of two numbers: the amount of the parent isotope K-40 in the sample and the amount of the daughter isotope (Ca-40 or Ar-40) in the sample that is radiogenic, in other words, the substance that originates from the decay of the parent isotope after the formation of the rock.
    A photo of a volcano eruption and lava flowing down the sides.

    © beboy/ShutterStock, Inc.

    An igneous rock is solidified magma

    The amount of K-40 in a sample is easy to calculate. K-40 comprises 1.17% of naturally occurring potassium, and this small percentage is distributed quite uniformly, so that the mass of K-40 in the sample is just 1.17% of the total mass of potassium in the sample, which can be measured. But for several reasons it is complicated, and sometimes problematic, to determine how much of the Ca-40 in a sample is radiogenic. In contrast, when an igneous rock is formed by volcanic activity, all of the argon (and other) gas previously trapped in the rock is driven away by the intense heat. At the moment when the rock cools and solidifies, the gas trapped inside the rock has the same composition as the atmosphere. There are three stable isotopes of argon, and in the atmosphere they occur in the following relative abundances: 0.063% Ar-38, 0.337% Ar-36, and 99.60% Ar-40. Of these, just one, Ar-36, is not created radiogenically by the decay of any element, so any Ar-40 in excess of 99.60/(0.337) = 295.5 times the amount of Ar-36 must be radiogenic. So the amount of radiogenic Ar-40 in the sample can be determined from the amounts of Ar-38 and Ar-36 in the sample, which can be measured.

    Assuming that we have a sample of rock for which the amount of K-40 and the amount of radiogenic Ar-40 have been determined, how can we calculate the age of the rock? Let P(t) be the amount of K-40, A(t) the amount of radiogenic Ar-40, and C(t) the amount of radiogenic Ca-40 in the sample as functions of time t in years since the formation of the rock. Then a mathematical model for the decay of K-40 is the system of linear first-order differential equations

    where

    (a) From the system of differential equations find P(t) if P(0) = P0.

    (b) Determine the half-life of K-40.

    (c) Use P(t) from part (a) to find A(t) and C(t) if A(0) = 0 and C(0) = 0.

    (d) Use your solution for A(t) in part (c) to determine the percentage of an initial amount P0 of K-40 that decays into Ar-40 over a very long period of time (that is, t → ). What percentage of P0 decays into Ca-40?

  2. Potassium–Argon Dating     The knowledge of how K-40 decays can be used to determine the age of very old igneous rocks.

    (a) Use the solutions obtained in parts (a) and (c) of Problem 5 to show that

    (b) Solve the expression in part (a) for t in terms A(t), P(t), λA, and λC.

    (c) Suppose it is found that each gram of a rock sample contains 8.6 × 10−7 grams of radiogenic Ar-40 and 5.3 × 10−6 grams of K-40. Use the equation obtained in part (b) to determine the approximate age of the rock.

Mixtures

  1. Consider two tanks A and B, with liquid being pumped in and out at the same rates, as described by the system of equations (3). What is the system of differential equations if, instead of pure water, a brine solution containing 2 pounds of salt per gallon is pumped into tank A?
  2. Use the information given in FIGURE 2.9.5 to construct a mathematical model for the number of pounds of salt x1(t), x2(t), and x3(t) at time t in tanks A, B, and C, respectively.
    Three tanks measuring 100 gallons each are labeled A, B, and C. 6 pipes are connected to the tanks and the rate of flow is as follows. Into A: pure water at four gallons per minute. From B to A: mixture at two gallons per minute. From A to B: mixture at six gallons per minute. From C to B: mixture at one gallon per minute. From B to C: mixture at 5 gallons per minute. Out of C: mixture at four gallons per minute.

    FIGURE 2.9.5 Mixing tanks in Problem 8

  3. Two very large tanks A and B are each partially filled with 100 gallons of brine. Initially, 100 pounds of salt is dissolved in the solution in tank A and 50 pounds of salt is dissolved in the solution in tank B. The system is closed in that the well-stirred liquid is pumped only between the tanks, as shown in FIGURE 2.9.6.

    (a) Use the information given in the figure to construct a mathematical model for the number of pounds of salt x1(t) and x2(t) at time t in tanks A and B, respectively.

    (b) Find a relationship between the variables x1(t) and x2(t) that holds at time t. Explain why this relationship makes intuitive sense. Use this relationship to help find the amount of salt in tank B at t = 30 min.

    Two tanks measuring 100 gallons each are labeled A and B. Two pipes are connected to the tanks and the rate of flow is as follows. From B to A: mixture at three gallons per minute. From A to B: mixture at two gallons per minute.

    FIGURE 2.9.6 Mixing tanks in Problem 9

  4. Three large tanks contain brine, as shown in FIGURE 2.9.7. Use the information in the figure to construct a mathematical model for the number of pounds of salt x1(t), x2(t), and x3(t) at time t in tanks A, B, and C, respectively. Without solving the system, predict limiting values of x1(t), x2(t), and x3(t) as t → .
    Three tanks are labeled A, B, and C. Their respective volumes are 200, 150, and 100 gallons. Four pipes are connected to the tanks and the rate of flow is as follows. Into A: pure water at 4 gallons per minute. From A to B: mixture at 4 gallons per minute. From B to C: mixture at 4 gallons per minute. Out of C: mixture at 4 gallons per minute.

    FIGURE 2.9.7 Mixing tanks in Problem 10

Predator–Prey Models

  1. Consider the Lotka–Volterra predator–prey model defined by

    where the populations x(t) (predators) and y(t) (prey) are measured in the thousands. Suppose x(0) = 6 and y(0) = 6. Use a numerical solver to graph x(t) and y(t). Use the graphs to approximate the time t > 0 when the two populations are first equal. Use the graphs to approximate the period of each population.

Competition Models

  1. Consider the competition model defined by

    where the populations x(t) and y(t) are measured in the thousands and t in years. Use a numerical solver to analyze the populations over a long period of time for each of the cases:

    (a) x(0) = 1.5, y(0) = 3.5

    (b) x(0) = 1, y(0) = 1

    (c) x(0) = 2, y(0) = 7

    (d) x(0) = 4.5, y(0) = 0.5

  2. Consider the competition model defined by

    where the populations x(t) and y(t) are measured in the thousands and t in years. Use a numerical solver to analyze the populations over a long period of time for each of the cases:

    (a) x(0) = 1, y(0) = 1

    (b) x(0) = 4, y(0) = 10

    (c) x(0) = 9, y(0) = 4

    (d) x(0) = 5.5, y(0) = 3.5

Networks

  1. Show that a system of differential equations that describes the currents i2(t) and i3(t) in the electrical network shown in FIGURE 2.9.8 is

    A circuit diagram. Inductor L is in series with 2 parallel paths. The first one consists of resistor R subscript 1. The second path consists of resistor R subscript 2 and capacitor C connected in series. Current i subscript 1 from source E flows through inductor L. It then branches into i subscript 2 and i subscript 3. i subscript 2 flows through the path with resistor R. i subscript 3 flows through the path with resistor R subscript 2 and capacitor C.

    FIGURE 2.9.8 Network in Problem 14

  2. Determine a system of first-order differential equations that describe the currents i2(t) and i3(t) in the electrical network shown in FIGURE 2.9.9.
    A circuit diagram. Resistor R subscript 1 is in series with 2 parallel paths. The first one consists of inductor L subscript 1 and resistor R subscript 2 in series. The second path consists of inductor L subscript 2 and resistor R subscript 3 connected in series. Current i subscript 1 from source E flows through resistor R subscript 1. It then branches into i subscript 2 and i subscript 3. i subscript 2 flows through the path with L subscript 1 and R subscript 2. i subscript 3 flows through the path with L subscript 2 and R subscript 3.

    FIGURE 2.9.9 Network in Problem 15

  3. Show that the linear system given in (18) describes the currents i1(t) and i2(t) in the network shown in Figure 2.9.4. [Hint: dq/dt = i3.]

Additional Mathematical Models

  1. SIR Model     A communicable disease is spread throughout a small community, with a fixed population of n people, by contact between infected individuals and people who are susceptible to the disease. Suppose initially that everyone is susceptible to the disease and that no one leaves the community while the epidemic is spreading. At time t, let s(t), i(t), and r(t) denote, in turn, the number of people in the community (measured in hundreds) who are susceptible to the disease but not yet infected with it, the number of people who are infected with the disease, and the number of people who have recovered from the disease. Explain why the system of differential equations

    where k1 (called the infection rate) and k2 (called the removal rate) are positive constants, is a reasonable mathematical model, commonly called a SIR model, for the spread of the epidemic throughout the community. Give plausible initial conditions associated with this system of equations. Show that the system implies that

    Why is this consistent with the assumptions?

  2. (a) In Problem 17 explain why it is sufficient to analyze only

    (b) Suppose k1 = 0.2, k2 = 0.7, and n = 10. Choose various values of i(0) = i0, 0 < i0 < 10. Use a numerical solver to determine what the model predicts about the epidemic in the two cases s0 > k2/k1 and s0 ≤ k2/k1. In the case of an epidemic, estimate the number of people who are eventually infected.