The graph above and the table below show the number of lynx trapped each year in the Mackenzie River District of Canada from 1821 through 1934. These records give some measure of the lynx population.
The Number of Lynx Trapped in the Mackenzie River District
Year Count | Year Count | Year Count | Year Count
----------------|---------------|---------------|--------------
1821 269 | 1850 361 | 1879 201 | 1908 345
1822 321 | 1851 377 | 1880 229 | 1909 382
1823 585 | 1852 225 | 1881 469 | 1910 808
1824 871 | 1853 360 | 1882 736 | 1911 1388
1825 1475 | 1854 731 | 1883 2042 | 1912 2713
1826 2821 | 1855 1638 | 1884 2811 | 1913 3800
1827 3928 | 1856 2725 | 1885 4431 | 1914 3091
1828 5943 | 1857 2871 | 1886 2511 | 1915 2985
1829 4950 | 1858 2119 | 1887 389 | 1916 3790
1830 2577 | 1859 684 | 1888 73 | 1917 674
1831 523 | 1860 299 | 1889 39 | 1918 81
1832 98 | 1861 236 | 1890 49 | 1919 80
1833 184 | 1862 245 | 1891 59 | 1920 108
1834 279 | 1863 552 | 1892 188 | 1921 229
1835 409 | 1864 1623 | 1893 377 | 1922 399
1836 2285 | 1865 3311 | 1894 1292 | 1923 1132
1837 2685 | 1866 6721 | 1895 4031 | 1924 2432
1838 3409 | 1867 4254 | 1896 3495 | 1925 3574
1839 1824 | 1868 687 | 1897 587 | 1926 2935
1840 409 | 1869 255 | 1898 105 | 1927 1537
1841 151 | 1870 473 | 1899 153 | 1928 529
1842 45 | 1871 358 | 1900 387 | 1929 485
1843 68 | 1872 784 | 1901 758 | 1930 662
1844 213 | 1873 1594 | 1902 1307 | 1931 1000
1845 546 | 1874 1676 | 1903 3465 | 1932 1590
1846 1033 | 1875 2251 | 1904 6991 | 1933 2657
1847 2129 | 1876 1426 | 1905 6313 | 1934 3396
1848 2536 | 1877 756 | 1906 3794 |
The interesting thing about this data is its almost cyclic behavior -- with
periods of high lynx population alternating with periods of lower lynx
population. This kind of behavior has been observed in other natural
habitats and in habitats in the laboratory shared by two species, one of which
is a predator and the other of which is its prey.
This will be a very wide-ranging module. We begin by looking at one model
of predator-prey interaction and we see that the mathematical techniques we
have developed so far are not sufficiently powerful to enable us to determine
the longterm behavior of these models. This is an important lesson for
students at any level. The real problems in which we are most interested
are often difficult and complex. In a particular class we may be able to
make some progress on these problems but we often cannot solve them
completely. We may learn more about them in subsequent classes or we may
have to wait until new mathematics is developed. You or some of your students
may help find that new mathematics.
In this situation the mathematical difficulties cannot be ignored -- in fact,
they are a clue that the model we are using has some problems. This is
a common situation. We are using numerical techniques to estimate solutions
to systems of differential equations (that is, continuous dynamical systems).
These techniques do not find perfect solutions. To someone interested
primarily in mathematics this situation is distressing. But to someone
interested in the study of populations this situation is enlightening.
Real populations in real habitats are subject to "errors" in the same way
that numerical methods are subject to errors. The fact that the errors
in numerical methods cause problems is a clue to the fact that in real
habitats we are likely to see problems as well.
We will look at some simpler models that produce cyclic behavior and we will
see that these simpler models can be used to describe simple spring-and-mass
systems. Our intuition about these simple physical systems can help us
understand the behavior of ecological systems. This is one of the wonderful
things about mathematics. It often allows us to make analogies -- sometimes very
precise analogies -- between seemingly different situations.
A Model of Predator-Prey Interaction
We begin by looking at the following model of one kind of predator-prey
interaction.
q' = (1 - 0.005 p) q
In the first equation, the factor (0.005 q - 1) determines the relative growth rate of the species p. This species is the predator and in this model it relies very heavily on the prey population represented by the letter q for food. In fact, in the absence of the prey the predator population decreases because p' is negative. It is only when q > 200 that p' becomes positive and the predator population increases.
In the second equation, the factor (1 - 0.005 p) determines the relative growth rate of the species q. This species is the prey. In the absence of p it would grow very rapidly but when p is above 200 the prey population drops. In the absence of p this is an exponential model for the population growth of q and since the relative growth rate of q is positive this model cannot be very realistic. It may not be a bad model, however, in the presence of p.
Draw a rough direction field or lake diagram similar to the drawings you made in the last two modules for this predator-prey model. Based on this diagram what possible long-term behavior might occur?
The emphasis in the question above is the bold-faced words -- possible and might. You won't be able to determine what long-term behavior will occur, just the possibilities. Answer the question above as completely as you can before following the link below to find out the answer.
Click here to open a window with a Java applet that we will use to explore this model further. Arrange these two windows as usual -- so that they overlap and you can move back-and-forth between them by clicking on the exposed portion of the inactive window to make it active.
This applet simulates the predator-prey model we've been studying. Begin by clicking on the graph at the point (300, 300) to establish this point as the initial condition. You should see a black dot appear, marking this point. Now click the red button in the middle of the bottom row of buttons. This will cause the applet to simulate the model for a period of time using a fast but not terribly accurate numerical method. Notice that as the "cork" is carried by the current, its path appears to spiral out. If this simulation is accurate then the populations of p and q will oscillate with larger and larger swings.
To see what would happen if this model continued, click the brown button (the fourth button from the left). The applet will continue the simulation for an additional period of time.
Now click the grey button at the right to clean the graph. Then click the blue button at the left. This will run the same simulation using a slower and more accurate numerical method. Notice that the cork still appears to spiral out but the spirals do not expand as quickly as before.
In fact, the exact solutions of this continuous dynamical system (that is, this system of differential equations) are closed loops. We have not discussed mathematical techniques that will enable us to prove this fact. Notice, however, that it does require proof. Our experimentation is inconclusive. If this model was a good model then we would see exactly the kind of behavior we saw in lynx data from the MacKenzie River District with the population oscillating between periods of relatively low population and periods with relatively high population.
This model, however, has one unpleasant feature that is revealed by our mathematical difficulties. The most important adjective for real world populations and equilibriums is the adjective stable. You may want to click here to open a window with some Java applets from the last module that we will use to illustrate the idea of stability for equilibrium points. Arrange these two windows as usual -- so that they overlap and you can move back-and-forth between them by clicking on the exposed portion of the inactive window to make it active.
The equilibrium point in the middle of the graph for competitive interaction is a good example of a stable equilibrium point. Notice if you click anywhere near this equilibrium point the long term behavior approaches the equilibrium. This is an important property for equilibrium points in the real world, especially for equilibrium points in biological systems. Real world systems are subject to all sorts of noise -- events, like hurricanes or droughts -- that are not part of our models but that can have an impact in the real world. These events can knock the system out of equilibrium but if the equilibrium is stable then if it is not knocked too far away from the equilibrium it will return to the equilibrium.
The equilibrium point in the middle of the graph for aggressive interaction is a good example of an equilibrium point that is not stable. If the initial populations of the two species are exactly equal then the long term behavior will approach this equilibrium but if anything happens that gives either species the slightest advantage that species will eliminate the other one. Try clicking close to this equilibrium. Unless your click is very precise the long term behavior will move away from the equilibrium.
When you are done with this applet, close its window.
| The graph at the right shows a different kind of "equilibrium" -- a dynamic equilibrium. Click anyplace on the graph and watch what happens. The "cork" moves around and around in a circle repeating the same path over and over again. This circle is an equilibrium cycle. After the action stops, click some place outside the equilibrium cycle and observe what happens. This time the "cork" swirls around in a spiral getting closer and closer to the equilibrium cycle. After the action stops, click someplace inside the equilibrium cycle and watch what happens. This time the "cork" swirls in an expanding spiral getting closer and closer to the equilibrium cycle. This is an example of a stable equilibrium cycle. Even if the system is knocked away from this equilibrium behavior, it returns to the same behavior. |
Spring and Mass Systems
Click here to open a window with a Java applet that we will use to explore a new physical situation that can help us understand the predator-prey model above. Arrange these two windows as usual -- so that they overlap and you can move back-and-forth between them by clicking on the exposed portion of the inactive window to make it active.
At the top of this applet is a sketch showing a horizontal grey track and a vertical grey wall. A spring is attached to the grey wall and to a brown cart on wheels that is free to move back-and-forth along the track. The cart is way over at the right of the track and the spring is stretched. It is pulling the cart back toward the left. Click on the brown graph and observe what happens.
The spring pulls the cart to the left and the cart begins to accelerate. After a while the spring is compressed and begins to push the cart toward the right. This slows the cart down and after a while the cart stops for an infinitely brief moment. Then it begins to move toward the right. After a while the spring is stretched and begins to pull the cart toward the left again. This slows the cart down until it stops for an infinitely brief moment at its initial position. This cycle repeats itself because there is no friction.
This system can be described by the differential equation
where k is a positive constant that is determined by the mass of the cart and the strength of the spring. The variable x measures the position of the cart with x = 0 representing the point at which the spring is relaxed and exerts no force. The minus sign in this equation represents the fact that when the end of the spring is to the right of the spring's relaxed position the spring is stretched and pulls the cart back toward the left. Similarly, when x is negative the spring is compressed and pushes the cart in the positive direction toward the right.
If the system has some friction then this equation is replaced by the equation
where p is a positive constant determined by the mass of the cart and how much friction there is in the system. The minus sign by p in this equation represents the fact that the force of friction acts in the opposite direction of the velocity.
By introducing a new variable v for velocity, we can rewrite this second order differential equation as a system of two first order differential equations -- that is, as the kind of two-dimensional dynamical system we have been studying.
v' = -k x - p v
The two graphs in the applet track the behavior of this physical system using both perspectives above. The green graph looks at the two variables x and v. The position, x, of the spring is measured along the horizontal axis. The graph is set up so that the arrow on the cart points down at the corresponding value of x. The relaxed position for the spring is at the center of the green graph. The vertical axis of the green graph is used for the velocity, v. The center of the graph represents v = 0. Click once again on the brown graph. The cart starts at the right edge of the graph with velocity zero. Notice as the cart moves to the left, its velocity is negative. At first it is going faster and faster until it passes the relaxed position. Then it begins to slow up and pauses infinitely briefly at the extreme left of its travels.
The brown graph shows position (in black) and velocity (in red) as a function of time. Notice how the green graph and the brown graph change as the cart moves back-and-forth.
You can experiment with different initial positions for the cart by clicking along the track above the green graph. Notice each time you click the position of the cart changes and the new initial condition is marked on the green graph. You can change both the initial position of the cart and its initial velocity by clicking at the appropriate point on the green graph. If you click on the track then the initial velocity is zero, as if you held the cart at that position and then released it. If you click on the green graph the initial velocity might be nonzero, as if you pushed or pulled the cart as you let it go. Try various different initial conditions.
Because the system has no friction, the Law of Conservation of Energy implies that the system will continue to oscillate forever with the same amplitude. In practice, systems like this always have some friction. Click on the scale at the far right at about the middle of the blue part of the graph. This will introduce some friction. Click on the brown graph to start the simulation and observe the results. Notice that this time the oscillations are damped and their amplitude decreases as friction removes energy from the system in the form of heat. Experiment with different amounts of friction and different initial conditions. As long as the friction is set by clicking in the blue part of the scale at the extreme right it is ordinary friction. You can also experiment with antifriction or negative friction by clicking in the pink portion of the friction scale. Antifriction is a fictional concept. Mathematically, it is just a change of sign for the constant p but physically it would represent some mechanism that would add energy to the system. What happens with antifriction?
The remainder of this module studies the numerical techniques we are using to estimate solutions for these kinds of systems but this is a very rich subject for study in both mathematics and physics. The exercises below ask you some questions that you might use to emphasize the connections between mathematics and physics.
You may want to look at the help modules for force and acceleration or for work and energy as you work through these questions.
The best way to write the original second order differential equation describing a spring-and-mass system is shown on the left hand side below. The corresponding system of two first order differential equations is shown on the right hand side.
The constant m is the mass of the cart and for purposes of this series of questions will be measured in grams. The variable x will be measured in centimeters. The constant k is called the spring constant and is measured in units of force per centimeter because the product -kx is the force produced by the spring when the cart's position is x. The constant p is called the coefficient of friction and is measured in units of force per (centimeters per second) because the product -pv is the force produced by friction when the velocity of the cart is v.
Euler's Method for Estimating Solutions of Initial Value Problems
Now we want to look one way of estimating solutions to initial value problems like the problems we have been studying. We will use the initial value problem
as an example. This is a good example to use because we understand it completely. It represents a spring-and-mass system with no friction and m = k = 1. From your work above you know that solutions of this dynamical system are circles and that the solution of the initial value problem above is a circle of radius 1.
|
It is also easy to draw the direction field for this dynamical system. See the
figure at the right. The dot marks the point (0.6, 0.8). At this
point x' = 0.8 and v' = -0.6. The direction arrow based at the dot
goes 0.8 units to the right and 0.6 units down. Notice that this arrow makes
a right angle with the line from the origin to the dot. There are two ways
to see this. This arrow is tangent to the solution through the point (0.6, 0.8)
and since the solution curves are circles, tangents are perpendicular to radii.
Another way to see this is by noticing that at any point
the direction arrow points in the direction and thus the slope of the direction arrow is the negative reciprocal of the slope of the line from the origin to the point. This implies they are perpendicular. |
|
We use Euler's Method to estimate solutions of initial value problems of the form
on an interval of the form [0, T]. The basic idea is simple. We start at the initial point (a, b) and we compute the "current" at that point
This current tells us in which direction and how fast the "cork" is moving. We can simulate this motion numerically by moving to the point
but there is a problem. If we move too far in this direction we may move quite far away from the real solution because we are moving along the tangent line rather than along the actual curve. Thus, we want to take a relatively small step. We accomplish this numerically by moving to the point
where h is a relatively small number. At this point we repeat the same procedure.
We describe this procedure symbolically as follows. We begin by dividing the time interval [0, T] into n pieces and we let h = T/n. Then we compute
for i = 0, 1, 2, ... (n - 1). The table below shows the results of applying Euler's method to our sample initial value problem on the time interval [0, 6.2832] with 100 steps of length h = 0.062832. You should verify the first few of these computations by hand.
old x old v new x new v
1.0000 0.0000 1.0000 -0.0628
1.0000 -0.0628 0.9961 -0.1257
0.9961 -0.1257 0.9882 -0.1882
0.9882 -0.1882 0.9763 -0.2503
0.9763 -0.2503 0.9606 -0.3117
0.9606 -0.3117 0.9410 -0.3720
0.9410 -0.3720 0.9176 -0.4312
0.9176 -0.4312 0.8905 -0.4888
0.8905 -0.4888 0.8598 -0.5448
0.8598 -0.5448 0.8256 -0.5988
0.8256 -0.5988 0.7880 -0.6507
0.7880 -0.6507 0.7471 -0.7002
0.7471 -0.7002 0.7031 -0.7471
0.7031 -0.7471 0.6562 -0.7913
0.6562 -0.7913 0.6064 -0.8325
0.6064 -0.8325 0.5541 -0.8706
0.5541 -0.8706 0.4994 -0.9055
0.4994 -0.9055 0.4425 -0.9368
0.4425 -0.9368 0.3837 -0.9646
0.3837 -0.9646 0.3231 -0.9887
0.3231 -0.9887 0.2609 -1.0090
0.2609 -1.0090 0.1975 -1.0254
0.1975 -1.0254 0.1331 -1.0379
0.1331 -1.0379 0.0679 -1.0462
0.0679 -1.0462 0.0022 -1.0505
0.0022 -1.0505 -0.0638 -1.0506
-0.0638 -1.0506 -0.1299 -1.0466
-0.1299 -1.0466 -0.1956 -1.0384
-0.1956 -1.0384 -0.2609 -1.0262
-0.2609 -1.0262 -0.3253 -1.0098
-0.3253 -1.0098 -0.3888 -0.9893
-0.3888 -0.9893 -0.4509 -0.9649
-0.4509 -0.9649 -0.5116 -0.9366
-0.5116 -0.9366 -0.5704 -0.9044
-0.5704 -0.9044 -0.6272 -0.8686
-0.6272 -0.8686 -0.6818 -0.8292
-0.6818 -0.8292 -0.7339 -0.7863
-0.7339 -0.7863 -0.7833 -0.7402
-0.7833 -0.7402 -0.8298 -0.6910
-0.8298 -0.6910 -0.8732 -0.6389
-0.8732 -0.6389 -0.9134 -0.5840
-0.9134 -0.5840 -0.9501 -0.5266
-0.9501 -0.5266 -0.9832 -0.4669
-0.9832 -0.4669 -1.0125 -0.4051
-1.0125 -0.4051 -1.0380 -0.3415
-1.0380 -0.3415 -1.0594 -0.2763
-1.0594 -0.2763 -1.0768 -0.2097
-1.0768 -0.2097 -1.0900 -0.1421
-1.0900 -0.1421 -1.0989 -0.0736
-1.0989 -0.0736 -1.1035 -0.0045
-1.1035 -0.0045 -1.1038 0.0648
-1.1038 0.0648 -1.0997 0.1341
-1.0997 0.1341 -1.0913 0.2032
-1.0913 0.2032 -1.0785 0.2718
-1.0785 0.2718 -1.0614 0.3396
-1.0614 0.3396 -1.0401 0.4063
-1.0401 0.4063 -1.0146 0.4716
-1.0146 0.4716 -0.9849 0.5354
-0.9849 0.5354 -0.9513 0.5973
-0.9513 0.5973 -0.9138 0.6570
-0.9138 0.6570 -0.8725 0.7144
-0.8725 0.7144 -0.8276 0.7693
-0.8276 0.7693 -0.7793 0.8213
-0.7793 0.8213 -0.7277 0.8702
-0.7277 0.8702 -0.6730 0.9160
-0.6730 0.9160 -0.6154 0.9582
-0.6154 0.9582 -0.5552 0.9969
-0.5552 0.9969 -0.4926 1.0318
-0.4926 1.0318 -0.4278 1.0627
-0.4278 1.0627 -0.3610 1.0896
-0.3610 1.0896 -0.2925 1.1123
-0.2925 1.1123 -0.2226 1.1307
-0.2226 1.1307 -0.1516 1.1447
-0.1516 1.1447 -0.0797 1.1542
-0.0797 1.1542 -0.0072 1.1592
-0.0072 1.1592 0.0657 1.1597
0.0657 1.1597 0.1385 1.1555
0.1385 1.1555 0.2111 1.1468
0.2111 1.1468 0.2832 1.1336
0.2832 1.1336 0.3544 1.1158
0.3544 1.1158 0.4245 1.0935
0.4245 1.0935 0.4932 1.0668
0.4932 1.0668 0.5603 1.0358
0.5603 1.0358 0.6254 1.0006
0.6254 1.0006 0.6882 0.9613
0.6882 0.9613 0.7486 0.9181
0.7486 0.9181 0.8063 0.8711
0.8063 0.8711 0.8610 0.8204
0.8610 0.8204 0.9126 0.7663
0.9126 0.7663 0.9607 0.7090
0.9607 0.7090 1.0053 0.6486
1.0053 0.6486 1.0460 0.5854
1.0460 0.5854 1.0828 0.5197
1.0828 0.5197 1.1155 0.4517
1.1155 0.4517 1.1438 0.3816
1.1438 0.3816 1.1678 0.3097
1.1678 0.3097 1.1873 0.2363
1.1873 0.2363 1.2021 0.1617
1.2021 0.1617 1.2123 0.0862
1.2123 0.0862 1.2177 0.0100
|
The figure above shows the same results
The figure below shows similar results
|
This module covered a lot of territory. We began by looking at some data from the Mackenzie River District and saw apparently cyclical behavior for the lynx population. Then we looked at one model of predator-prey interaction. This model appeared to predict the same kind of cyclical behavior we saw in the data but there were some difficulties. Our numerical simulations of this model were inconclusive. This lead us to examine another situation -- spring-and-mass systems -- in which we also see cyclical behavior (if there is no friction) and our knowledge of physics helps us understand what is happening. Finally, we looked at Euler's Method, the numerical method used in our Java applets in this module, to try to understand what is going on mathematically.
This module has lots of mathematical messages -- in particular, the difficulty in using just numerical and graphical methods to analyze cyclical behavior. But, there is another message as well. The difficulty we observed using numerical estimates is a tip-off that these simple models may not be very good for describing biological models in which we expect lots of "noise." More realistic and more complex biological models are likely to involve stable cycles rather than the unstable cycles we saw in this module.
The CAS laboratories for this module set-up Euler's Method. You should use your CAS window to repeat the investigations we did above for spring-and-mass systems and for our predator-prey model.
