In the last few sections, we discussed first-order models of various systems and studied the types of interactions that could be modeled and the nature of the solutions of these models. Of the several indicated generalizations that could be made, this section will consider adding another state variable, so that the effects of two interacting variables can be used and studied. This will greatly increase the class of systems we can model and the class of solutions that result. In addition, a very interesting set of classical problems fall into this class with interesting solutions and interpretations.

To illustrate the general problem, consider a system that contains populations of two different types with distinctly different characteristics. Assume these two populations have a strong effect on each other, as well as being influenced differently by their environment, so that modeling them by a single total population would not yield useful results. We must, therefore, have two separate state variables to describe the systems, and this could perhaps be done in the following way.

Here the rate of change of population p1p1 is assumed to depend on both the populations p1p1 and p2p2; and likewise, the rate of change of p2p2 is assumed to depend on p1p1 and p2p2, but in perhaps a different way.

Many types of interactions could be considered. It might be that p1p1 and p2p2 compete for the same source of food or resources; it might be that p1p1 is a prey of the predator p2p2; or it could be that they both contribute to the welfare of the other. These assumptions would be implemented in the choice of f̲f̲ and g̲g̲ to describe the particular case. The best known classical models of these types were proposed by Lotka (1925) and Volterra (1926). Later, Gause (1934) did further experimental and interpretative work. Most of this type of work was done in population ecology.[6], [1].

Consider the particular for for the two-variable model to be

This might be a simple model of two competing populations, where aa and cc are the net rate of increase that would occur if the other population did not exist. The coefficients bb and dd model the negative effects of interaction on the rates of change as a measure of how often one encounters the other.

To simplify the mathematics, a change of variables will be made. Consider the rearrangement of Equation 3 into

Now let x=dcp1x=dcp1 and y=bap2y=bap2 then, Equation 5 becomes

Note that xx and yy are related to p1p1 and p2p2 by simple constant multipliers or scale factors, and therefore, the nature of the solution of Equation 7 is the same as Equation 3, but now there are only two parameters, aa and cc, to consider. In fact, by allowing a change of scale of the time variable, it is possible to reduce the number of parameters to one, but we will not do that. The problem of solving the coupled equation of Equation 7 or, more generally, of Equation 1 can be approached three ways. In some cases, an analytical equation for the solution can be found. This is always true if the equations are linear, but unfortunately, almost never true if they are nonlinear. Another approach was the phase plane where one solution is plotted as a function of the other, with time as an implicit variable. Very important characteristics of the solution can often be determined by phase plane methods without actually finding the solution. Finally, the equations can be numerically solved by simulation on a digital computer using Euler's method or some other more efficient algorithm.

The pair of equations in Equation 1 can be reduced to a single equation by eliminating the time variable tt. This can be done by simply dividing one by the other to give

The solution of this equation is examined in the p1,p2p1,p2 plane, which is called the phase plane.

As an example, consider the competition model in Equation 7 in the phase plane

Solutions in the phase plane are

The trajectory starting at AA gives the value of xx and yy with time tt, an implicit variable, indicated by the values shown. If a different initial mixture of populations had been assumed, e.g., BB, then a different trajectory would result. Indeed, any initial mixture is a point on the phase plane, and the trajectories indicate how they evolve in time. The more conventional time solutions are shown for the initial of AA by

and for BB by

Note the relation of the phase plane plots and the time plots. This particular problem will later be examined in greater detail.

One might wonder why this peculiar representation of the solutions is the form of one variable considered as a form of the other. This phase plane approach, although a bit unnatural at first, proves to be a very powerful tool. It is used by many in the literature [6] [10] [1] [2] and is a standard mathematical tool. [8] [9] It is worthwhile developing this concept before analyzing several physical systems.

Note that the phase plane contains all possible time plane plots for various mixtures. It can be shown that if the system has unique solutions, then the phase plane trajectories cannot cross. This means that a few key trajectories can be constructed which will make obvious what all other trajectorie will have to be. For example, in the above competition model, the initial mixture always determines who the eventual winner will be. Any initial mixture to the right of the line from the origin to CC results in yy increasing without bound and xx becoming extinct. Initial mixture to the left gives the opposite result.

