Ex 6.1 - Simple Surface Water Model

1. In Stella you can clamp your state variables to make sure that they never become negative. For example in this model the Surface_Water is non-negative. Note that the "non-negative" option is checked in the variable definition box. Since by default the variables get clamped in Stella this may be sometimes somewhat confusing and may hide some of the errors, when the variable is actually negative, but you do not see it. It is good practice to make sure that your processes (flows) are described properly, and do not deplete state variables beyond levels that are intended. Uncheck the non-negativity in this model and see what realy happens to the Surface_Water. Redefine the flows in the model to make sure that Surface_Water does not go negative. 


Obviously soil saturation is going to be a major determinant of the rate of surface water infiltration, but since we don't have a soil component in our model (yet, I hope), I just made infiltration a function of surface water: infiltration = .5 * surface_water

 
2. Let us supposed that the whole area got paved. How do we describe this in the model? What happens to Surface_Water? Does the result look plausible? Are there any other processes that we may be missing?

Assuming that the pavement is totally impermeable, surface water simply accumulates as the integral of rainfall. Obviously this is unrealistic, surface water will leave the site by flowing to lower elevations.

Ex. 5.3 - Trophic Chains

How does the model dynamics change if we consider the mortality process in the last model with the uptake of resource limited by the biomass in the first trophic level (N = u₀T₁)? Consider even and odd trophic chains and try to make some generaliztions. 

Well, the models' behaviors are interesting and quite divergent. Generalizations are going to be tough. Let's look at the three-trophic-levels model first. With death coefficient (d) = .05, the top predator dies off rapidly and the other two assume a periodicity that is reminiscent of the earlier predator-pray model.

Three levels, d = 0.05

If we decrease the mortality to d = .001, we get very different behavior. The first and third tropic levels decline, while the middle level enters a stable oscillation. I was curious how the system would collapse as levels 1 and 3 became extinct, so I extended the model to run to t = 1000, but the dynamics are unchanged. Levels 1 and 3 continue to decline, but are never eliminated, so level 1 continues to act as a conduit of "N" for level two.

Three levels, d = 0.001

Turning now to the four-level model, with d = .05, the overall behavior of the system is rather similar to the three-level model with d = .05 (the first set of graphs). There are some quantitative differences, but we see the same oscillations for the two lowest trophic levels (blue and red), and the highest (green) rapidly collapses.

4-levels, d = 0.05

In contrast, the four-level system with d = 0.001 exhibits very different behavior than the three-level system with d = 0.001. Whereas levels 1 and 3 declined in the three-level system, all four levels are stable in this system. Levels 1 and 3 oscillate roughly in-phase (lower-left scatter plot, note roughly direct proportionality), while levels 2 and 4 oscillate together, out-of-phase relative to levels 1 and 3.

4-levels, d = 0.001

Generalizations here are tough, but I'll try. One is that alternate trophic levels seem to move together, regardless of other details of the system, and adjacent trophic levels often move opposite each other. These are intuitive when we think in terms of top-down control: If there's a boom in coyotes, the will be a crash in rabbits, which will cause a boom in lettuce. In contrast, if we think in terms of bottom-up control, we wouldn't expect such behavior: a boom in lettuce will yield a boom in rabbits which will yield a boom in coyotes. Perhaps we can conclude that, at least for the system as we've formulated it, predator-based population controls dominate resource-based population controls. I wonder if I can turn down predation and turn up T1's resource use to make it not so... Indeed, just turning up u0, the first trophic level's source breaks the T1-T3 correlation. It also makes the system much flashier. Compare the above lower-left scatter plot, where when T1 is large, T3 is large and vise-versa, with the following lower-left scatter plot, where T1 and T3 appear to be inversely related. The only difference between them is the increased flow to the first trophic level (from 0.1 to 0.5).

And indeed, turning u0 up even further (to 1.0) strengthens the inverse relationship between T1 and T3. So perhaps we can conclude that when reourses at the bottom of the trophic web are abundant, they tend to govern the system, and when they are scarce, predation governs the system. That seems intuitive enough.

Ex. 5.2 - Predator-prey model sensitivity analysis

Run sensitivity analysis for this model with respect to the Half-saturation parameter SS. What do you observe when SS varies in the [0.5, 1.5] interval? How can you explain that?

Run 1 in red is SS = 0.5; run 5 in purple is SS = 1.5

As SS increases, the system converges toward the stable, non-trivial equilibrium more rapidly; that is, increasing SS dampens the extremes of rabbit and wolf populations as they approach their equilibrium values. The equilibrium populations are unaffected. This can be explained on several levels.

