APPM 2360 – Lab 1
Fish Population Modeling

Lab Goals and Instructions:

In this lab, we will model the population of various species of fish in a suburban lake.  This lab is designed to demonstrate the usefulness of basic numerical and qualitative techniques.  There are several questions listed at the end of the lab that must be addressed in the report.  Labs may be done in groups 3 or less.  One report must be turned in for each group.  Labs must include a cover page with each student's name and student number as well as their recitation section, Professor and TA. The lab is to be turned in during lecture on Friday, February 22, 2008.  Late labs will not be accepted under any circumstances.

Model

Populations can be modeled well by differential equations. Suppose we live in a fairly suburban area with a lake stocked with fish which people may try to catch. Since the lake is in a suburban area, we may assume, for simplicity, that people are the only predators of the fish. We want to estimate the population of fish in that lake because we want to make sure that the people do not catch all of the fish.

To do this, the logistic population model with harvesting is introduced:

 

If we rearrange the variables a bit and group some of the parameters, we can give the same model in the following form:

       

where y = y(t)  is the fish population as a function of time and a=r, c=r/L, H are constant parameters. This differential equation is a good start in modeling the fish population because it takes into account the fact that the lake provides only finite resources. However, it assumes a constant rate, H, of harvesting.  It makes more sense to assume that the amount of fish people catch is dependent upon the number of fish in the lake (the more fish, the more readily they are caught). Hence, a better model (not constant) for the harvesting  would be:

 

  


where H(y) represents how harvesting is related to y(t). Time t is in units of days and the population y is in units of hundreds of fish. The parameters a and c can be determined by biologists for a given species of fish. Two new constant parameters, p and q, are introduced in the predation term H(y). These represent the people's effectiveness at catching the fish.

This model assumes that a fixed number of fishing licenses have already been sold and that they do not restrict the number of fish a license holder may catch in a day. In this sense, we cannot control the people, but we can still choose the type of fish we decide to put in the lake.

Since the number of people fishing at any given time and the number of fish caught per person are not predictable, H(y) represents the harvesting effect on average. We are interested in whether the fish population is sustainable. This means that we are more interested in the long-term behavior of our population.

Since different species of fish will have different values of the parameters a and c and we have yet to decide which species we will use, we want to explore the long-term behavior of solutions of (1) for a range of the parameters a and c. By this analysis we will be able to predict what will happen to the populations of various different fish species over time.  Will they grow out of control?  Will they die out completely?  Will they reach some steady-state population?

Show that when the harvest rate H(y) is constant, a maximum sustainable harvest Hmax is a2/4c, which occurs when the population is half the carrying capacity.  (Hint:  Set  y’ = 0 and use the quadratic formula to find Hmax).  Solve the differential equation using separation of variables, assuming that H remains constant.  Now consider the case where H is a function of y as defined in equation (1).  Show that the limiting value of H(y) as the population y approaches infinity is a constant.  How does making this harvesting value non-constant affect the model compared to previously when we assumed it did not depend on y? 

It is difficult (but possible) to determine the exact solution of (1).  However this solution is a bit complicated and quite difficult to interpret, so instead of using the exact solution we will investigate (1) graphically.

Suppose we are considering stocking the lake with different types of fish; for all the types the parameter a = 0.65, but c ranges in value from 0.04 to 0.12 throughout the various species. We will fix p = 1.2 and q = 1.

Let c = 0.04 and define

 

so we can rewrite (1) as y'(t) = f(y). Use any graphing tool available to you to plot f as a function of y. Use this plot to determine (graphically) the equilibrium solutions of (1). [Hint: an equilibrium occurs when all time-derivatives are equal to zero simultaneously.]  Make sure you use a suitable scale when plotting so that you don't miss any equilibrium values.

Now change c incrementally by 0.04 until c = 0.12, finding the equilibria (by plotting, as above) at each stage.