There are several procedures that aid in the construction and interpretation of phase plane trajectories. There are special points on the plane known as equilibrium points or singular points that are important. If both dxdtdxdt and dydtdydt in (Reference) are zero, then xx and yy are constants and the system is in equilibrium. This means that at these points both the numerator and denominator of Equation 9 are zero. For the competition model of Equation 10, there is a singular point at x=0x=0 and y=0y=0, and another at x=1x=1 and y=1y=1. Singular points may be stable or unstable depending on whether small perturbations away from the point tend back to it or go away from it. Both points mentioned above are unstable.

A particular informative way of finding the singular or equilibrium points is to consider what are called partial equilibrium lines in the phase plane. The curve of all possible solutions of the equition

is called the partial equilibrium curve for population p1p1. This is understood by considering the first equation in Equation 1 alone. The equation Equation 11 implies dp1dt=0dp1dt=0, therefore, one side of this curve dp1dtdp1dt will be positive and on the other side it will be negative. If a particular f=0f=0 curve was given by

for any given fixed p2p2, p1p1 would move to the f=0f=0 partial equilibrium curve. This curve would, therefore, give the effects of p2p2 on the equilibrium values of p1p1. In other words, for a system controlled by the first equation of Equation 1 if p2p2 is given, the f=0f=0 curve will give the equilibrium value psub1psub1 approaches. .pp In fact, however, p2p2 is not fixed, but must obey the second equation in Equation 1. If this equation is examined separately, we have a second curve called the partial equilibrium curve for p2p2 given by

A similar analysis of this equation shows the effects of p1p1 on the equilibrium values of p2p2, and can be visualized by the following illustration of a g=0g=0 curve.

If these two curves are considered simultaneously, then not only are the singular points determined by the intersections, but the stability of the points and the nature and direction of the trajectories can be estimited by the signs of p ˙ 1 p ˙ 1 and p ˙ 2 p ˙ 2 in the various regions. For these illustrated curves of f=0f=0 and g=0g=0, we have This determines the singular point, and the directions show that it is stable. Applying this to the Lotka-Volterra competition model of Equation 7 for the partial equilibrium curves gives

or

and

or

In the phase plane, these are

Another tool that is very useful and is related to the preceding discussion is the method of isoclines.[8](Reference) Here we find curves in the phase plane where all the trajectories that cross that curve have the same slope. The partial equilibrium curves are two isoclines. The f=0f=0 curve implies from Equation 9 that the slope of all trajectories along that curve is zero. The slope of all trajectories along the g=0g=0 curve is infinite. If we find the isocline for a slope of mm, this is done from Equation 9 by setting

For the competition model with a=c=1a=c=1, we have

Solving for xx as a function of yy gives

The m=0m=0 isocline is y=1y=1. The m=∞m=∞ isocline is x=1x=1. The m=1m=1 isocline is

and m=-1m=-1 gives

In the phase plane the isocline looks like

Note how the isoclines aid one in sketching or visualizing the phase plane solution trajectories.

This should be enough detail on this approach to allow application to the various two-variable models that can be so interesting.

We will now return to the competition model of Equation 3 and examine it in more detail. Consider a situation where the uninhibited growth rate of population p1p1 is 10%. This implies a=0.1a=0.1 in Equation 3. Assume that the negative effects of p2p2 are such that 100 members of p2p2 cancel the positive effects of one member of p1p1. >From the first equation of Equation 3, we have

If a=0.1a=0.1, then b=0.001b=0.001. We also assume that p2p2 has the same self-growth rate, and p1p1 affects p2p2 in the same way that p2p2 affects p1p1. This gives c=0.1c=0.1 and d=0.001d=0.001. The model becomes

Using Euler's method to convert these differential equations to difference equations, we see that

These were programmed on a Tektronix 31 programmable calculator with an automatic plotter to give the phase plane output shown in Figure G. The trajectories were generated by running the simulation with various initial populations. For example, the lowest trajectory was run with an initial population of p1=25p1=25 and p2=50p2=50. The next one used p1=30p1=30 and p2=35p2=35.

