Nonlinear Systems of Differential Equations—Predator-Prey ...



Nonlinear Systems of Differential Equations—Consumer-Resource Models

Nonlinear, autonomous systems of ordinary differential equations are of the form

[pic]

where each of the functions fi on the right-hand side are real-valued functions in n variables. Most of the time, we will restrict the analysis to systems of two variables. We will focus on equilibria and stability.

Equilibria and Stability

Consider the system of two autonomous differential equations

(1) [pic]

The first step is to find the equations of the zero isoclines, which are defined as the set of points that satisfy

[pic]

Each equation results in a curve in the x-y space. Point equilibria occur where the two isoclines intersect (Figure 1). A point equilibrium [pic] of (1) therefore simultaneously satisfies the two equations

[pic]

We will call point equilibria simply “equilibria.”

[pic]

Figure 1: Zero isoclines corresponding to the two differential equations. Equilibria occur where the isoclines intersect.

The analytical approach to stability relies on analyzing the effects of small perturbations. We say that the equilibrium [pic] is locally stable if the system returns to [pic] after a small perturbation, and unstable otherwise. Mathematically, this can be analyzed through linearizing the right-hand side of each of the two differential equations in (1) about the equilibrium.

The system (1) can be written in matrix form

[pic]

where the right-hand side is a vector-valued function that maps a point in R2 (the two-dimensional real plane) into a point in R2. To linearize a vector-valued function, we need to linearize each component separately. Linearizing a function of two variables about a specific point means to find the tangent plane at this point (this, of course, is not always possible). The equation of a tangent plane of f(x,y) about [pic] is given by

[pic]

We thus find for the linearization of the vector-valued function [pic]

[pic]

Now, considering [pic] and [pic] as the perturbations, then with [pic] and [pic], we find

(2) [pic]

The matrix

[pic]

is called the Jacobi matrix.

The system (2) is a linear system of two equations, and we can use the results from linear systems of two differential equations to determine the stability of the equilibria. Namely,

• The equilibrium is a node if both eigenvalues of the Jacobian evaluated at the equilibrium are real, distinct, nonzero, and are of the same sign. The node is locally stable if the eigenvalues are negative, and unstable if the eigenvalues are positive.

• The equilibrium is a saddle if both eigenvalues of the Jacobian evaluated at the equilibrium are real and nonzero but have opposite signs. A saddle is unstable.

• The equilibrium is a spiral if both eigenvalues of the Jacobian evaluated at the equilibrium are complex conjugates with nonzero real parts. The spiral is locally stable if the real parts of the eigenvalues are negative, and unstable if the real parts of the eigenvalues are positive.

Task 1:

Show that (3,2) is an equilibrium of

[pic]

and determine its stability.

Task 2:

Consider the following set of two differential equations:

[pic]

(a) Graph the zero isoclines into the vector field (see handout in class).

The vector field was created in Matlab using

>> [x,y] = meshgrid(-2:.1:4,-2:0.1:4);

>> M=x.*(10-2*x-y);

>> N=y.*(10-x-2*y);

>> quiver(x,y,M,N,2)

>> grid on

Note that the multiplication is “.*” since we are multiplying componentwise.

(b) Show that [pic] is an equilibrium, and determine its stability using the analytical approach.

Task 3:

Find all equilibria of

[pic]

and use the analytical approach to determine the stability of each equilibrium. Graph the zero isoclines to confirm that you found all equilibria.

Euler Method for Systems of Differential Equations

To numerically solve a system of differential equations, we can use the Euler method. This is the simplest approximation scheme. It has disadvantages, like slow convergence and instabilities. There are far superior numerical schemes (like the Runge-Kutta method), but the Euler scheme suffices for our purposes here.

The Euler scheme is directly based on the definition of the derivative. Namely,

[pic]

To approximate the system of differential equations

[pic]

we approximate the left-hand side using the definition of the derivative

[pic]

and then rewrite this as

[pic]

which can be coded up in a spreadsheet. The time step h needs to be small enough so that the approximation is sufficiently good but not too small so that rounding errors accumulate too quickly.

Consumer-Resource Models

Model 1: Lotka-Volterra Predator-Prey Model

There are a large number of models that deal with consumer and resource interactions. The oldest such model is the Lotka-Volterra model, which describes the interaction between a predator and its prey. If x denotes the abundance of the prey and y the abundance of the predator, then the model is given by the following set of differential equations:

(3) [pic]