By trial and error, find the approximate values of c where the number of equilibria changes (be sure that you are within 0.005 of the exact value). Be careful: within this range there are two c values where the number of equilibria change. The point where the number of equilibrium points changes is called the bifurcation value of c.

Now plot vector fields of (1) for the same increments of c as above. Again, you may use any vector field program available to you (see end of this lab instruction).

Use any ODE solver to plot solution curves to the differential equation for the same increments of c as above.  For each value of c, choose three representative initial values of y and overlay the solution curves corresponding to each of these on a single graph.  Make sure you use an appropriate time interval so that you can observe the long-term behavior of the solutions.  (One recommended tool for this is the ode45 command in MATLAB.  To learn more about this function, type “help ode45” in the MATLAB command window or search for it on the web.  You may also use similar functions in Mathematica or MVT).

Describe the long-term behavior of solutions for various initial conditions, paying close attention to the differences in behavior for values of c on either side of the bifurcation value. Note that you may need to take a (very) large time interval to see the true behavior of solutions.  Compare what you see here with how solutions behave with constant harvesting Hmax that you found previously. 

Along with reporting on the analyses above, your write-up should answer the following questions and address the following issues:

  1. Classify (1) and interpret (in words) each term of f(y).  By classify, it is meant that you should explain as much as you can about the form of the ODE.  What is its order?  Is it linear or nonlinear?  What makes it nonlinear?  Is it autonomous or non-autonomous and what makes it so?  Constant or variable coefficients?  Etc.

  2. When interpreting H(y), consider these questions: what happens to H(y) in the limiting cases as y gets very large or very small (we assume, logically, that we cannot have negative fish population)? What does this represent physically? Does this make sense? [Hint: it may help to remember the units of y' (and, therefore, H(y)).]

  3. What must be the units of the parameters a and c? (Remember what the units of y and t are.)

  4.    a.  

    5. Use separation of variables to solve (1) analytically when H(y) = 0 (i.e. p = 0).
       b.  Set up the integrals that one gets when separation of variables is done on the model that includes harvesting. Suggest an analytical technique that could be used to evaluate the integrals.

  5.    a.  Define in words what an equilibrium solution is in general.  Explain what this means in the case of the model at hand. 
       b.      Give an equivalent mathematical definition and explain how you used this to find the equilibria in the function plots.

  6. Interpret the function plots and vector fields:
       a.  Be sure to label and title all your graphs.
       b.  What is happening to the population when f(y) > 0? When f(y) < 0? When f(y) = 0?
       c.   Identify these behaviors on the vector fields. Make sure that your function plots and vector fields are consistent. Mark the equilibrium values on both. Also, make sure to scale your arrows on the vector fields appropriately, so they do not overwhelm the plot. 
       d.  What are the possible long-term behaviors? Explain how these behaviors depend on the initial conditions (population at time t = 0) and the value of the parameter c?
       e.  Your lab report should include a representative selection of the graphs to illustrate your conclusions.

  7. Suppose carp have c = 0.04, catfish have c = 0.08, and (delicious) bass have c = 0.12. Based on your results, which would be the best fish with which to stock the lake at the beginning of the season? (There may be various issues to consider -- be sure to explain the motivation behind your choice.)

  8. Briefly discuss the following issues:


    9.    a.  Are there any weaknesses in the model used? Explain.
       b.  What are some additional effects that the fish population model should take into account?
       c.  Can the model be improved and if so, how?
       d.  For the purposes of this study, would the researcher benefit greatly from obtaining a closed-form solution to (1)?  Why or why not?

Interesting link

Mathematical Visualization Toolkit. This is a link where you can find all the necessary mathematical tools (function plotter, ODE-solver, etc.) needed for this lab.  You are also welcome (and encouraged) to use MATLAB or Mathematica to complete your graphical analysis for this lab. 

This lab was originally written by David Trubatch and then modified by several billion people, the most important being Brock Mosovsky.  It is now the collaborative work of many great mathematical minds, and as such, you should feel privileged to perform these insightful analyses under their direction.