If a different situation is considered where one population has a growth rate of 20% and the other 5%, but the interactions are still at a ratio of 100, the equations become

The solutions for this case are shown in Figure H. Here the results of the different rates are rather startling. The trajectory number 1 starts at p1=10p1=10 and p2=1p2=1, yet p2p2 overcomes p1p1. Trajectory number 2 starts at p1=16p1=16 and p2=1p2=1, and p2p2 still wins; but when the initial values are p1=17p1=17 and p2=1p2=1, trajectory number 3 shows p1p1 wins For p2=500p2=500 and p1=200p1=200 or 240, trajectories numbers 4 and 5 show p2p2 wins; but with p2=500p2=500 and p1=250p1=250 or 300, trajectories 6 and 7 show p1p1 wins. This exemplifies the very large difference a four-to-one growth rate ratio can make, and how critical the outcome depends on the initial values. It also illustrate the power of the phase plane in describing the model.

In the basic competition model described by Equation 3, and when normalized, described by Equation 7, we see that even if the interactive terms are very small, one population always grows without limit and the other becomes extinct. This describes a "survival of the fittest" model, but the unlimited growth and no possibility of coexistence seems unreasonable.

The next level of complication is the addition of a limit to growth in the same manner that the exponential was changed to a logistic. A crowding or self-competition term is added to the simple competition model. Consider now

Using the normalizing procedure that was used before on Equation 3 reduces the number of parameters from six to four:

where

Consider the partial equilibrium curves for this model.

On the phase plane, this becomes

It is obvious that the character of this system depends on the relative values of a,a,b,b,KK and LL, and indeed these are from rather different possible systems.

We will first consider the case illustrated above where both limiting factors are relatively small.

Note that as KK and LL approach zero, the system approaches the previously studied system. For this case, there are three possible equilibrium or singular points. There is an unstable point at the intersection of the two partial equilibrium curves, and a stable point at y=0y=0, x=aKx=aK, and another stable point at y=cLy=cL, x=0x=0. In this case, as before, one or the other population always wins, depending on initial conditions, and the remaining population dies to zero. There is now a limit reached by the winner and indeed, the time plot of the winning population looks very similar to a logistic. For example, for particular initial xx and yy, we have

The phase plane trajectories are illustrated for the normalized variables xx and yy in Figure J. The terms aa and cc are set equal to one, with KK and LL set equal to one-half. The winning population approaches 2 as its equilibrium value, and the loser becomes extinct.

The second case to consider has strong self-limiting factors relative to the interactive terms

The partial equilibrium curves in the phase plane are

Here the signs of the derivatives in the various regions of the phase plane show that there is only one stable equilibrium point at the intersection of the two curves. Here is a case of stable co-existence predicted for a competition model. The phase plane trajectories are shown in Figure K for the normalized variables where a=c=1a=c=1 and K=L=2K=L=2.

The third case is not symmetric. It allows one population to have a stronger self-limiting feature, and the other a stronger interactive term. This is given by

The partial equilibrium curves in the phase plane are

The equilibrium point at x=0x=0 and y=cLy=cL is the only stable point. For this case, yy always wins for any non-zero initial values, and xx always becomes extinct. Figure L illustrates the phase plane trajectories for the case where a=c=1a=c=1, K=2K=2, and L=12L=12.

The fourth case is similar to the third, but the roles of the two populations are reversed. The results are similar with xx always winning and approaching an equilibrium point of x=aKx=aK and y=0y=0. The phase plane trajectories look like Figure L with the axes reversed.

The use of these competition models can be very interesting in what they say about the effects of the various growth, interactive, and limiting parameters. Applications can be made in short time spans to competing populations in population ecology, or over longer time spans to biological evolution. There are many other possibilities of economic models or international models that could be pursued, but we now turn to a very different type of interaction to be modeled.

If the relation between two populations is not one of competition but one of one population preying on the other, a very different dynamic situation results. We first consider the simple Lotka-Volterra model where Equation 1 becomes