where the parameters b, c, e, and d on the right-hand side are positive constants. The parameter b is the prey birth rate, c is the consumption rate, e is the conversion fraction of prey into new predators, and d is the predator death rate. The equilibria can be obtained algebraically by setting the right-hand sides equal to 0 and solving for x and y. We find the trivial equilibrium (0,0) and a nontrivial equilibrium [pic]. The corresponding Jacobi matrix is given by

[pic]

If we evaluate this at the trivial equilibrium (0,0), we find

[pic]

Therefore, the two eigenvalues are [pic] and [pic]. Both eigenvalues are real. Since one eigenvalue is positive and the other is negative, we conclude that the trivial equilibrium is a saddle and thus unstable.

If we evaluate the Jacobi matrix at the nontrivial equilibrium[pic], we find

[pic]

To compute the eigenvalues, we solve

[pic]

which yields the two complex conjugate eigenvalues [pic] and [pic]. The real parts of these two eigenvalues are both equal to 0. The linear stability analysis is thus inconclusive.

It turns out that the equilibrium is neutral and this system of equations exhibits neutral oscillations. This will be explored numerically in the computer lab exercises.

Model 2: Density-dependent Growth of the Prey

It is possible to stabilize the predator-prey model by including density-dependent growth of the prey in the form of logistic growth. This takes the form

(4) [pic]

where all parameters are positive. In the absence of the predator, the prey dynamics reduce to logistic growth in Equation (4). Namely, If we set [pic], then [pic].

The system of equations (4) has the nontrivial equilibrium (0,0), which is always unstable. In addition, it has the prey only equilibrium (K,0), which is locally stable provided [pic]. If [pic], an additional nontrivial equilibrium in the first quadrant appears, which is locally stable. If [pic], the prey only equilibrium is unstable.

Model 3: Saturating Functional Response

The functional response of a predator is defined as the per predator predation rate. In this case, the functional response is simply cx, a linear function of prey density. In nature, one often observes functional response curves that saturate with increasing prey density. This happens, for instance, if predators need to spend time handling prey.

A common form for a saturating functional response is

(5) [pic]

This model is known as Rosenzweig’s model. It has the trivial equilibrium (0,0) and the prey only equilibrium (K,0). If [pic], there is an additional equilibrium in the first quadrant, which may be locally stable or unstable. The interesting property of this model is that when increasing the carrying capacity K, the equilibrium becomes unstable. This is known as Rosenzweig’s paradox of enrichment. (You can read more about this in the attached paper.)

Task 4:

Graph the functional response curve [pic] for c=0.81 and a=0.25. The parameter c is called the saturation level since [pic], and a is the half-saturation constant since [pic]. Label a and c on your graph.

Computer Lab (In Class)

In three separate worksheets (Tabs: “LV,” “LV with DD,” and “Rosenzweig”), we will code up the three different predator-prey models using the Euler method. Use the following parameters:

[pic]

Use initial conditions [pic] and [pic]. As the time step, choose h=0.007.

Your worksheet for the first model (Tab “LV”) looks as follows:

|  |A |B |C |D |E |F |G |H |

|2 |1.2 |0.25 |0.81 |0.8 |0.25 |0.007 |0.5 |0.25 |

|3 | | | | | | | | |

|4 |Time |Producer |Consumer | | | | | |

|5 |0 |0.5 |0.25 | | | | | |

|6 |0.007 |0.503491 |0.25013 | | | | | |

|7 |0.014 |0.507007 |0.250263 | | | | | |

|8 |0.021 |0.510546 |0.250401 | | | | | |

|9 |0.028 |0.51411 |0.250542 | | | | | |

Row 1 contains the parameter names and Row 2 the corresponding values. In A5, enter 0. Cell A6 is “=A5+$F$2”. Drag this cell down so that the final time is about 32. In B5 and C5, copy “=$G$2” and “=$H$2”. In B6, write “=B5+$F$2*($A$2*B5-$C$2*B5*C5)” and in C6, write “=C5+$F$2*($D$2*$C$2*B5*C5-$E$2*C5)”. Drag both cells down to about row 5000.

Investigate the behavior of the three models as the carrying capacity K changes. Simulate the three models for K=0.1, 0.25, 0.75, 1.0, and 2.0, and describe what you see.

Literature Cited

Rosenzweig, M.L. (1971) Paradox of enrichment: destabilization of exploitation

ecosystems in ecological time. Science 171: 385-387.

-----------------------

x

y

[pic]

[pic]

Equilibrium

................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download