Mathematically, since SS is in the denominator of the predation function, as SS increases, predation slows, which effectively stabilizes the system.

Thinking about the system ecologically, if wolves consume rabbits very rapidly—that is, the slope of V vs. x is steep, which is a consequence of lower values of SS—then an increase in rabbit populations will produce a rapid upward swing in wolf population, because predation will increase their numbers before wolf mortality stabilizes their population. As a result of the increased number of wolves, the rabbit population will plummet, which will lead to a wolf die-off, and the cycle will repeat. On the other hand, when SS is large, predation is slow, so the wolf population reacts more gradually to fluctuations in the rabbit population, which effectively stabilizes both populations.

Exercise 5.1 - Modeling preditor-prey populations

1. Can you think of any examples of other systems that demonstrate the kind of behavior that we found in the predator-prey model?

Parasite-host dynamics are likely similar, where parasite population increases until the host’s capacity is reduced, at which point the parasite population declines until it finds a new host or the original host recovers to the point where the parasite can again expand.

2. In some predator-prey systems the prey can take refuge to hide from predators and avoid being consumed. Usually there is only a certain fixed number of individuals that the refuge can house. When the population of prey is large there is not enough refuge for all and the prey that could not find a place to hide gets consumed like in the standard predator-prey formalism. However when there are just a few preys their consumption slows down because they can find enough refuge places to hide. Build a predator-prey model with refuge and describe the dynamics that you observe. What equations did you modify and how? How can you explain the effect of refuge on the overall system dynamics?

To accomplish this, I made predation an exponential function:

V = y*a*e^b*x

Where V is predation, y is # wolves, x is # rabbits, and a and b are parameters that represent, respectively, the magnitude of the refuge available to the rabbits and how easily rabbits beyond the refuge are taken by wolves.

In general, this has the effect of stabilizing the system’s equilibrium. Effectively, the rabbit population is limited to the size of the refuge and the wolf population consequently limited by the size of the rabbit population. This causality can be seen reflected in the fact that the wolves’ curve reaches its asymptote slightly later than the rabbits’ curve.

Exercise 2.3 - Wizzles and Woozles

This is exercise 3.3 from the text, which is 2.2 online, and 2.3 in my nomenclature.

1. Consider a population of wizzles, which unlike woozles can multiply only when their numbers are larger than a certain minimal value a. If there are less byukies than a, they cannot find a partner for mating and the population dies off. The carrying capacity of the island where wizzles live is also A. Obviously A > a. What model would describe the population dynamics of wizzles on this island?

2. What are the equilibria in the byukies model? Are they stable? 

The solution that comes to mind isn't very satisfying, but it works. I've added an IF statement to the growth rate, such that if the A > a, the growth rate is positive, and if A < a, the growth rate is negative. The first image shows the model still works when the Ao > a; the second image shows the stable equilibrium at zero if the Ao < a.

Exercise 2.2 - Modeling a Toilet

A note on my nomenclature: This is exercise 3.2 from the text; it also happens to be 3.2 online, but online it's part of the art of modeling chapter, and in the text it's part of the math chapter. Following the online order, I'm calling the math chapter ch. 2, but since I'm working from the text for the math chapter (and only the math chapter) I'm calling this ex. 2.2. Whatever.

The model, as I downloaded it, needed to be broken before it could be fixed. First, I deleted the "/dt" from the outflow equation. Second, I changed the flush's if statement to be "if Use>0.995..." because I wasn't getting many flushes with the 0.999 threshold. Third, I deleted the seed from the random function ["RANDOM(0,1,14)" -> "RANDOM(0,1)"] in the parameter "use" to get different random flushings on each run. Doing that gives the following with dt = 0.5:

Q: Can you reformulate the model in such a way that the stock would be drained out properly under any computational method and any time step DT? As you may have noticed, when we switch to other methods, the outflow starts to deplete the stock, but then, apparently, does not have enough time to take it all out. How can this be fixed?

To get the tank to empty no matter what the dt, I just cranked up the outflow. Since tank is defined as a non-negative reservoir, this has no consequences. I made the outflow 1000 x tank, which should work up to a dt of 0.001. The result is:

Math! Simluation Modeling Exercises 2.1

1. Consider a population [in] which cells divide into 3 every time step. What mathematical model will describe this population? 

x(n) = 3x(n - 1)
= 3[3x(n-2)]
= 3{3[3x(n - 3)]} ...
= a3^n