This represents a system where p1p1 is the population of the prey that has a growth rate aa when there are not predators. The parameter bb is the negative effect of predation on p1p1, and the product p1p2p1p2 models the frequency of encounter of the two. The population p2p2 is that of the predation who would die out at a rate of cc if there were no prey. The coefficient dd gives the positive effect of the prey on the predator, and again p1p2p1p2 models the frequency of encounter.

As was done before, if the populations are normalized by xdcp1xdcp1 and y=bap2y=bap2, then Equation 42 becomes

The phase plane equation is

Using the method of isoclines by finding the curves where the slope is constant shows a remarkable result. First, consider the partial equilibrium curves

On the phase plane, this is

The solution trajectories in the phase plane are shown in Figure M. It can be shown [8] that for any positive values of the parameters in Equation 42, the solutions in the phase plane are closed nested curves that enclose the singular point. The closed trajectories are called limit cycles, and they give rise to periodic or cyclic function when displayed as a function of time. The example used assumed an unlimited growth rate of the prey population p1p1 to be 5% per year. The death rate of the predator with no prey is set at 10% per year. The interactive terms are set to be equal to the self rates

when p1=100p1=100 and p2=200p2=200. This gives

There are several interesting features of the solutions. For any initial mixture of population, a limit cycle passes through it. The resulting oscillations have amplitude and frequency that depend on the starting condition, and oscillations neither grow or decay. Unfortunately, the use of Euler's method destroys the exact form of the solutions. Note that the trajectories did not exactly close on the left side of Figure M. They can be made to approximate the exact solution of Equation 51 by choosing TT very small – but that slows down the calculations and can sometimes cause other numerical errors.

Time plots of these solutions are shown in Figure N for initial values of p1=200p1=200 and p2=50p2=50. In Figure P the initial values are p1=50p1=50 and p2=200p2=200. Compare the initial values, maximum and minimum values, with the phase plane trajectories. Note that the period of the oscillation in Figure N is 91 years, and in Figure P, 110 years. The slight increase in amplitude of the oscillations is due to the Euler algorithm, not the model. The time interval TT was set at 0.2 years.

There are both interesting theoretical and practical aspects to this model. Serious error can occur when one of the populations is small. Minor variations which are assumed to average out with large numbers, do not. In many experimental verifications of this model, one of the populations will die out at a minimum rather than regenerate. Also, the model is rather sensitive to small errors. The addition of small terms to the basic equation Equation 42 causes great change in the character of the solution. This model has been studied in detail by population ecologists. [6] [7] [8] [10] [1]

We will next examine the effects on the simple predator-prey model of adding a crowding term as was done on the competition model and the logistic model. Consider the model

where the coefficients ee and ff describe the negative effects of crowding and competition within the population or perhaps cannibalism. These equations can be normalized as done before to the form

The effects of the added terms are rather dramatic. The partial equilibrium curves are shown in the following phase plane.

This assumes that aK1aK1. The singular points are denoted by a circle. The singular points at the origin and at x=aKx=aK, y=0y=0 are unstable, while the one at the intersection of the two curves is stable. The equations were programmed with a=c=L=1a=c=L=1 and K=0.5K=0.5. The trajectories in the phase plane are shown in Figure Q. Compare these results with the derivative signs and singular point locations found above. Note that if K=L=0K=L=0, the two partial equilibrium curves become vertical and horizontal, giving the same results as found earlier in Figure M. Solutions of this model are shown as a function of time in Figure R for initial values of x=1x=1 and y=2y=2, in Figure S for x=0.3x=0.3 and y=2y=2, and Figure T for x=2.5x=2.5 and y=0.3y=0.3.

Note the relations of these time curves to the phase plane trajectories in Figure Q. In all cases, there are "overshoots" and "undershoots" as the populations interact, but they finally settle down to a constant co-existence that is the same for any initial condition.

