PHYS 220 — Lab #2 [Solution]

  1. One population model (owing to Verhulst) is given by dN/dt = aN - bN², where N is the number of individuals in a population; a and b are positive constants.
    1. Which of the constants represents a birth rate? a
    2. Which of the constants represents a death rate? b
    3. Explain your answers. Since aN is positive, it governs the increase in population over time; similarly, bN² governs the decrease.


  2. Write a MATLAB function to calculate the Verhulst derivative at a particular point in time. A working, commented code for this is posted at verhulst.m. Note the ordering of its call arguments: independent variable, dependent variable (stored as a vector), and other parameters.


  3. Write a MATLAB function or script that calculates a solution to the Verhulst population ode, above, using either your own your Euler step function or my generic one (odeStep_euler.m) and a function handle to your derivative function. Solve the equations for t ∈ [0,10], a = 0.5, b = 0 and N(0) = 1e3. Plot your results along with a comparison to the exact result N(t) = Noeat. Briefly describe your results, including a physical description of what you observe.
    For this solution, I'll make use of the paradigm used by odeStep_euler.m. My derivative function is called verhulst.m. A sample code to "drive" all of this (often called a "wrapper") is given by run_verhulst.m. Note that you can simply change parameters in this script to accomplish multiple parts of the lab & not have to waste time cutting & pasting lines to the command line. Below you'll see a plot of the results for this part of the lab. Note that the smaller the time step (we use Δt, Δt/2 and 2Δt) the better we approximate the exact solution. Note also that we do not connect the dots since we're only calculating a solution at discrete points. The exact curve is OK, since we have an exact solution for all t. in the plot we witness exponential growth (exactly) and see that the errors begin to diverge, no matter what our time step. This is global divergence.



  4. Change a to -a and re-run your simulation. Comment on the change in the errors you notice from the solution with a>0. Does the a parameter still represent the birth rate in this scenario? For these parameters, we've just changed the problem to exponential decay. Note that the solutions are now globally convergent.



  5. Modify your function or script for t ∈ [0,1] with the parameters a = 10, b = 2 and N(0) = 10. Describe your results. Since with the non-zero death rate you no longer have an expression for the exact result, briefly discuss a strategy for estimating the errors your code is producing. Note here I had to cut my time step by an order of magnitude – more on this topic when we get to partial differential equations. here we see something approximating exponential decay. Although you weren't asked to find it, I've plotted the exact result for comparison. If you didn't know what the exact result was, you could estimate your errors by examining the differences between solutions using different time steps – this is the basis for looking at adaptive step sizes. We'll cover this in greater detail as the course develops.



  6. Finally, run your code one last time with the parameters a = 1e6, b = 5e-4 and N(0) = 2e9. Describe your results. Is your code working correctly? No need for a plot here – you should have a constant value N(t) = N(0), independent of time step or time. The reason is that aN(0) = bN(0)², so the differential equation encoded in verhulst.m gives you dN/dt = 0, a constant result.