This assignment is due at the start of lecture on Monday, 9 Feb; however, you ought to be able to finish during the lab period. Your solution should include all codes and plots, as well as written responses where appropriate.

Since we'd like our ode solvers to be generic, we need a way to pass the name of a function in MATLAB as if it were an ordinary variable. This is done via function handles (if you're familiar with C or C++, these are known as pointers to functions). Recall that the Euler step we covered in class assumes a derivative can be written as dy/dx = f(x,y); however, in practice we often have derivatives that are functions of other parameters we'd like to change in order to study a dynamical system. A good example would be simple exponential decay: dN/dt = -N/τ. We'd like to change the time constant τ without re-writing the code. Here's how that can be done in MATLAB:

• assume you have written a function (you can download this function directly: expDec.m, just as it was demonstrated in class) that calculates a derivative for exponential decay. It is called as:

  dNdt = expDec(t,N,tau);

• the function handle would be assigned as follows:

  foo = @ (t,N) expDec(t,N,tau)

• this syntax is telling MATLAB that foo refers to a function whose variables are t and N, with an additional parameter tau that is already defined
• try this series of statements from the command line to see how this syntax works:

  tau = exp(1), t = pi, N = 1e9
  foo = @ (t,N) expDec(t,N,tau)
  dNdt_func = expDec(t,N,tau)
  dNdt_handle = foo(t,N)
  dNdt_func - dNdt_handle
  

• notice that there's no difference in calling the derivative function via its 3-input standard call expDec and the generic 2-argument call foo via the handle — this is very convenient for an Euler or Runge-Kutta method, which generically assumes 2-input calls

For this lab you may use your own codes or make modifications to the codes I've provided on the course code site. A more-than-complete example for exponential decay is given by run_expDec.m

1. Physicists study models of biological populations because they often exhibit rich dynamics. 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. This 1^st-order ode describes the population as a function of time.
1. Which of the constants represents a birth rate?
2. Which of the constants represents a death rate?
3. Explain your answers.
2. Write a MATLAB function to calculate the Verhulst derivative at a particular point in time. This function should be compatible with the generic ode solvers I presented in class. Hint: since we're setting this function up to be called by a generic ode solver, make sure your function has the appropriate inputs, even if they are not all used in the function! Recall, we're going to assume a generic derivative has the form dy/dx = f(x,y).
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) = N_oe^at. Briefly describe your results, including a physical description of what you observe.
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?
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.
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?