Worksheet 4: Ordinary Differential Equations
Ordinary Differential Equations
Ordinary differential equations are equations of the form
(1)
On the left hand side is the derivative of the dependent variable x with respect to the independent variable t. On the right hand side, there is a function that may depend on both x and t. The independent variable t often represents time. In contrast to discrete time equations of the form , where time t is discrete (), the independent variable t in Equation (1) is a continuous variable, that is, it takes on real values, for instance, . In addition, we prescribe the initial value at time 0, namely x0. (The initial condition could be stated for some other time but time 0 is quite commonly used.) A differential equation is called ordinary if it involves only one independent variable. We begin with autonomous differential equations. These are ordinary differential equations where the right hand side of Equation (1)does not explicitly depend on the independent variable t. They arethus of the form
(2)
where f(x) is some function. Examples of this type of equation are, for instance, exponential growth or logistic growth,. In both examples, the independent variable t represents time, and the dependent variable N represents population size. Both differential equations model the growth of a population. We discuss both of them in the following.
Exponential Growth
The differential equation for exponential growth with parameter r, where, and initial condition is given by
(3)
To understand where this equation comes from, we look at the definition of a derivative
We see from this that the following approximation for small holds:
from which we conclude that
or
This now is in the form of discrete time exponential growth where each time step is of length .
We can also relate the discrete time and continuous time growth by starting with discrete time growth, which is given by
At time 0, we assume that the population size is equal to . If we measured population size every h units of time then every unit of time where we think of h as being a number between 0 and 1, we would write for the dynamics of this population
Note that the growth parameter has a subscript that indicates the length of the time step. We can compute how much the population size changes from time step to time step, namely
If we divide both sides by h and let h tend to 0, we arrive at
(4)
The left-hand side converges to the derivative of N with respect to t (provided the limit exists), and the right-hand side converges to if the limit exists and we denote
(5)
Equation (4) together with Equation (5) then yields the ordinary differential equation describing exponential growth as given by Equation (3).
To see why the differential equation in (3) models exponential growth, we solve it. That is, we want to find a function that satisfies Equation (3), including the initial condition. To solve Equation (3), we use separation of variables followed by integration.
Carrying out the integration, we find
where C1 is the constant of integration. The following algebraic steps simplify the expression
We can now incorporate the initial condition
We thus find as the solution of (3)
We see that the solution is the exponential function. If , then N(t), the population size at time t, will be 0 as well. This is thus a constant solution of (1). Furthermore, if r<0, then, for any initial population size, N(t) will approach 0 as time tends to infinity. If r>0, then N(t) will go to infinity when N(0)>0.
We can use the symbolic toolbox in Matlab to solve differential equations. Enter into the Command Window
clear
syms xtr
x=dsolve('Dx=r*x','x(0)=x0','t')
If you want to evaluate the solution, you need to convert the string x into a function. This is done using the inline function
xsol=inline(char(x),'t','x0','r')
To find the solution when , type into the Command Window
t=2;x0=0.5;r=2;
xsol(t,x0,r)
To plot the function, you need to evaluate the solution for different values of t and then use the plot function. First, however, you need to tell Matlab that you want a vector function:
xsol=inline(vectorize(char(x)),'t','x0','r')
We then generate values for the time vector:
t=0:0.1:10;
We can evaluate the function using the same command as before
xsol(t,x0,r);
Plotting is done as follows
plot(t,xsol(t,x0,r));xlabel('t');ylabel('x(t)');title('Exponential Growth');
Logistic Growth
The differential equation for logistic growth with the two parameters r>0 and K>0
(6)
can also be solved explicitly using separation of variables followed by integration. It is more complicated, though: Integration involves a partial fraction decomposition. The solution of (6) is then given by
provided . If , then for all later times. We see that if , then for all later times. There are thus two constant solutions for this equation, namely 0 and K. Furthermore, N(t) will approach K as time tends to infinity for any positive initial population size.
Numerical Methods for Solving Differential Equations
Many differential equations cannot be solved exactly. Numerical methods have been developed to approximate solutions. Numerical analysis is a field in mathematics that is concerned with developing approximate numerical methods and assessing their accuracy, for instance for solving differential equations. We will discuss the most basic method, the Euler method, explain the basic ideas behind this numerical method, and how to implement it on a spreadsheet using the example of exponential and logistic growth.
Approximating the Solution of a Differential Equation
An initial value problem can be written as
(7)
We assume that f and are continuous in a rectangle , which contains the point . This then guarantees that (7) has a solution in some interval where .
If it is not possible to solve an initial value problem analytically, numerical methods are employed. The idea is to use a discrete set of points, , (for instance, equally spaced points , where h is a positive constant), and then to find values yn so that yn approximates the solution at the points xn, that is,
Euler’s Method
We consider the initial value problem
(8)
and assume it has a solution in some interval where . The simplest scheme for approximating the solution is Euler’s method, which can be implemented in a spreadsheet. Euler’s method is not the best method in the sense of accuracy but is clearly the simplest. It is directly based on the definition of the derivative. Recall that the definition of a derivative is
provided the limit exists. When h is sufficiently small, we can write
We say that the error in approximation is of order . This means that the error between the exact solution and the approximation tends to 0 like a constant times as h tends to 0.
This procedure leads to the approximation scheme
(9)
This method is a first-order method since the error in each step is of order and the total number of steps is of order , resulting in a global error of h.
Euler’s Method for Autonomous Differential Equations
An autonomous differential equation is of the form
and we can write the initial condition as before, namely. The approximation scheme simplifies to
(10)
Task 1:
Look at Equation (10) and the graph below: If and , find and in the following graph. Where in the graph is the error of approximation we made in the first iteration of the approximation scheme?
Figure 1: Euler’s method
Exponential Growth and the Euler Approximation
We illustrate the iteration scheme for exponential growth. We choose a specific example, namely we consider the initial value problem
(11)
We know what the solution of (11) is, namely . We numerically solve (11) using Euler’s method with h=0.1 in the time interval [0,0.5], and then check how well this method performs. We have . Then
and
We summarize this in the following table. If h=0.1, then
x / y / Exact / Difference0 / 5 / 5 / 0
0.1 / 6 / 6.107014 / 0.107014
0.2 / 7.2 / 7.459123 / 0.259123
0.3 / 8.64 / 9.110594 / 0.470594
0.4 / 10.368 / 11.1277 / 0.759705
0.5 / 12.4416 / 13.59141 / 1.149809
The third column contains the exact values, . The last error contains the absolute error after each step, computed as |y-Exact|. We see that when h=0.1, the numerical approximation is not very good after five steps. If we repeat the same approximation with a smaller value for h, say h=0.01, the following table results for the first five steps:
x / y / Exact / Difference0 / 5 / 5 / 0
0.01 / 5.1 / 5.101007 / 0.001007
0.02 / 5.202 / 5.204054 / 0.002054
0.03 / 5.30604 / 5.309183 / 0.003143
0.04 / 5.412161 / 5.416435 / 0.004275
0.05 / 5.520404 / 5.525855 / 0.005451
Doing five steps only gets us to x=0.05. We can do more steps until we reach x=0.5. We find that
x / y / Exact / Difference0.5 / 13.45794 / 13.59141 / 0.133469
Choosing a smaller value for h resulted in a better approximation at x=0.5 but also required more steps. One source of error in the approximation comes from the approximation itself. Another source comes from rounding errors when we implement the approximation on a computer. It is therefore not necessarily the case that smaller values of h always improve the approximation.
Logistic Growth
Logistic growth is used to describe the growth of a population whose growth rate declines with increasing population size. Its mathematical formulation is in terms of a differential equation
(12)
where is the population size at time t. The parameter r is the intrinsic rate of growth and the parameter K is called the carrying capacity. The population dynamics evolve in continuous time, that is, changes occur at any time not just at discrete time steps. This modeling scheme is appropriate when a population does not have discrete breeding times. An example is the growth of Paramecium aurelia in a test tube containing some culture medium with bacteria as food (Gause 1934).
Figure 2: Population size of Paramecium aurelia as a function of time. The data is from Gause (1934). A logistic growth curve is fitted to the data. (Source: Case, 2000, page 104)
With the initial condition N(0)=N0, we found earlier that the solution is
(13)
Task 2:
This task refers to logistic growth and uses Equations (12) and (13). Use a spreadsheet to plot the following graphs for N0=5, r=2 and K=20: (a) the growth rate dN/dt as a function of the population size N [use Equation (12)], (b) the per capita growth rate dN/Ndt as a function of the population size N [divide both sides of Equation (12) by N], and (c) the solution N(t) as a function of t [use Equation (13)].
Task 3:
This task refers to logistic growth and uses Equation (13).Use a spreadsheet to plot for r=2 and K=20 the solution N(t) in Equation (13) as a function of t for the following initial conditions: (a) N0=1, (b) N0=20, and (c) N0=25. Plot all three graphs in one coordinate system. Describe what you see.
Euler’s Method
The differential equation for logistic growth can be solved analytically. Now, however, we will set up Euler’s method to approximate the solution of the logistic equation
with , and then compare the approximate solution to the exact solution in the computer lab exercises. As explained above, Euler’s method is based on the definition of the derivative. If we choose a sufficiently small time step h, we can write
or, after rearranging the approximation,
In the case of the logistic growth equation, we have , , and Euler’s scheme is
(14)
This can be used iteratively and thus implemented into a spreadsheet. Let’s see how the first few steps work. Start with specifying the initial population size N0 and the two parameters r and K. Then choose a step size h sufficiently small. The approximation can be used to compute N(h)
(15)
This gives an approximation for N(h). Using N(h), we can compute N(h+h)=N(2h), then N(2h+h)=N(3h), and so on.
Theoretically, the accuracy of this method improves with decreasing h. However, as mentioned earlier, there are two important limitations: Decreasing h increases the computing time and the cumulative errors due to rounding.
Task 4:
When h=1, Equation (15) becomes
(16)
Compare this to the recursion of the discrete time logistic equation
(17)
We discussed (17) in Worksheet 3. How is r in (16) related to R in (17)?
Equilibria and Stability
The two examples, exponential growth and logistic growth, introduced two important concepts: constant solutions and long-term behavior. Constant solutions are called point equilibria and the long-term behavior tells us something about the stability of these point equilibria.
Point equilibria of a differential equation are defined as the solutions of . Point equilibria are the constant solutions of the differential equation.
Task 5:
Find all equilibria of (a) exponential growth and (b) logistic growth.
Suppose that is a point equilibrium of the differential equation . We say that is locally stable if the solution returns to after a small perturbation. Otherwise we call unstable. There is a graphical and an analytical way to determine stability.
Graphical Analysis
For the graphical analysis, we graph f(x) as a function of x. The point equilibria are the points of intersection of the graph with the horizontal axis. To determine stability, we check how behaves in a neighborhood of each point equilibrium (Figure 3).
Figure 3: Locally stability can be analyzed graphically. The arrows indicate whether x will increase or decrease.
Analytical Analysis
The analytical analysis is based on linearization about an equilibrium. The linearization of f(x) about the equilibrium is
A perturbation about the equilibrium can be written as . Now, since , we find
for z sufficiently small. Comparing this with exponential growth, it follows that is locally stable if , and unstable if . When , the linear stability analysis is inconclusive. Looking at Figure 3 again, we can confirm that this characterization agrees with the graphical analysis. We call an eigenvalue.
Homework (Hand in on ______)
Step 1
Set up a spreadsheet to code up Euler’s method for exponential growth
for step size h=0.1 and time interval [0,5]. This is Equation (11) and your result should agree with what we obtained in the worksheet earlier. Paste the table in your report.
Step 2
Set up a spreadsheet to code up Euler’s method for logistic growth for the parameters N(0)=5, r=2, and K=20, and step size h=0.1. Your spreadsheet should look like
A / B / C / D / E1 / Logistic Growth: Euler's Method
2
3 / Parameters
4 / r / K / N(0) / h
5 / 2 / 20 / 5 / 0.1
6
7
8 / Time / N(t)--Euler / N(t)--Solution / relative error / maximum error
9 / 0 / 5 / 5 / 0 / 0.015421746
10 / 0.1 / 5.75 / 5.786715119 / 0.006344725
11 / 0.2 / 6.569375 / 6.642399462 / 0.010993687
12 / 0.3 / 7.451683121 / 7.55733683 / 0.013980283
13 / 0.4 / 8.386743932 / 8.517935115 / 0.015401759
14 / 0.5 / 9.360717981 / 9.507337728 / 0.015421746
Enter the parameters as numbers in row 5. Enter “0” in cell A9. Enter “=$C$5” in cells B9 and C9. The relative error is computed as “=ABS((C9-B9)/C9)” and this goes into cell D9. Write “=A9+$D$5” in cell A10 and drag this cell down to row 109. The approximate solution based on Euler’s method goes into cell B10: “=B9+$D$5*$A$5*B9*(1-B9/$B$5)” and the exact solution goes into cell C10: “=$B$5/(1+($B$5-$C$5)*EXP(-$A$5*A10)/$C$5)”. Drag both cells down to row 109. Finally, to compute the maximum error, write “=MAX(D9:D109)” into cell E9.
Use Euler’s approximation to solve the logistic equation for this set of parameters numerically. Plot both the numerical and the analytical solution into one coordinate system for the time interval [0,10].
Step 3
By changing the step size by increments of 0.1 from 0.1 to 1.0 (i.e., 0.1, 0.2, 0.3, …, 1.0), you can investigate how the accuracy of the solution is affected by the choice of the step size. (Note that you showed in Task 4 that if h=1, the model becomes the discrete logistic growth equation.) In your lab report, provide the graphs for the parameters N(0)=5, r=2, and K=20, and h=0.3, 0.7, and 1.0 for the time interval [0,10]. Repeat the same but now set r=2.8 and extend the time interval to [0,30]. Report any interesting observations.
Step 4
(a) You have now become familiar with difference and differential equations. Compare the two modeling approaches and comment on similarities and differences addressing the following issues.
- Provide one biological example for each modeling approach that illustrates the difference between discrete time and continuous time modeling.
- How do you find equilibria in each of the two modeling approaches?
- How do you analyze stability of equilibria in each of the two modeling approaches?
- How do you implement each modeling approach in EXCEL, including whether the solution is exact or an approximation?
(b) When you look back over the models we have discussed in more detail, what aspects of the models surprised you (if any)?
(c) Reflecting back on difference and differential equations, what do you still find unclear?
Homework (Hand in ______)
Problem 1 (Stoichiometry):
(a) The nutrient content of a consumer (for instance, the percent nitrogen of the consumer’s biomass) can range from reflecting the nutrient content of its food to being constant. The former is referred to as absence of homeostasis, the latter as strict homeostasis. A model for homeostatic regulation is provided in Sterner and Elser (2002). The model relates consumer’s nutrient content (denoted by y) to its food’s nutrient content (denoted by x) as follows:
where is a constant. Solve the differential equation and relate θ to absence of homeostasis and strict homeostasis.
(b) When phosphorus content of Daphnia was plotted against phosphorus content of its algal food on a log-log plot, a straight line with slope 1/7.7 resulted (see Sterner and Elser, 2002; data from DeMott et al., 1998). Find a differential equation that relates phosphorus content of Daphnia to the phosphorus content of its algal food.
Problem 2 (Logistic Growth):
Assume that the size of a population evolves according to the continuous time, logistic growth equation with intrinsic rate of growth r=1.5. Assume that the carrying capacity K=100.
(a) Find the differential equation that describes the rate of growth of this population.
(b) Write a Matlab code to solve the differential equation symbolically.
(c) Find all equilibria, and discuss the stability of the equilibria using the graphical approach.
(d) Find the eigenvalues associated with the equilibria (i.e., find where denotes the equilibrium), and use the eigenvalues to determine stability of the equilibria. Compare your answers with your results in (b).
Problem 3 (A Simple Model of Predation):
Suppose N(t) denotes the size of a population at time t. The population evolves according to the logistic growth equation, but, in addition, predation reduces the size of the population so that the rate of change is given by
(18)
The first term on the right-hand side describes the logistic growth; the second term describes the effect of predation.
(a) Set