The model is changed by removing the limiting term on population yy by setting L=0L=0. This causes the yy partial equilibrium curve to become horizontal; the resulting phase plane trajectories are shown in Figure U. The results are similar to those in Figure Q, but there is more oscillation before the final equilibrium is reached. If a limiting factor is made large by setting L=4L=4, the phase plane trajectories of Figure V result, giving very little oscillation.

A rather different situation results if the parameters are such that ak1ak1. In this case, the intersection of the partial equilibrium is in the second quadrant which has no physical meaning. In the first quadrant where populations are positive, the equilibrium point at x=aKx=aK, y=0y=0 is the only stable one. This was programmed for ac=1ac=1 and K=2K=2. The phase plane trajectories are shown in Figure W, and the time solution for initial values of x=1.5x=1.5 and y=1y=1 in Figure X.

For these conditions, the predator dies out and the prey self-limits in a manner similar to the logistic.

By choosing more complicated interaction functions for the f(p1,p2)f(p1,p2) and g(p1,p2)g(p1,p2), it is possible to obtain other types of solutions. For the simple case with no limiting in Equation 42, the partial equilibrium curves were vertical and horizontal straight lines and the trajectories were closed. With limiting added in Equation 53, the curves remained straight, but were no longer ecessarily vertical or horizontal, and the solution trajectories were no longer closed, but would either spiral or smoothly move to an equilibrium point. Although not illustrated here, it is possible to use a model of limiting similar to (Reference) that will cause a single stable limit cycle to occur, that all trajectories starting outside of it would spiral in to it, and all starting inside of it would spiral out to it. This would give a steady-state oscillation as a time function. Perhaps this type of model could be used to explain some of the cyclic variations that occur in business and economics. Much more work could be done on both the mathematics and interpretation of this predator-prey type system, but we will move on to others now.

Here we assume a simple system consisting of a population yy that depends on, and consumes, a resource that cannot be replaced. The equations are somewhat similar to the predator-prey model, but the prey could grow and the resource here cannot. If xx is the amount of resource, and it is distributed in such a way that the consumption by yy is modeled by the product xyxy, the normalized equations are:

The partial equilibrium curves in the phase plane are

The phase plane trajectories are shown in Figure Y. This uses a=c=1a=c=1.

Note that the resource monotonically decreases while the population may or may not initially increase, but in any case, it ultimately dies out. An interesting result of this model is that there is some resource left after the population is gone. This is caused by the assumption that the distribution is such that consumption is governed by the product xyxy. If it is assumed that the resource is easily accessible, a better resource model might be

The nature of the solutions of this system is left to the reader.

We will now move into rather different systems to see how models might be applied. The history of application of dynamic models to problems such as national armament is fairly new, but perhaps older than most realize. [5][4]

A simple linear model is the following

where the state variables xx and yy are measures of the arms level of two nations. The coefficients aa and cc are measures of confidence or expense that cause a decrease in military expenditures; bb and dd are the effects of the opponent's arms level on ones military build-up. The constants ff and gg represent the minimum level that would be maintained even if the opposition disarmed.

There are two general cases possible. If the situation is such that acac>bdbd, then the partial equilibrium curves look like

The signs of the derivative show that the singular point is stable. This states that the arms race stops at a stable level for this case. If, on the other hand, the conditions are such that bdbd>acac, the curves look like.

Here, there is no stable point, and the armament levels of both nations increase without limit.

This model is linear so that analytical expressions for the solution can be found, and computer simulation is unnecessary. On the other hand, the model is too simple to be realistic, and a more reasonable one would be nonlinear. Again, while these models can be interesting, they leave out too many other state variables to be used for more than gaining insights.

While it is certainly easier to model and measure quantitative variables such as population, food, or money, it is also possible to apply dynamic modeling techniques to more subjective variables involved with attitudes and feelings. These variables must be quantified in some way that is obviously going to be somewhat subjective. Even though this process is difficult and subject to challenge, it must be done if more complete models of social systems are to be developed. This becomes apparent when, in trying to chose the state variables for a system, it is necessary to know how a group of people feel about another variable to predict their actions. One accepted example is the practice of assigning a monetary value to the good will of a company.