2. Formulate the model with triple growth in terms of a differential equation and solve it. Is the difference between the discrete and continuous models larger or smaller than in the double growth model?

 x(t + dt) = x(t) + 2x(t)dt
[x(t + dt) - x(t)] / dt = 2x(t)
dx/dt = 2x
dx/x = 2dt ... Integrating gives:
ln x = 2t + C
x(t) = ae^(2t)

Now that we have the continuous and discrete models (ie., continuously compounded vs. compounded once per "t"), I plugged the equations into Excel to see how different they are. And they're hugely different! And the difference is much greater for the tripling equations: there's a million-fold difference between the two computation methods by the 16th generation!

time	2x discrete	2x cont.	3x discrete	3x cont.	2x ratio	3x ratio
1	2.00E+00	2.72E+00	3.00E+00	7.39E+00	1.36E+00	2.46E+00
2	4.00E+00	7.39E+00	9.00E+00	5.46E+01	1.85E+00	6.07E+00
3	8.00E+00	2.01E+01	2.70E+01	4.03E+02	2.51E+00	1.49E+01
4	1.60E+01	5.46E+01	8.10E+01	2.98E+03	3.41E+00	3.68E+01
5	3.20E+01	1.48E+02	2.43E+02	2.20E+04	4.64E+00	9.06E+01
6	6.40E+01	4.03E+02	7.29E+02	1.63E+05	6.30E+00	2.23E+02
7	1.28E+02	1.10E+03	2.19E+03	1.20E+06	8.57E+00	5.50E+02
8	2.56E+02	2.98E+03	6.56E+03	8.89E+06	1.16E+01	1.35E+03
9	5.12E+02	8.10E+03	1.97E+04	6.57E+07	1.58E+01	3.34E+03
10	1.02E+03	2.20E+04	5.90E+04	4.85E+08	2.15E+01	8.22E+03
11	2.05E+03	5.99E+04	1.77E+05	3.58E+09	2.92E+01	2.02E+04
12	4.10E+03	1.63E+05	5.31E+05	2.65E+10	3.97E+01	4.98E+04
13	8.19E+03	4.42E+05	1.59E+06	1.96E+11	5.40E+01	1.23E+05
14	1.64E+04	1.20E+06	4.78E+06	1.45E+12	7.34E+01	3.02E+05
15	3.28E+04	3.27E+06	1.43E+07	1.07E+13	9.98E+01	7.45E+05
16	6.55E+04	8.89E+06	4.30E+07	7.90E+13	1.36E+02	1.83E+06
17	1.31E+05	2.42E+07	1.29E+08	5.83E+14	1.84E+02	4.52E+06
18	2.62E+05	6.57E+07	3.87E+08	4.31E+15	2.50E+02	1.11E+07
19	5.24E+05	1.78E+08	1.16E+09	3.19E+16	3.40E+02	2.74E+07
20	1.05E+06	4.85E+08	3.49E+09	2.35E+17	4.63E+02	6.75E+07

3. Suppose you have $10,000 on your savings bank account with an annual interest rate of 3%. Build a model to calculate your interest earnings in 5 years. What will be the difference in your earnings if the interest is calculated monthly instead of annually? How much does the bank make on your account by calculating the interest monthly, instead of doing it continuously?

Continuous interest accumulates by the formula:

A = Ao x e^(rt)

Where A is the final amount, Ao is the initial amount, r is the interest rate, and t is time elapsed.

So after five years:

A5 = $10,000 x e^(0.03 x 5) = $11,618.34

If the interest is compounded in discrete steps, the operative formula is:

  (thanks, Wikipedia!)

where n is the number of compounding periods per year, in this case 12, so:

A5 = $10,000 x (1 + 0.03/12)^(12 x 5) = $11,616.17

So the bank pockets $2.17 every five years by compounding interest monthly rather than continuously.

Is this discrepancy so much smaller than for the doubling and tripling functions because of the low interest rate? If, instead of a savings account, this were a delinquent credit account that accumulated interest at 28%, after five years, the $10k would be:

Monthly compounding: $39.904.99
Continuous compounding: $40,552.00

Indeed, the difference after five years at 28% (~$550) is about 5% of the initial amount, where the difference after five years at 3% (~$2) is about 0.02% of the initial amount. Since these are exponential functions, as more time elapses, the difference will become more pronounced, as it did in the doubling and tripling exercises.

So the take home message from these exercises seems to be: the greater the rate of change, the more important the time increment.