As an example, we will consider the dynamics of feelings between two populations in terms of their friendliness or hostility. [2] Let xx be the measure of friendliness of population p1p1 toward population p2p2, and yy the friendliness of p2p2 toward p1p1. Negative friendliness is considered hostility. The equations are

The determination and interpretation of ff and gg is a bit more difficult here. Recall from Section B the definition of partial equilibrium curves. The f(x,y)=0f(x,y)=0 curve in the phase plane gives the equilibrium values of xx for a fixed yy. For this model, f(x,y)=0f(x,y)=0 gives the degree of friendliness xx that will be approached by p1p1 for an independently set amount of friendliness yy of p2p2 toward p1p1. Consider the following case:

Here, if yy is neutral, then xx become neutral. If yy is friendly, then xx is friendly. As yy becomes more friendly, xx increases to a point and then levels off at a maximum amount of friendliness. As yy becomes hostile, xx responds likewise and finally levels off at a maximum amount of hostility.

Now consider the g(x,y)=0g(x,y)=0 curve which is the p2p2 response to p1p1. This is shown by

Population p2p2 is naturally more friendly than p1p1. It is friendly even if p1p1 is neutral, as is shown by point AA. Only after p1p1 becomes fairly hostile does p2p2 begin to return with hostility as shown by point BB.

Given these relations and considering the signs of the derivatives, it is seen that the singular point at CC is a stable equilibrium point.

Consider a different set of characteristics where a population would initially return a large amount of hostility a hostile opponent, but after a certain level, would submit or surrender by having a reactionary response to a very hostile situation. This might be described by

Another characteristic is to have almost no response up to a certain level, then to react suddenly as described by:

Many interesting models can be posed and the resulting solutions examined in the phase plane. In some cases, the results are insensitive in the sense that small changes in the partial equilibrium cause only little change in the solutions. Other cases are very sensitive.

A very important aspect of this approach to modeling is the dynamic description. When the trajectories or time solutions are found, not only are the equilibrium points found, but how they are reached is predicted. In some cases, the effects of a .... more important than the final value. Also, sometimes the time necessary to achieve a certain condition is as important as the condition itself.

Reading, in this day and time, the 1798 essay by T.R. Malthus, one is struck by both the insight and the naite of this very influential statement. [4] Malthus saw the possible dire consequences of an increasing population in a finite environment. His predictions of doom were based on the assumption that the world population increases according to a geometric sequence, while the food increases according to an arithmetic sequence. Stated in our terms, he assumed population increases exponentially and food linearly. The fact that world food production has more or less kept up with the population has lead critics to discount Malthus and his followers as irrelevant doomsday prophets. In fact, some feel this pessimistic view is not irrelevant, it is dangerous. The strong critics, known as technological optimists, assume that the factors that have prevented Malthus' predictions from coming true so far will continue to so do. [3] In fact, they claim that growth is not the problem, it is the solution. Growth has given us the highest material standard of living by our abilities of growing faster than our problems.

Our purpose in this treatise is not to take sides, but to suggest a different way of looking at our situation that will give more understanding and insight. Malthus based his prediction on observed and assumed growth patterns of population and food. If one looks at the underlying models that might support these assumptions, the population might come to be modeled by an equation similar to (8). One is hard-pressed to explain his assumed food growth, and that has indeed been the flaw in his assumptions.

The technological optimists likewise have implied models to obtain their predictions. The costs of continued economic and technical growth must be considered in any realistic model. A bit of reflection shows that the question is not whether to use a model or not, but to determine what kind of model to use and whether it will be examined and debated explicitly.

Very interesting results have been obtained in the preceding sections on applying two-variable models to various systems of interest to us. It is a valuable exercise to try to better model world population and food than Malthus did. Very soon one discovers (or should discover) that the systems are too complicated to be described by any two-state variables. Unfortunately, when using higher order models, many of our analytical methods no longer work. We lose the powerful tool of the phase plan. Second-order systems are useful to gain insight into simple systems with relatively simple interactions, but now we will have to rely primarily on computer simulation to give solutions of the resulting third and higher order models.