Qualitative analysis of differential equations
Alexander V Panfilov
TL;DR
This paper introduces fundamental qualitative analysis methods for differential equations, aimed at readers with minimal mathematical background, to facilitate understanding in mathematical biology.
Contribution
It provides an accessible introduction to qualitative analysis techniques for differential equations tailored for non-experts in mathematics.
Findings
Simplified methods for analyzing differential equations
Applicable to biological modeling contexts
Enhances understanding for non-mathematicians
Abstract
Here I introduce basic methods of qualitative analysis of differential equations used in mathematical biology for people with minimal mathematical background.
Click any figure to enlarge with its caption.
Figure 1
Figure 1
Figure 2
Figure 2
Figure 3
Figure 4
Figure 4
Figure 5
Figure 5
Figure 6
Figure 11
Figure 12
Figure 13
Figure 14
Figure 15
Figure 16
Figure 17
Figure 18
Figure 19
Figure 20
Figure 21
Figure 22
Figure 23
Figure 24
Figure 25
Figure 26
Figure 27
Figure 28
Figure 29
Figure 30
Figure 31
Figure 32
Figure 33
Figure 34
Figure 35
Figure 36
Figure 37
Figure 38
Figure 39
Figure 40Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsEducational Technology and Optimization · advanced mathematical theories
QUALITATIVE ANALYSIS OF
**DIFFERENTIAL EQUATIONS **
Alexander Panfilov
Contents
-
3.1 Differential equations of one variable and their solutions
-
3.2 Qualitative methods of analysis of differential equations of one variable
-
4.4 Phase portraits for complex eigen values: spiral, center
-
5.3 Determinant-trace method for finding the type of equilibrium
-
6 Graphical methods to study systems of differential equations
-
11.0.1 Solution of the initial value problem for a linear system
Chapter 1 Preliminaries
1.1 Basic algebra
1.1.1 Algebraic expressions
Algebraic expressions are formed from numbers, letters and arithmetic operations. The letters may represent unknown variables, which should be found from solutions of equations, or parameters (unknown numbers) on which the solutions depend.
Below, we review examples several basic operations which help us to work with algebraic express of ions.
One of the most basic algebraic operations is opening of parentheses, or simplification of expressions. For that we use the following rule:
[TABLE]
note, that here means , etc., as in algebra the multiplication is often omitted.
Example (open parenthesis):
Sometimes we use parentheses to factor expressions:
Example (factor expression):
In many cases we also need to work with fractions:
Example (the same denominator):
Example (different denominators):
Example (fractions simplifications):
To divide a fraction by another fraction we just need to multiply it by its inverse : ,
Example (division):
Also note that:
1.1.2 Limits
We call a limit of the function when approaches , if the value of get closer and closer to when takes values closer and closer to . We write it formally as:
[TABLE]
In many cases, finding the limit is trivial: we just need to substitute the value of into our function:
[TABLE]
Example:
Functions which have such property are called continuous and most of the functions used in biology are continuous. However, there are several important exception.
The first case, which will be the most important for us, is finding of limit of the function when . Finding such limits is important as it gives an asymptotic behaviour of our system when the size of a population becomes very large. Unfortunately, there is no such number ’’ which we can substitute into our function to find a limit using formula (1.2). For functions without parameters, we can guess the limit by substituting large values to (1.2), e.g. , but what to do for functions with parameters?
Let us discuss this problem for a special class of functions, which are the most relevant to our course, the so called rational functions , where and are polynomials. In that case we can always find the limit using the following property of the power function:
[TABLE]
where is an arbitrary constant and .
To prove it note that if approaches (becomes larger and larger), the power function with also becomes larger and larger and therefore will be closer and closer to zero, thus in accordance with the definition
To find the limit using this rule we need to do the following: (1) find the highest power of in our expression , (2) divide each term in our function by in that power, and (3) find the limit of each term using property (1.3). Let us consider three typical examples:
Example (find the limit): .
The highest power is , division gives: . The limits of the individual terms are:
Example (find the limit): , .
Similar steps give us: . This expression does not have sense as we cannot divide by zero and we do not have a finite limit for this function.
Example (find the limit): , .
.
Another non-trivial situation occurs when the denominator of our function becomes zero for some value of , for example for . In this case the formula (1.2) for limit cannot be used and other more careful analysis is necessary. If using of calculator we substitute some numbers into our function around point we will find the following: if becomes closer and closer to from the left, e.g. the function value becomes larger and larger, while if becomes closer and closer to from the right, e.g. the function value is negative and its absolute value also becomes larger and larger. We can formally write it as , while . However, in a strict sense, as there is no real number for which approaches for close to (from either side) thus the limit here does not exist.
We will use limits for drawing our functions and will see that limits at infinity give us horizontal asymptotes of our graphs, while blow up of functions for some (as in the last example) give us vertical asymptotes.
1.1.3 Equations
An equation is a mathematical relationship involving unknown variables. These unknowns are usually expressed by letters ’x’, ’y’, however in biology we use many other letters (e.g. ’N’, ’P’. ’T’, ’V’, etc.), which maybe somewhat related to the name of the species they describe. Solving equations means finding unknown(s) such that after substitution in the equation the left and right hand sides will be equal to each other. For example: equation has a solution , as .
The usual way to solve equations which have unknown variables in the first power only (linear equations), is to isolate the unknowns:
[TABLE]
We can achieve that by using the following rules of equation algebra: (1) we can multiply, or divide both sides of the equation by the same number, and (2) we can move numbers/expressions from one to the other side of the equation, by changing their sign. The proof of these rules is trivial. Indeed if two expressions are the same , then if we multiply (or divide) both of them by the same number , they still will be the same . Similarly if , we can add to the both sides of the equation, which will not change the equality, but we get: , or . Thus we see, that we were able to move from the right hand side to the left hand side of our equation, but it changed its sign as a result.
Example(solve equation): .
Solution:
For a quadratic equation we use the ’abc’ formula, which gives us the solutions as .
Example (solve equation): .
Solution: from the ’abc’ formula .
Equation may also contain parameters. A parameter is an unknown number (constant) that may have any value. It is different from the unknown variable, as the parameter is just a constant on which our solution depends.
Example (solve equation): , where in unknown variable and are parameters.
Solution: By multiplying both sides by we get , or , or , thus . We see that the solution depends on 3 parameters and if someone provides us with their values we will be able to find the solution by substituting the parameter values into the final formula.
1.1.4 Systems of equations
To solve a system of two linear equations we express one variable via the other and substitute it into the other equation.
Example (solve the system of equations): \left\{\begin{array}[]{l}2x+y=5\\ x+y=3\end{array}\right.
Solution: From the second equation we find , so we substitute into the first equation: , now substitute this value to and find , thus the solution is .
Unfortunately, there are no general rules to solve a system of nonlinear equations. The usual practical way is to start with a more simple equation, try to obtain from it as much information as possible and then substitute it to the other equation. It is also very helpful to factor expressions in order to simplify them.
Example (solve the system of equations): \left\{\begin{array}[]{l}2n-2n^{2}-2np=0\\ np-2p=0\end{array}\right.
Solution: From the second equation by factoring we find . The product is zero only if one of the multipliers is zero, thus we have two possibilities or . If we substitute into the first equation we find , or , thus for we have two solutions , or ; now substitute into the first equation: , , , thus for we found . Overall, we found the following three solutions of the given system .
Systems may also contain parameters.
Example (solve the system): \left\{\begin{array}[]{l}an-an^{2}-bnp=0\\ np-kp=0\end{array}\right., where are variables and , are the parameters.
Solution: We proceed similarly as in the previous case. From the second equation: , thus we have two cases or . After substituting into the first equation we get: , , thus , or ; after substituting into the first equation we get: , , , thus . Therefore, we found three solutions: . It is easy to see that if we substitute the parameter values to these formulas we obtain the solution of the previous problem. Note also, that for systems with parameters we need to be careful as not all operations are allowed for arbitrary parameter values. In our example in order to obtain the solution we had to make several divisions by parameters , and . However we can always do that as the parameters are positive numbers () and thus they cannot be equal to zero.
Finally note, that we can solve systems of three and more equations similarly, by subsequent substitutions from one equation to another, etc..
1.2 Functions of one variable
In science the relationships between quantities are normally expressed using functions. The simplest type of functions are functions of one variable. The function of one variable is a rule that allows us to find the value of a variable (number) from a single variable (number) . We denote it as . Below are examples of the most important functions:
- power functions , for example
[TABLE] 2. polynomials, , for example:
[TABLE] 3. rational functions :
[TABLE] 4. trigonometric functions , , :
[TABLE] 5. exponential and logarithmic function
[TABLE] 6. etc.
The derivative of a function at point is given by the following limit:
[TABLE]
The derivative shows the rate of change of a function at a point and has many important applications:
- •
If x(t) is the distance traveled by a car as a function of time , then gives the velocity of the car.
- •
If is the size of a population as the function of time, then gives the rate of growth of the population.
Geometrically, the derivative gives the slope of the tangent line to the graph of the function at the point (fig.1.1).
A graph of a line tangent to the function at point (fig.1.1) is given by the following equation:
[TABLE]
Equation (1.10) is also known as a linear approximation of function at point :
Let us check formula (1.10) by approximating the function at . We find: , hence . At this approximate formula gives . The exact value is . So the error is just 0.6%. However, if while the exact value is . So we see, that the approximate formula works good if is close to only.
Functions with parameters. Functions may depend not only on variable(s) but also on parameters. We have already seen the following example of the function that depends on three parameters :
[TABLE]
Equation (1.11) describes a general quadratic polynomial. If we choose, for example we will get the function given by equation (1.5). Studying functions with parameters allows us to obtain results for whole classes of functions. We will frequently use functions with parameters in our course. This is because biological models usually depend on many (up to hundreds) parameters and in many situations the exact values of these parameters are unknown. One of practical difficulties in working with parameters is that use of calculators is very limited, because calculators cannot do calculations with unknown quantities. The most valuable methods to study functions with parameters are direct algebraic computations and analysis of the obtained formulas. In this course we will widely use the graphical methods of representation of function. Let us start with review of the basic function graphs.
1.3 Graphs of functions of one variable
Example of graphs. We usually represent functions using graphs. To do that we plot the value of the variable along the -axis and the value of the function along the -axis. Let us start first by listing typical graph shapes that are important in this course.
Tips on graphs
Let us list important rules that may help to plot graphs of function with parameters.
The graph of the function can be obtained by a vertical shift by units of the graph of .
Example: In function , parameter just shifts the graph of by units above.
Important points of the graph are points at which the graph crosses the -axis (-intercept), given by , and points where the graph crosses the -axis (zeros of the function), given by . Note, that some graphs do not cross the or the axis and thus do not have -intercepts or zeros. For example graph of function (Fig.1.10a) does not have finite zeros or -intercepts.
Example: For function , the -intercept is . Zeros can be found from , which gives , or , or , thus zero is given by the formula , which is always valid as .
Another important graph feature are asymptotes. To find a horizontal asymptote we need to compute the . For functions without parameters, you can try to compute this limit using calculator by filling in a large numbers 10000, 20000, etc and looking if the function approaches some constant value. For functions with parameter, you can try to fill in some ’reasonable’ parameter value and try to find similarly if the asymptote exists, however the best way here is to find the limit using our plan from section 1.1.2. A vertical asymptote is usually a point where a denominator of a fraction is zero. Not all graphs have asymptotes, for example graph of function does not have any vertical or horizontal asymptotes. However, even if the asymptotes are absent it is still useful to understand behavior of the functions at large and show it in the graph.
Example: For function we can find as , thus this graph has a horizontal asymptote . The vertical asymptote here is at point where , or line .
Several features of the graph can be found from the derivative of the function: a function grows if its derivative , decreases if and has a local extremum (maximum or minimum) if . We do not necessarily need to compute these feature for each graph, but it may be helpful for some functions.
In many applications we will be interested in points of intersection of graphs of two functions and . Because at the intersection point functions are equal to each other, such points can be found from the equation .
The above mentioned tips are represented graphically in fig.1.13.
Finally let us formulate the main rules for graphing functions with parameters.
Plan for graphing functions with parameters
- 1
Try to simplify the function and determine if it belongs to a known class of functions with graphs from (Fig.1.2-1.11). 2. 2
Computer trail:
- (a)
Put parameters to ’reasonable’ values and plot the graph using a calculator. 2. (b)
Collect qualitative information such as : number of zeros, existence of vertical and horizontal asymptotes. 3. (c)
Vary parameter values to see how this changes the shape of the graph. 3. 3
Algebraic approach (note, not all steps may be possible):
- (a)
Find -intercept (), and zeros of the function (). 2. (b)
Find horizontal asymptote from the limit and vertical asymptote(s) (for rational function they are at the points where the denominator becomes zero ()) ( fig.1.13). 3. (c)
Find other special points (e.g. maximum, minimum, etc), if they are important determinants of the graph shape. 4. (d)
Draw the graph and indicate how the graph shape changes for different parameter values.
Example Plot the graph of the function . Find how the graph depends on the parameters and
Solution.
- 1
We do not need to simplify the function. The function equation has some similarities with graph classes listed above, but does not coincide exactly with any of them. 2. 2a
Let us put and plot the graph using calculator (fig.1.14a). 3. 2a
The graph (fig.1.14a) has the following characteristic features: the -intercept here is , we see one zero of the function . If we fill in large values of , we find that , and , thus we expect to have a horizontal asymptote , we do not see any vertical asymptotes and function has an extremum point (maximum). However, will these features persist for other parameter values? In order to answer that let us perform an algebraic study. 4. 3a
-intercept is . Zeros of the function are given by , which has only one solution , as for all . 5. 3b
The function does not have vertical asymptotes as the denominator cannot be zero ( for all and ). To find the horizontal asymptote let us compute , thus the horizontal asymptote is the -axis. 6. 3c
Important point here is the location of the maximum of the function. Let us find it. For that let us find the points where the derivative of the function is zero. The derivative of the function is , thus the expression is zero if , or . For we have just one solution . The value of the function at this extreme point is , thus the maximum is at . 7. 3d
Let us draw the graph now. Because the graph always goes through the origin. Then the graph will reach the maximum at (Fig.1.14a, symbol ’1’) and then approaches the axis. If we put all this information together we obtain a qualitative graph shown in fig.1.14b. We see that it qualitatively coincides with the calculator sketch. Now let us find how the graph shape depends on the parameter values. We see that the graph has a bell-shape, with a single maximum at . The location of this maximum depends on the parameter only, but the maximal value of the function increases if increases. The solid and the dashed line in Fig.1.14c illustrate how the graph shape changes if increases while we keep constant. Alternatively, if we keep the value constant but increase the value of the parameter (the solid and the dot-dashed line in fig.1.14c), the location of the maximum shifts to the right, and the maximal value of the function decreases.
1.4 Implicit function graphs
As we know the relation between two variables and can be expressed explicitly in terms of a function that gives us the value of if we know the value of . It is also possible that the relation between and is given implicitly as an equation. Such relations are called implicit functions, and their graphs are implicit function graphs. One of the most effective methods to plot such graphs is to try to solve that implicit equation and rewrite it as one or several explicit functions. In some cases the relations between and can be plotted directly. Let us consider two examples:
Example: Draw a graph of the function(s) given by equation:
Solution: We can either rewrite it as two explicit functions and draw the two graphs given by this equation. Alternatively, we can note that gives a square of the distance from the point to the point , thus equation gives the points located at a distance from the origin. That is a circle with a radius with the center at (Fig.1.15a). We will use this graph later in our course in chapter 4 to plot fig.4.8a.
Example: Draw a graph of the function(s) given by equation: , where are the variables and are the parameters.
Solution: Let us factor the equation:
[TABLE]
The product of two numbers is zero if one of these numbers is zero, therefore this equation is equivalent to:
[TABLE]
Graphing of is trivial. In order to graph let us rewrite it as , or . The horizontal asymptote of this graph is: . The vertical asymptote occurs if the denominator of the fraction is zero, i.e. at . However, because and this asymptote will be outside the range of our function. Additionally note that this function is similar to the graph of Fig.1.10c, but it is shifted upward by . Thus the function graph here contains two branches and that are plotted in Fig.1.15b
1.5 Exercises
Exercises for section 1.1
Perform the indicated operations:
- (a)
2. (b)
2. 2.
Find limits:
- (a)
2. (b)
3. 3.
Solve the equation for the specified variable:
- (a)
find in: 2. (b)
find in: 3. (c)
find in: 4. (d)
find in: 5. (e)
find in: 4. 4.
Solve the system of equations for the specified variables:
- (a)
find in: \left\{\begin{array}[]{l}x-2y=-5\\ 2x+y=10\end{array}\right. 2. (b)
find in: \left\{\begin{array}[]{l}ax+by=0\\ cx+dy=-b\end{array}\right. 3. (c)
find in: \left\{\begin{array}[]{l}x(1-2x)+xy=0\\ 4y-xy=0\end{array}\right. 4. (d)
find in: \left\{\begin{array}[]{l}4x-xy-x^{2}=0\\ 9y-3xy-y^{2}=0\end{array}\right. 5. (e)
find in: \left\{\begin{array}[]{l}b(1-{R\over k}-d-aN)R=0\\ (R-\delta)N=0\end{array}\right.,\;\;a,b,d,k,\delta\neq 0;
Exercises for section 1.2
Find the derivative of :
- (a)
2. (b)
3. (c)
4. (d)
Exercises for section 1.3
Without plotting the function find the following information about their graph: find the -intercept and zeros; find horizontal or/and vertical asymptotes (if they exist). (Proof of non-existence of asymptotes is not required).
- (a)
function given by 2. (b)
function given by 7. 7.
Sketch graphs of the following functions:
- (a)
2. (b)
3. (c)
. Find how the shape of the graph for depends on the value of the parameter .
Exercises for section 1.4
(a) Sketch qualitative graphs of the following implicit functions. (b) Find how special points of these graphs (intercepts, zeros, asymptotes) depend on parameters. (c) If graph contains several lines find their intersection points. Note, all parameters represent positive** numbers.**
- (a)
on the plane 2. (b)
on the plane 3. (c)
on the plane 4. (d)
on the plane 5. (e)
on the plane 6. (f)
on the plane 7. (g)
on the plane 8. (h)
on the plane
Additional exercises
Perform the indicated operations:
- (a)
2. (b)
10. 10.
Solve the system of equations for the specified variables:
- (a)
find in: \left\{\begin{array}[]{l}rA(1-{A\over K}-AB)=0\\ AB-dB-BC=0\\ BC-fC=0\end{array}\right.\;\;r,K,d,f>0 11. 11.
Find the derivative of :
- (a)
2. (b)
3. (c)
4. (d)
5. (e)
6. (f)
7. (g)
8. (h)
9. (i)
12. 12.
Solutions of differential equations
- (a)
Show that function , where is an arbitrary constant, is a solution of the differential equation: . For that, compute derivative of this function and substitute this derivative and the function itself to the equation and show that the left hand side of the equation equals to the right hand side. 2. (b)
Show, using the same steps that the function , where is an arbitrary constant and is a parameters, is a solution of the differential equation 13. 13.
Assume that is an unknown function of . For listed below find the following derivatives: and .
- (a)
2. (b)
3. (c)
Find the expression for for an arbitrary . 14. 14.
Without plotting the function find the following information about their graph: find the -intercept and zeros; find horizontal or/and vertical asymptotes (if they exist).
- (a)
function given by 2. (b)
function given by 15. 15.
Sketch graphs of the following functions:
- (a)
2. (b)
3. (c)
4. (d)
. Find how the shape of the graph depends on the value of the parameters and . 5. (e)
, Find how the shape of the graph depends on the value of the parameter and . 6. (f)
, find for which values of the parameter the graph touches the -axis. Tip: draw graph for and think about how affects this graph. 16. 16.
Sketch qualitative graphs of the following pairs of implicit functions on the same graph. Find all intersection points.
- (a)
2. (b)
Chapter 2 Selected topics of calculus
In this chapter we introduce several new notions on calculus and algebra which are important for our course.
2.1 Complex numbers
Complex numbers were introduced for the solution of algebraic equations. It turns out that in many cases we can not find the solution of even very simple quadratic equations. Consider the general quadratic equation:
[TABLE]
The roots of (2.1) are given by the well known ’abc’ formula:
[TABLE]
where
[TABLE]
What happens with this equation if ? Does the equation have roots in this case?
Complex numbers help to solve such kind of problems. The first step is to consider the equation
[TABLE]
Let us claim that (2.4) has a solution and denote it in the following way:
[TABLE]
where
[TABLE]
Here is the basic complex number which is similar to for real numbers. Using it we can denote solutions of other similar equations. For example if
[TABLE]
Similarly the equation , has solutions . Although we call a complex number, it is quite different from usual real numbers. Using complex numbers we cannot count how many books we have in the library, for example. The only meaning of is that , and is just a designation of a root of the equation .
Now we can solve equation (2.1) for the case . If , then and
[TABLE]
Example. Solve the equation
Solution.
[TABLE]
or
We see, that solution of this equation has two parts, one part is just a real number ’-1’, which is the same for and and the other part, is times another real number ’3’ which has opposite signs for and . This is a general form of representation of complex number. Any complex number can be represented in the form:
[TABLE]
where is called ** the real part** of the complex number , and is called ** the imaginary part** of . The notation for the real part is and for the imaginary part is . In our example and
We can work with complex numbers in the same way as with usual real numbers and expressions. The only thing which we need to remember, is that .
To add two complex numbers we need to add their real and imaginary parts. For example
[TABLE]
Similarly, multiplication by a real number results in multiplication of the real and imaginary part by this number
[TABLE]
Multiplication of two complex numbers is just an exercise in multiplication of two expressions
Similarly
[TABLE]
Now we can check that is a solution of the equation in example (2.8). In fact: , i.e. left hand side of this equation after substitution of equals zero and thus is the root of this equation.
One more definition. The number is called the complex conjugate to the number and is denoted as . Complex conjugate numbers have the same real parts, but their imaginary parts have opposite signs.
Roots of a quadratic equation with negative discriminant are complex conjugate to each other. It follows from the formula (2.7)
[TABLE]
hence:
[TABLE]
Finally consider two more basic operations. If , then, is called the absolute value, or modulus of . Note, that , as .
We use this trick to introduce division of two complex numbers
[TABLE]
So, to divide two complex numbers we multiply the numerator and the denominator by a number which is the complex conjugate to the denominator, and we get the answer in the usual form.
Example
[TABLE]
2.2 Matrices
From a very general point of view a matrix is a representation of data in the form of a rectangular table. An example of a matrix composed of numbers is given below:
[TABLE]
This matrix has two rows and three columns. We will call this a matrix of the size . In general matrix size is defined as . Even if you did not have matrix algebra in school, you probably know at least one matrix object, that is a vector. Indeed, a vector is an object which is characterized by its components: two numbers in two dimensions or three numbers in three dimensions. In matrix algebra vectors can be represented in two forms: as a column vector, i.e. as matrix (preferred representation), or as a row vector, i.e. as matrix. For example a vector with the -component and the -component can be represented a column or as a row vector as:
[TABLE]
Using matrices we can perform the same operations on large blocks of data simultaneously. For example, if we need to multiply all 6 numbers of matrix in (2.12) by , we can write it as which will mean:
[TABLE]
This operation is called multiplication of a matrix by a number. For a general matrix this can be written as:
[TABLE]
Similarly, addition of matrices is adding the numbers that have the same location. This operation is defined only for two matrices of the same size:
[TABLE]
For general matrices it can be written as:
[TABLE]
Multiplication of matrices is not so trivial. In general matrix multiplication is defined as the products of the rows of the first matrix with the columns of the second matrix. Thus, to fund the element in row and column of the resulting matrix we need to multiply the th row of the first matrix by the th column of the second matrix. Thus we can multiply two matrices only if the number of columns in matrix equals the number of rows in matrix .
For a product of two matrices this gives:
[TABLE]
From this it follows that multiplication of a matrix by a column vector is given by:
[TABLE]
The last equation is useful for representation of linear systems as can be seen from the following example. Assume we have a system of linear equations:
[TABLE]
we can write the coefficients at and in the left hand side as a square matrix:
[TABLE]
We also have two numbers in the right hand side which we can write as a column vector:
[TABLE]
Now if we write and as a column vector:
[TABLE]
we can represent system (2.18) using matrix multiplication (2.17) as:
[TABLE]
Indeed, from (2.17) we get \left(\begin{array}[]{lr}1&-2\\ 2&1\end{array}\right)\left(\begin{array}[]{c}x\\ y\end{array}\right)=\left(\begin{array}[]{c}1*x-2*y\\ 2*x+1*y\end{array}\right), that proves this result.
Another important matrix operation is the determinant of a square 2x2 matrix, which for the matrix A=\left(\begin{array}[]{lr}a&b\\ c&d\end{array}\right) is defined as:
[TABLE]
The determinant of a matrix has many important applications in algebra. For example using determinants it is possible to find solution of system of linear equations (e.g. system (2.18)) in the form of so-called Cramer’s rule, which was published by Gabriel Cramer as early as in mid-18th century. Cramer’s rule is briefly formulated in exercise 7 at the end of this chapter.
Now let us consider one of the most important problems in matrix algebra: the eigen value problem.
2.3 Eigenvalues and eigenvectors
Let us start with a definition:
Definition 1
A nonzero vector and number are called an eigen vector and an eigen value of a square matrix if they satisfy equation:
[TABLE]
Eigen vectors are not unique, and it is easy to see that if we multiply it by an arbitrary constant we get another eigen vector corresponding to the same eigen value. Indeed by multiplying (2.21) by we get:
[TABLE]
therefore, we can say that is also an eigen vector of (2.21) corresponding to eigen value .
For example, for matrix A=\left(\begin{array}[]{lr}1&2\\ 2&1\end{array}\right), number and vector \mathbf{v}=\left(\begin{array}[]{c}1\\ 1\end{array}\right) are an eigen value and eigen vector as:
[TABLE]
If we multiply \mathbf{v}=\left(\begin{array}[]{c}1\\ 1\end{array}\right) by any number, e.g. , , or etc., we will get new eigen vectors \mathbf{v}=\left(\begin{array}[]{c}2\\ 2\end{array}\right), \mathbf{v}=\left(\begin{array}[]{c}25\\ 25\end{array}\right) of this matrix for . You can check it in the same way as we did in (2.23) for a vector \mathbf{v}=\left(\begin{array}[]{c}1\\ 1\end{array}\right).
Finding eigen values and eigen vectors is one of the most important problems in applied mathematics. It arises in many biological applications, such as population dynamics, biostatistics, bioinformatics, image processing and many others. In our course we will apply it for the solution of systems of differential equations, which we will consider in chapter 4.
Let us consider how to solve the eigen value problem for a 2x2 matrix A=\left(\begin{array}[]{lr}a&b\\ c&d\end{array}\right). For that we need to find and \left(\begin{array}[]{c}v_{x}\\ v_{y}\end{array}\right) satisfying:
[TABLE]
We can rewrite it as a system of two equations with three unknowns :
[TABLE]
If we collect all unknowns at the left hand side we will get the following system:
[TABLE]
This system always has a solution , however it is not an eigen vector, as in accordance with the definition the eigen vector should be nonzero. In order to find non-zero solutions let us multiply the first equation by , the second equation by and subtract them. Multiplication gives:
[TABLE]
Subtraction of the equations results in:
[TABLE]
[TABLE]
as we get:
[TABLE]
This is a quadratic equation with unknown and for each particular coefficients we can find two solutions: and using the ’abc’ formula. Thus we found that the eigenvalue problem for a matrix (2.24) has solutions for the eigen values . In general, for a matrix that the eigen value problem has solutions for .
Equation (2.29) is very important in our course and it has a special name: characteristic equation. In most of the courses on mathematics this equation, however, is written in a slightly different matrix form. To derive it let us recall the definition of the determinant of a matrix given in section 2.2:
det\left|\begin{array}[]{lr}a&b\\ c&d\end{array}\right|=ad-cb. Similarly the determinant of matrix \left(\begin{array}[]{lr}a-\lambda&b\\ c&d-\lambda\end{array}\right) is:
[TABLE]
which coincides with the left hand side of characteristic equation (2.29) and thus the characteristic equation can be rewritten as:
[TABLE]
Let us use this approach to find the eigen values of matrix from example (2.23). We get the following characteristic equation:
[TABLE]
From the ’abc’ formula:
[TABLE]
therefore we found two eigen values and .
Now, let us find eigen vectors. For that let us substitute the found eigen values to the original equation (2.26) and solve it for and . Let us do it first for a particular example (2.23) for which we have found eigen values and . For eigen vector corresponding to eigen value we obtain:
[TABLE]
Both equations give the same solution . This means that if , then and a pair \left(\begin{array}[]{c}1\\ 1\end{array}\right) satisfies the system and thus gives an eigen vector of problem (2.23). We can also use any other value for . For example, if we use then will be and we get another eigen vector \left(\begin{array}[]{c}2\\ 2\end{array}\right), etc. In general any , and give an eigen vector. We can express it by the following formula:
[TABLE]
where is an arbitrary number. Formula (2.33) gives all possible solutions of eq.(2.32). It also illustrates a general property of eigen vectors which we have proven in (2.22), that if we multiply an eigen vector by an arbitrary number will get also an eigen vector of our matrix. Using this property we can formulate an easy way to write a formula for all eigen vectors. For that we take any found eigen vector and multiply it by an arbitrary number . Note, that if for problem (2.32) we use another found eigen vector \left(\begin{array}[]{c}2\\ 2\end{array}\right), we can write an answer as \left(\begin{array}[]{c}v_{x}\\ v_{y}\end{array}\right)=k\left(\begin{array}[]{c}2\\ 2\end{array}\right). At the first glance this formula is different from (2.33). However, it is easy to see that both formulas give the same result: this is because in (2.33) is an arbitrary constant and any vector given by the formula (2.33) with \left(\begin{array}[]{c}1\\ 1\end{array}\right) can be obtained using the formula k\left(\begin{array}[]{c}2\\ 2\end{array}\right) for another value of . Thus the answer to our problem: to find eigen vectors of matrix (2.23) for eigen value , can be written as \left(\begin{array}[]{c}v_{x}\\ v_{y}\end{array}\right)=\left(\begin{array}[]{c}1\\ 1\end{array}\right), or \left(\begin{array}[]{c}v_{x}\\ v_{y}\end{array}\right)=\left(\begin{array}[]{c}2\\ 2\end{array}\right), or etc. These vectors give particular solutions of this problem. We can also write a formula for all solutions as \left(\begin{array}[]{c}v_{x}\\ v_{y}\end{array}\right)=k\left(\begin{array}[]{c}1\\ 1\end{array}\right), or \left(\begin{array}[]{c}v_{x}\\ v_{y}\end{array}\right)=k\left(\begin{array}[]{c}2\\ 2\end{array}\right), or etc. As we discussed above all these answers will be correct and equivalent.
Similarly we find the eigen vector corresponding to the other eigen value :
Substitution:
[TABLE] 2. 2.
Relation between and :
[TABLE] 3. 3.
Eigen vector: use e.g. , thus
[TABLE]
The general form is \mathbf{v}=k\left(\begin{array}[]{c}-1\\ 1\end{array}\right), where is an arbitrary number.
Note, that in both cases in order to find eigen vectors we could use the first equation only (see equations (2.32) and (2.34)), and the second equation in both cases did not provide us any new information. It is not a coincidence, and this property is the basis for the following express method for finding eigen vectors:
Express method for finding eigen vectors
Let us derive a formula for finding the eigen vectors of a general system (2.25). We assume that we have found eigen values and from the characteristic equation (2.31). To find the corresponding eigen vectors we need to substitute the found eigen values into the matrix and solve the following system of linear equations (2.26):
[TABLE]
It is easy to check that if we use for and the values and it gives the solution of the first equation:
[TABLE]
If we substitute these expressions into the second equation we get:
[TABLE]
To prove that this expression is also zero, note that is zero in accordance with the characteristic equation (2.29). Therefore and give a solution of (2.35) which is an eigen vector corresponding to the eigen value . Similarly we find the eigen vector corresponding to the the eigen value .
However, this approach does not work if in (2.26) both and . In this case we can use the second equation and find an eigen vector as and . Indeed:
[TABLE]
As in the previous case it is easy to show that this vector satisfies the other (first) equation as: due to (2.35).
The final formulas are:
[TABLE]
or
[TABLE]
where are the elements of the matrix A=\left(\begin{array}[]{lr}a&b\\ c&d\end{array}\right).
Either (2.36) or (2.37) can be used to find eigen vectors. (Both answers will be valid.) If, however, one of the formulas gives a zero eigen vector, we should use the other one to obtain a non-zero vector.
Let us apply these formulas for the system (2.23) with matrix A=\left(\begin{array}[]{lr}1&2\\ 2&1\end{array}\right) and eigen values . The eigen vectors can be found from (2.36) as:
[TABLE]
and from (2.37) as:
[TABLE]
We see that the vectors differ from the vectors found earlier, but it is easy to find that they are equivalent. For example, if we multiply the first vector by we find -{1\over 2}\left(\begin{array}[]{c}-2\\ -2\end{array}\right)=\left(\begin{array}[]{c}1\\ 1\end{array}\right), thus the same vector which we found earlier in (2.33). We also see that formulas (2.38) and (2.39) give equivalent result. Indeed, first vectors obtained form (2.38) and (2.39) are the same. For second vectors note that: -1\left(\begin{array}[]{c}-2\\ 2\end{array}\right)=\left(\begin{array}[]{c}2\\ -2\end{array}\right).
2.4 Functions of two variables
A function of two variables describes the rule of finding the value of function , if we know the values of the variables and . For example, the area of a right-angled triangle with the sides , and is given by the following function of two variables: . Another example is the rate of growth of a prey population in a typical ecological predator-prey model: , where is the prey population and is the predator population. The graph of the function of one variable is a line on the -plane. To sketch the graph of the function of two variables , we must use a three dimensional space : the -plane for the values of the independent ’ input’ variables , and the third axis for the function ’output’ value . In such a representation the graph will be a surface in a three dimensional space. Fig.2.1 shows a graph of the function plotted by a computer.
Derivatives. The next step is the definition of the derivative of . The main idea of finding the derivative of is to fix one variable at a constant value, say . After that we will get a function of one variable only (). Now, we can find the derivative of , as the usual derivative of a function of one variable . For example, . Let us fix . We get the following function of one variable: . We can easily find the derivative now: .
This type of derivative is called the partial derivative of with respect to at . We denote it as
[TABLE]
We can find such a derivative at , or at any other value of . In fact for an arbitrary , , and
[TABLE]
Here as we replaced by a constant and the derivative of a constant is zero. Similarly, , and , as is a constant and the derivative of . It is generally accepted to make all these differentiations without explicitly replacing by . We just should keep in mind, that for such a differentiation we treat as a constant. Thus, to find the derivative of with respect to we just write:
[TABLE]
keeping in mind that is considered as a constant and not a variable during this differentiation.
This expression is called the partial derivative of with respect to and is denoted as .
Similarly, we can introduce a partial derivative of with respect to : . To compute it, we fix (treat as a constant) and make the usual differentiations with respect to . In our example it gives:
[TABLE]
Here , , and as is fixed.
Example. Find and for
Solution , as for we fix , and . Similarly, , as and hence is treated as a constant.
The geometrical representation of a partial derivative is clear from fig.2.2. To compute we fix , i.e. assume that has some value . The condition geometrically gives a horizontal line on the plane fig.2.2a, or a line parallel to the -axis. In 3D this line gives a curve on the 3D surface in graph fig.2.2b, which is a 1D function. The partial derivative with respect to for this particular will give us the slope of the tangent line to this 1D function. Thus (see fig.2.2) gives the slope of the tangent line in the direction of the -axis or the rate of change of in the direction. Similarly, computing we fix , i.e. assume that has some value . It gives us a vertical line on the plane fig.2.2a, or a line parallel to the -axis. Thus gives the slope of the tangent line in the direction of the -axis, or the rate of change of in the direction. If we consider as a mountain gives the slope of the mountain if we climb in the -direction and gives the slope of the mountain if we climb in the -direction.
Note, that in general at each point on a surface we can draw a tangent line in any direction, and partial derivatives and give the slopes of two of these possible tangent lines. Note, that the slope of a tangent line any direction can be obtained as a combination of these two slopes.
Linear approximation Let us derive a formula for approximating a functions of two variables . Let us assume that we know and its partial derivatives at some point and that we want to find the value of a function at the close point (fig.2.3).
Let us move to the point in two steps. Let us first move from the point to the point , i.e. in the -direction, and then from to , i.e. in the -direction. Let us apply the formula for approximation of a function of one variable in formulation (1.10) at each part of this motion. Because on the first part we move along the direction the change of the function () will be given as the product of the rate of change of the function in the direction () times the distance between the points ():
[TABLE]
Similarly, on the second part of our motion, we move along the -axis, and the change of the function here () will be given as the product of the rate of change of the function in the direction () times the distance between the points ():
[TABLE]
If we add equations (2.40) and (2.41) and solve it for we find the following formula which gives the approximation for a function of two variables:
[TABLE]
This expression is called a linear approximation, as the independent variables are in the first power only, we do not have terms , or , or etc.
Example Find the linear approximation for the function at the point
Solution. We use the formula (2.42) with and .
; at
, at
Finally, .
At the approximate formula gives: . The exact value of
2.5 Exercises
Exercises for section 2.1
Find all roots of the given equations
- (a)
2. (b)
Exercises for section 2.2
Write the following linear systems in a matrix form . Find the determinant of matrix .
- (a)
\left\{\begin{array}[]{l}2x-4y=3\\ x+y=1\end{array}\right. 2. (b)
\left\{\begin{array}[]{l}ax+by=0\\ cx+dy=-b\end{array}\right.
Exercises for section 2.3
Find eigen values and eigen vectors of the following matrices:
- (a)
\left(\begin{array}[]{lr}-2&1\\ 1&-2\end{array}\right) 2. (b)
\left(\begin{array}[]{lr}1&4\\ 1&1\end{array}\right) 3. (c)
\left(\begin{array}[]{lr}-1&5\\ -1&3\end{array}\right)
Exercises for section 2.4
Find partial derivatives of these functions. After fining derivatives evaluate their value at the given point (if asked).
- (a)
and for at 2. (b)
for at 3. (c)
and for at and at 4. (d)
and for . 5. (e)
and for 6. (f)
and for 7. (g)
and for 8. (h)
and for
Additional Exercises
Perform the indicated operations:
- (a)
2. (b)
3. (c)
4. (d)
6. 6.
Find all roots of the given equations
- (a)
2. (b)
7. 7.
Cramer’s rule on an example. Cramer’s rule for system of two linear equations:
[TABLE]
allows us to find solutions from determinants of matrices. First we need to find a determinant of the main matrix
[TABLE]
Then we need to find the determinant of a matrix formed by replacing the -column values of the matrix with the answer-column values \left(\begin{array}[]{c}E\\ F\end{array}\right) as detD_{x}=det\left|\begin{array}[]{lr}E&b\\ F&d\end{array}\right| and similarly for the -column: detD_{y}=det\left|\begin{array}[]{lr}a&E\\ c&F\end{array}\right|. The solutions of the system will be given by the ratios of these determinants as:
[TABLE]
.
- •
Find solutions of the following system using the Cramer’s rule:
[TABLE]
- •
Find also solution of (2.43) by usual method. Show that Cramer’s rule gives a correct result. 8. 8.
Find eigen values and eigen vectors of the following matrices:
- (a)
\left(\begin{array}[]{lr}-1&6\\ 2&-2\end{array}\right) 2. (b)
\left(\begin{array}[]{lr}2&1\\ 7&-4\end{array}\right) 9. 9.
Find a linear approximation for the function at the given point.
- (a)
at
Chapter 3 Differential equations of one variable
Differential equations are equations that contain a derivative of an unknown function. As we know derivative gives a velocity of a process and differential equations occur when we describe various processes via their velocities. Differential equations are widely used for modeling in a variety of disciplines: in mathematics, physics, chemistry, economics, engineering, medicine, life sciences, etc. Development of methods of study of differential equations is the main subject of this course.
In this chapter we will introduce differential equations, give first definitions, show how to solve simple differential equations analytically and consider a few examples. Then we will develop qualitative methods for the analysis of differential equations of one variable and will apply them for biological models.
3.1 Differential equations of one variable and their solutions
3.1.1 Definitions
Let us construct a first differential equation. Consider a motion of a car with a contact velocity, for example . If we denote the distance traveled by the car at time as we can write , as the velocity is the derivative of the distance with respect to time. Thus we can write the following differential equations for this process:
[TABLE]
If the car travels with an acceleration, then the velocity will linearly increase with time. If we assume that the acceleration is , then the velocity in the course of time will be given by and we get a differential equation as:
[TABLE]
Let us consider a biological example. If is the population size of a species at time , then the rate of change of the population size is:
[TABLE]
Let us assume that the birth and the death terms are proportional to . This assumption is quite reasonable, as it means, that if, for example, we know the growth rate of a population of some insects on one tree, then the total growth rate of the whole population of insects in the forest will be proportional to the total number of trees. Thus each term in equation (3.3) will be proportional to and we get the following famous ’Malthus’ equation for population dynamics:
[TABLE]
where and are the rate constants for the birth and death processes, and we see that if , and , if .
Another model assumes that there is a maximum population size (called the carrying capacity) and that the the rate of growth of population depends on how close the population is to this maximum size. This yields the following differential equations:
[TABLE]
Now note, that in mathematical sense equations of the type (3.2) can be written as
[TABLE]
as the variable with respect to which we differentiate function is also present in the right hand side of our equation.
Alternatively the equations of the type (3.4) and (3.5) can be written as:
[TABLE]
as here the variable is not present at the right hand side and we have only the unknown function (for eqns(3.4),(3.5) ) there.
The latter equations are the most important for us in this course and they have a special name an ’autonomous differential equations’:
Definition 2
Equation
[TABLE]
is called an autonomous differential equation
Example
Before we find how to solve differential equations, let us discuss from a very general point of view what kind of solutions can we expect here. If we assume that a differential equation describes how the size of a population will change in time, then we may think about two types of problems. The first one it to find this size for a particular population at each time moment. For that we obviously would need to know the initial size of a population. We can also ask a more general question: to find the population size for an arbitrary initial size. This solution will be called the general solution of a differential equation. Because the general solution contains information on solutions for arbitrary initial conditions, it normally depends on an arbitrary constant. The differential equation with given initial condition is called an initial value problem:
Definition 3
The problem is called the initial value problem; Its solution is called the orbit or trajectory.
The initial value problem for most has a unique solution.
Now let us consider how we can solve differential equation.
3.1.2 Solution of a differential equation
e can easily solve an equations of type (3.6) using the method of separation of variables. The main idea of this method is to think about the derivative of an unknown function in as a fraction over . If we multiply both sides of this equation by we will get:
[TABLE]
If we integrate both sides of this equation and get:
[TABLE]
where is an anti-derivative of and an arbitrary constant. This is the general solution of differential equation (3.6). We see that this solution contains an arbitrary constant, as expected from the discussion in the previous section.
Let us apply this method to a few examples.
For the simplest differential equation we get:
[TABLE]
Therefore, the solution of this equation is equals any constant .
For the equation of motion of a car with a velocity (3.1) we get:
[TABLE]
Thus this solution shows us the position of a car as a function of time , and the arbitrary constant here accounts for the initial position of the car.
Finally, for the motion from the rest with an acceleration (3.2) we get:
[TABLE]
We obtained a formula which is well known to you from your school physics and here also accounts for the initial position of a car.
It turns out that we can apply the method of separation of variables also for an autonomous equation (3.8) (). However, here we will need to separate variables as and as a result we will not usually get the explicit formula for . However, in many cases we can do it after some transformations.
Let us solve equation for population dynamics (3.4).
[TABLE]
To find note that equation has a solution , thus we find and if we denote we will get:
[TABLE]
where is an arbitrary constant.
Let us apply it for the following initial value problem with and the initial population size of :
[TABLE]
The general solution here is given by . To find the solution of the initial value problem we note that at the population size was , i.e. we can write: , or we find that Hence, we have the following solution of the initial value problem (3.15): .
Similarly, we can solve the general initial value problem (3.4) for an arbitrary initial size of and find:
[TABLE]
thus here gives the initial size of the population.
Finally let us solve equation (3.5), for specific parameter values :
[TABLE]
This is a general solution. Let us find a particular solution, corresponding to the initial size of the population , for example. For that we find: , or , thus the population dynamics in the course of time will be given: .
The solutions of equation (3.4) and equation (3.5) are shown in fig.3.1.
We see that the size of the population for eq.(3.4) goes to infinity (fig.3.1a). Quantitatively the size of population increases in times each seconds. Indeed, from the particular solution we find In general, for arbitrary in (3.4) this characteristic time of change is given by , which follows from equation (3.16). For equation (3.5), the population approaches its carrying capacity value of (see fig.3.1b) and the characteristic time is also determined by the value of in the following sense: the difference between the current population size and its stationary value decreases in times over the time period . This follows from solution (3.17), which gives for this difference . Similar expression for an arbitrary value of is given by , which gives for the characteristic time .
This concludes our analytical study of differential equations. In the next chapter we will formulate another method of analysis of differential equations that does not require direct integration of these equations.
3.2 Qualitative methods of analysis of differential equations of one variable
In this section we will consider a general non-linear differential equation and develop an effective method for qualitative analysis of this equation without finding solutions analytically. In the next chapters this method will be extended to the systems of two differential equations.
3.2.1 Phase portrait
Let us start with equation (3.18) which we considered in section 3.1.
[TABLE]
We found that is a solution for this equation for and the initial population size of (see (3.15)). This solution was represented graphically in fig.3.1.
If the initial size of the population was different, for example or , etc., we get other solutions of equation (3.18) and if we represent these solutions graphically we will obtain the following curves () shown in fig3.2. Let us analyze them. An important characteristic of any line is its slope. It turns out that we can easily find the slope for the solutions of (3.18): . We see that the slope depends only on (the size of the population) and does not depend on other factors, for example on the initial conditions. For example, if the slope of the line representing solution at point is 4*3=12 for any initial condition. Geometrically this means if we graph several solutions (as in fig3.2) and determine slopes of these functions for given (at points of intersection of the dotted line with the solution curves) we will find that all slopes are the same ().
We can use this information and represent a qualitative picture of solutions of (3.18) using only one -axis. For that let us denote the slope of the solution on the -axis for each value of (see numbers 4,8,12,16, etc.). We see, however, that this is not very helpful for representation of the solution of our equation.
To improve it let us think about biological interpretation of the slope of the curve in fig3.2. The slope of the curve gives us the rate of change of the function and because fig3.2 shows how the size of population depends on time, the slope values (represented above the -axis) show the growth rate of population at given . The most important qualitative aspect of the dynamics here is the gowth of population. We can represent it in the following qualitative way:
i.e. we show the growth of the population by an arrow which is directed to the right. Of course, this is not a complete description of our system, but it gives a good idea about the behavior of our system. It shows that if we start at some initial value of , then will grow and the size of the population will be continuously increasing. Note that to obtain this result we have only used the direction of the arrows in the figure.
As we will see in the next section, such representation can be easily obtained for any autonomous differential equation from the graph of the function at the right hand size. Such representation is called the phase portrait:
Definition 4
The collection of all possible orbits of a differential equation together with the direction arrows is called the phase portrait.
3.2.2 Equilibria, stability, global plan
Let us consider two differential equations that arise in population ecology:
[TABLE]
and
[TABLE]
Let us sketch phase portraits for (3.19) and (3.20). We can do it without finding a solution. In general, to sketch a phase portrait of an equation we need to draw or arrows on the -axis. The arrow means growth of , i.e. . The arrow means decreasing of , or . Because we will graph function and then assign to that regions where the graph is above the -axis and to that regions where the graph is below the -axis. For equation (3.19) , the graph of is shown at the top panel of fig.3.3a and we draw the right arrow for , and the left arrow for (fig.3.3a bottom).
The interesting point here is . Here , thus the rate of change of here is zero and we cannot assign any direction for the arrow at this point. However, the dynamics of our system here is trivial: do not change in the course of time, or for all . This means, that if the initial size of the population was zero, it will be zero forever. Such points of a phase portrait are called equilibria. They occur at points where the rate of change is zero () For equation equilibria occur if , which is also used as a definition of an equilibrium.
Definition 5
A point is called an equilibrium point of , if
Finally the phase portrait of eq.(3.19) in fig.3.3a gives the following dynamics of : if the initial value of is to the right or to the left left of the equilibrium point , it will go to plus or minus infinity respectively.
Let us study equation (3.20). Again, our plan is: (fig.3.3b). Here we have an equilibrium point which is the root of the function , and the arrows (flow) for this case are shown in fig.3.3b. So, the dynamics of solutions of our equation will be as in the following figure:
i.e. for any initial condition, will eventually approach the equilibrium point.
If we compare the equilibria of equations (3.19) and (3.20), we can see that they are different. The variable diverges from the equilibrium point of equation (3.19). Such equilibrium points are called non-stable equilibria. On the contrary, the variable converges to the equilibrium point of equation (3.20). Such equilibria points of are called stable equilibria or attractors.
Now we can formulate a general plan for finding the phase portrait of .
Global plan.
Sketch the graph of . 2. 2.
Draw the phase portrait. For that transform the points where to equilibria points, the regions where to right headed arrows (), and the regions where to left headed arrows (). This gives the overall phase portrait.
Let us apply it to the logistic equation for population growth
[TABLE]
This equation describes growth of a population in a medium with limited resources. We can study (3.21) for arbitrary values of parameters . However for simplicity let us fix and . The equation becomes
[TABLE]
Let us find the phase portrait and the dynamics of the solutions of (3.22). First we use the global plan.
The right hand side of our equation is . The graph of this function is a parabola, opened below with the roots .
The construction of the phase portrait is clear from fig.3.5a.
The dynamics of the system for different initial conditions is shown in fig. 3.5b, for an initial size of the population , and in fig.3.5c for . We see that in the course of time the size of the population becomes , i.e. the stable equilibrium point is the only attractor of our system.
Sometimes, differential equations have several stable equilibria (attractors). For example, the model for the spruce bud-worm population (3.23) has the following phase portrait (fig.3.6).
[TABLE]
We see that there are two attractors: A1 and A2 which correspond to bud-worm populations of different size. We see that if the initial size of the population is , then the population eventually reaches A1; if , then population eventually reaches A2. These intervals are called basins of attraction.
Definition 6
The basin of attraction of a stable equilibrium point is the set of values of such that, if is initially somewhere in that set, it will subsequently move to the equilibrium point .
In the case of fig3.6, the basin of attraction of the equilibrium (A1) is the interval ; the basin of attraction of the equilibrium (A2) is the interval . It is very important to know basins of attraction of a system in order to predict its behavior.
3.3 Systems with parameters. Bifurcations.
One of the aims of modeling in biology is to predict the behavior of a system for different conditions. In that case we differential equations will contain parameters. Let us consider two examples. The first is a general linear equation with one parameter
[TABLE]
If we draw the graph of the right hand side function for different values of the parameter we find that we can have two possibilities depending on the sign of the parameter (fig.3.7): a non-stable equilibrium point at for (fig.3.7a) and a stable equilibrium point at for (fig.3.7b).
Another example of an equation with parameters is the logistic equation for population growth (3.21). This equation depends on two parameters and , where accounts for the growth rate and accounts for the carrying capacity. Let us consider a slight modification of equation (3.21) for a population which is subject to harvesting at a constant rate :
[TABLE]
where is the harvesting rate (an extra parameter).
Let us fix the parameters and (the same values as in equation (3.22)), and study only the effect of varying the harvesting parameter on, the dynamics of the population.
[TABLE]
When equation (3.26) coincides with equation (3.22) which was studied in fig.3.5. Now assume that the harvesting is not zero. Let us plot graphs of for (fig.3.8a).
The phase portraits for are shown in fig3.8b. We see that at the behavior of the system is qualitatively similar to the behavior of the system without harvesting (): the population eventually approaches the stable non-zero equilibrium. However the final size of the population in this case is slightly smaller than for the population without harvesting (fig.3.5a) At the situation is different. We do not have a stable and non-stable equilibrium anymore. The flow is always directed to the left and the size of the population decreases. In this simple model this means the extinction of the population. The important question here is, what is the maximal possible harvesting rate at which the population still survives. From the previous analysis it is clear that the critical harvesting is reached when the parabola touches the -axis (fig.3.8c). To find this critical value we note, that just shifts the parabola downward. Therefore, the situation of (fig.3.8c) occurs, when the shift equals the maximum of the parabola . To find the maximal value we find a point where the derivative
[TABLE]
So the maximal value of equals and therefore the maximal harvesting is
If there are no equilibria and the population will go extinct. If there is a stable and non-stable equilibrium and the population will go to the stable equilibrium. At we are at a boundary between these two qualitatively different cases. Such a qualitative change in system behavior is called a bifurcation. Bifurcations are studied in a special section of mathematics: theory of dynamical systems.
3.4 Exercises
Exercises for section 3.1
Assume that a population grows in accordance with the following equation:
[TABLE]
If the initial size of the population was , find what will be size at time . Find at what time the population will double its initial size. 2. 2.
A bacterial population doubles its size each 20 minutes. The growth of this population satisfies the differential equation . Find the value of in .
Exercises for section 3.2
Study the listed differential equations by answering the following questions:
- •
Draw the phase portrait.
- •
How many equilibria do we have here? At which ?
- •
For each equilibrium tell whether is it stable or non-stable
- •
What will be the final value of if . (e.g. converges to equilibrium, or goes to infinity, etc.)
- •
List attractor/attractors and determine their basin/basins of attraction.
- (a)
2. (b)
3. (c)
4. (d)
5. (e)
with the following graphs of :
- (f)
(this equation (Adler 1996) describes the dynamics of population of species which cannot bread successfully when numbers are too small or too large) 4. 4.
The following equation describes the production of a gene product with concentration :
[TABLE]
Here is the parameter accounting for a chemical which produces the gene product. The initial state of the system was: and the value of . At some moment of time the value of was slowly increased from to and then slowly decreased back to the value .
- (a)
The value of function at its local minimum is . *(Optional) Show this from function derivative, or using calculator. * 2. (b)
What will be the value of the concentration of the gene product at the end of the described process for ? 3. (c)
The same for ? 4. (d)
Show that there is a critical value of that separates different outcomes of this process. Find this critical value of . 5. 5.
Consider a model population with logistic growth which is subject to harvesting at a constant rate
[TABLE]
Find the maximal yield .
Additional Exercises
Assume that the growth of a mass of an animal can given by the following differential equation:
[TABLE]
where W is the weight in grams and time in weeks.
- (a)
Find the solution for the initial . 2. (b)
At what time the mass will reach half of the saturated value. 3. (c)
If we assume that the linear size of the animal is proportional to the cubic root of the mass, find at what time the object will reach half of its saturated linear size. 7. 7.
The dynamics of the ionic channels in the famous Hodgkin-Huxley model for a nerve cell is described by the following type of equations:
[TABLE]
where is a gating variable and are the parameters.
Find the steady state values of the gating variable and the characteristic time of approaching this steady state. 8. 8.
Consider a model where the harvesting is proportional to the size of the population:
[TABLE]
Find the maximal yield. 9. 9.
Compare the harvesting strategies (3.28) and (3.29). Which strategy is better. Why?
Chapter 4 System of two linear differential equations
Many biological systems are described by several differential equations. One of the most simple types of such systems is a system of two linear differential equations that on a general form can be written as:
[TABLE]
here and are unknown functions of time , and are constants (parameters).
System (4.1) by itself has many practically important applications, for example compartmental models in biology and pharmacology, electrical circuits in physics, models in economics, etc. System (4.1) will also be very important for study so-called nonlinear system of differential equations which is widely used in theoretical biology and will be considered in chapter 5.
This chapter we will introduce main definitions for linear systems (phase portrait and equilibria points) we will derive a formula for general solution of this system and classify possible solutions of this system and their phase portraits. These results will be used later to study models of biological processes.
4.1 Phase portraits and equilibria
Let us consider an example of system (4.1) with particular values for the parameters :
[TABLE]
Let us first solve this system on a computer. For that we need to choose initial values for and and let the computer find their dynamics in the course of time. Solutions for are shown in fig.4.1. We see, that in the course of time, and approach the stationary values
Let us represent this solution graphically. For a differential equation with one variable ) we presented the solutions in terms of a one-dimensional phase portrait using the -axis. For system with two variables, we need to use two axes to represent the dynamics. Let us consider a two dimensional coordinate system with the -axis for the variable and the -axis for the variable y. Such a coordinate system is called a phase space. Let us represent the trajectory from fig.4.1 on the -plane. The initial sizes of the populations were , thus we put this point () on the -plane. At the next moment of time we get other values for and and we also put them on the -plane and the and the coordinates of the next point, etc. Finally, we will get the line shown in fig.4.2a that starts at () and ends at (). To show how and change in the course of time we draw an arrow as in fig.4.2a. This trajectory is the first element of the phase portrait. If we start trajectories from many different initial conditions we will get the complete phase portrait of system (4.2) (fig.4.2b). Each trajectory represents a certain type of dynamics of , which can be easily shown on time plots similar to fig.4.1. The phase portrait in fig.4.2b give us the overall qualitative dynamics of our system: the variables and approach from all possible initial conditions.
The main aim of our course is to develop a procedure of drawing a phase portrait of a general system of two differential equations without using a computer which will allow us to study models of biological processes. In the 1D case the phase portrait consisted of two main elements: equilibria points and flows (trajectories) between them. Similar elements also compose the phase portrait of a system of two differential equations. Let us start with the first main element and define equilibria of the system.
In the 1D case equilibria were points where our system is stationary: placed at an equilibrium point the system will stay there forever. Mathematically equilibria for eq. (3.8) were determined as the points where , i.e. where . In the 2D case it is required that at the equilibrium point both variables and do not change their values, i.e. both and . For system (4.2) these conditions give the following the system of two algebraic equations for finding equilibria:
t
[TABLE]
that have the only solution . Therefore system (4.2) has an equilibrium point . As we see in fig.4.2b this equilibrium is an attractor for all trajectories.
For a general linear system (4.1) the equilibria will be given by:
[TABLE]
This system always has a solution and thus the general linear system (4.1) always has an equilibrium at the point . In the next sections we will find out how to sketch a phase portrait of (4.1) around this equilibrium. Our plan will be the following. We will first find the general analytical solution of this system and then will use it to draw the phase portraits.
4.2 General solution of linear system
Consider a general system of two differential equations with constant coefficients:
[TABLE]
The general solution of (4.5) can be written in the following form
[TABLE]
where are eigen values of the matrix A=\left(\begin{array}[]{lr}a&b\\ c&d\end{array}\right), and \left(\begin{array}[]{c}v_{1x}\\ v_{1y}\end{array}\right), \left(\begin{array}[]{c}v_{2x}\\ v_{2y}\end{array}\right) are the corresponding eigen vectors.
We will not derive the formula (4.6), we will show the main ideas behind the real derivation. For that let us consider first the one dimensional case, and then extend it to a two dimensional system.
The easiest way to find solutions of system (4.5) is by the method of substitution. Let us illustrate this method on example of 1D analog of system (4.5) which is 1D linear differential equation:
[TABLE]
We can easily solve (4.7) using the direct method of separation of variables and subsequent integration. However, let us find the solution using another method, the method of substitution. The main idea of this method is to look for a solution in some known class of functions which should be chosen in advance from some preliminary analysis. It was found that for linear systems this class is the class of exponential functions . Important questions such as: how was this class found and is this class unique etc, will not be discussed. Our aim here will be an illustration of the main components of the solution rather than comprehensive analysis of linear systems, which is a large special section of mathematics. Once the class of functions is given (in our case the class of exponential functions), we need to check under which circumstances it will satisfy the equations we are solving. Thus we will look for a solution of (4.7) of the form , where and are unknown coefficients. The main idea of the method of substitution is to find these unknown coefficients for a particular system. Let us substitute into (4.7). We find: , or:
[TABLE]
We can cancel and , and we get:
[TABLE]
Hence we found the following the solutions of (4.7):
[TABLE]
where is an arbitrary constant.
Now, let us use the same approach for the two dimensional system (4.5). It turns out that the class of functions in two dimensions will be the same as in one dimension, but because we have two variables, we need to introduce different constants for and , so our substitution will be
[TABLE]
where are unknown coefficients. Let us make this substitution for a particular example:
[TABLE]
Substitution for gives:
[TABLE]
we can cancel , (but not as in one dimensional case), and get:
[TABLE]
or in the matrix form:
[TABLE]
We see that to find the solution of (4.11) we need to solve the problem (4.14), which is the eigen value problem considered in section 2.3. To solve it we find eigen values from the characteristic equation (2.31): Det\left|\begin{array}[]{lr}1-\lambda&4\\ 1&1-\lambda\end{array}\right|=(1-\lambda)*(1-\lambda)-4=\lambda^{2}-2\lambda-3=0\;\;\lambda_{12}={2\pm\sqrt{2^{2}+3*4}\over 2}=1\pm{1\over 2}\sqrt{16}=1\pm 2. Hence
We use the express formula (2.36) for eigen vectors and get:
[TABLE]
Formula (4.15) gives just one eigen vector for each eigen value. We also know (see formula (2.33 in section 2.3) that all eigen vectors can be found by multiplication of this eigen vector by an arbitrary constant. If we denote by the constant for the first eigen vector and by the constant for the second eigen vector we will get the following solution of the eigen value problem (4.14):
[TABLE]
If we substitute these eigen vectors into the formula (4.10) we find the following solutions of (4.11)
[TABLE]
Now let us prove the following property of system (4.5):
If and are two solutions of (4.5), then is also a solution of (4.5).
Proof: As and are the solution this means that they satisfy (4.5), i.e.
[TABLE]
If we add equations for and we get: , which can be re-written as: . If we do the same for equations for and we will finally get:
[TABLE]
which explicitly shows that is a solution of (4.5).
If we apply this result for two found solutions (4.17) of (4.11) we can conclude that the sum of these two solutions is also a solution of (4.11):
[TABLE]
We proved formula (4.6) for a particular system. If we apply the same steps for a general system (4.5) we will get the general solution in the form (4.6).
So, we solved system (4.5). In the next sections we will find out how to draw its phase portraits.
4.3 Real eigen values. Saddle, node.
As we know, the general solution of (4.1) is given by (4.6). We can use this formula to sketch a phase portrait of this system. It turns out, that we can have several different types of equilibria depending on the eigen values . As we know are the roots of the characteristic equation (2.31), which is a general quadratic equation. Therefore, the roots can be real or complex numbers. In this section we consider the case of real roots.
So, assume that the eigen values and are real. This yields the following three cases:
Eigen values have different signs (, or ). 2. 2.
Both eigen values are positive () 3. 3.
Both eigen values are negative ()
Note, that we do not consider the case when . This situation is quite rare and is not considered in this course.
4.3.1 Saddle; , or
System (4.11), which we solved in section 4.2, had eigen values . Let us draw its phase portrait. The general solution of this system is given by (4.20)
[TABLE]
Because are arbitrary constants let us consider three simple cases, in which one of these constants, or both of them are zero.
-
If , then and do not depend on time. The trajectory is just one point (0,0), which is the equilibrium point of the system (4.11).
-
If any number, then \left(\begin{array}[]{c}x\\ y\end{array}\right)=C_{2}\left(\begin{array}[]{c}2\\ 1\end{array}\right)e^{3*t}. Because can change from 1 (at ) to any infinitely large number and is an arbitrary positive or negative number, this expression can be rewritten as:
[TABLE]
where is an arbitrary number from and is a vector \left(\begin{array}[]{c}-4\\ -2\end{array}\right). Thus expression (4.22) means multiplying of the eigen vector by an arbitrary number . In general, if we multiply a vector by a positive number we get a vector with the same direction but the length will be increased times. If we multiply the vector by a negative number, the direction of the vector will be changed to the opposite and the length will be changed by a factor . Because in (4.22) assumes all values from , this will give a straight line along this vector \left(\begin{array}[]{c}2\\ 1\end{array}\right) (fig.4.3a). To complete drawing the trajectory we need to show arrows indicating the motion of a point along the trajectory in the course of time. Because time dynamics is given by , in (4.22) will grow in the course of time, i.e. it will become either more positive, or more negative depending on its initial sign. Geometrically this means that a point will move apart from the origin of the coordinate system and we will get a picture as in fig.4.3b.
- The third case is any number, . The solution in this case is
[TABLE]
As in previous case we conclude, that all trajectories in this case will be located on a straight line along the vector \left(\begin{array}[]{c}-4\\ 2\end{array}\right) and we just need to show the direction of flow along this line. In this case, time dynamics is given by the function , which approaches zero when goes to infinity. Therefore, independent of initial conditions (independent of the value of ) we will approach the point , and the arrows will have the following direction (fig.4.4).
Finally let us draw the phase portrait for arbitrary and (fig.4.4b). Let us consider one trajectory for which , for example: . The solution in this case is given by:
[TABLE]
This trajectory starts as point . In the course of time will become smaller and smaller, while will grow. So, the first term in (4.24) 0.1\left(\begin{array}[]{c}-4\\ 2\end{array}\right)e^{-t} will be small compared to the second term 0.1\left(\begin{array}[]{c}-4\\ -2\end{array}\right)e^{3t}, and the dynamics at large will be described by an approximate formula : \left(\begin{array}[]{c}x\\ y\end{array}\right)\approx 0.1\left(\begin{array}[]{c}-4\\ -2\end{array}\right)e^{3t}, thus the trajectory will approach the line (4.22) presented in fig.4.3b. There will be similar behavior for any other trajectories: independently on staring points they will approach line of fig.4.3b from various directions. The qualitative picture will be as in fig.4.4b.
Such a phase portrait is called a saddle point. It has the following important features: (1) There is an equilibrium point at . (2) There are two lines associated with eigen vectors of our system (fig.4.3b,fig.4.4a). These lines are called manifolds. The manifolds in fig.4.3b and fig.4.4a are different. If we follow the trajectory along the manifold in fig.4.3b the distance to the equilibrium increases (see fig.4.4). On the contrary, if we follow the trajectory along the manifold in fig.4.4a the distance to the equilibrium decreases. The manifold from fig.4.3b is called a non-stable manifold. The manifold from fig.4.4a is called a stable manifold.
Conclusions of this study can be easily generalized. If we consider an expression:
[TABLE]
it will obviously determine a manifold (straight line) along the vector \left(\begin{array}[]{c}v_{x}\\ v_{y}\end{array}\right) and the stability of this manifold will be determined by . There are two main types of behavior of the function (fig.4.5). If approaches zero, when increases. If grows to infinity with increasing . Hence, if , eq.(4.25) will determine a stable manifold: will approach [math] in the course of time (as in fig.4.4a). If , then will diverge to infinity and we will get a non-stable manifold.
Conclusion 1
The equation (4.25) on a phase portrait gives a manifold in the form of a straight line. This line goes through the origin and is directed along the vector \left(\begin{array}[]{c}v_{x}\\ v_{y}\end{array}\right). This manifold is stable if and non-stable if .
Finally the formal definition of a saddle point:
Conclusion 2
Fig.4.4 shows the phase portrait of a saddle point. It occurs close to equilibrium, at which eigen values of the system are real and have different signs, i.e. , or . The phase portrait of a saddle point has two manifolds directed along the eigen vectors. One manifold is stable (corresponding to the negative eigen value of the system). The other manifold is non-stable (corresponding to the positive eigen value of the system).
4.3.2 Non-stable node;
Let us draw the phase portrait for the case when eigen values are real and are both positive. The general solution of the system is given by (4.6):
[TABLE]
From the previous analysis we immediately conclude, that the phase portrait in this case has the equilibrium point at and two unstable manifolds along the vectors \left(\begin{array}[]{c}v_{1x}\\ v_{1y}\end{array}\right) and \left(\begin{array}[]{c}v_{2x}\\ v_{2y}\end{array}\right) (fig.4.6a). Let us put that on the graph (fig.4.6b).
To complete the picture we need to add several trajectories which start between the manifolds. For such trajectories and both terms in (4.26) will diverge to plus or minus infinity. So, we get trajectories as in fig.4.6b. Such an equilibrium is called a non-stable node.
Conclusion 3
If the eigen values of system (4.1) are real and both positive () we have an equilibrium point called a non-stable node. To draw a phase portrait at this equilibrium we need to show two non-stable manifolds along the eigen vectors of system (4.1) and add several diverging trajectories between the manifolds.
4.3.3 Stable node;
The general solution in this case has the same form (4.26). The phase portrait will be similar to fig.4.6, but because both manifolds will be stable. So we get a picture fig.4.7a
If the trajectory starts between the manifolds it will also approach equilibrium as both terms in (4.26) will converge to 0, because (fig.4.7b). Such an equilibrium is called a stable node.
Conclusion 4
If the eigen values of system (4.1) are real and both negative () we have an equilibrium point called stable node. To draw a phase portrait at this equilibrium we need to show two stable manifolds along the eigen vectors of system (4.1) and add several trajectories converging to the equilibrium (0,0).
4.4 Phase portraits for complex eigen values: spiral, center
In the previous section we have studied the case when the roots of the characteristic equation (2.31) are real, and found three possible types of phase portrait (equilibria): saddle, stable node and unstable node. Here we will study the case when the roots of the characteristic equation (2.31) are complex.
4.4.1 General ideas on equilibria with complex eigenvalues
From section 4.2 we know that in order to find the type of phase portrait of linear system we need to solve the characteristic equation (2.31) which is a general quadratic equation, and because the parameters of the linear system are arbitrary the coefficients of the characteristic equation are also arbitrary. Therefore, it may happen that the discriminant of this quadratic equation will be negative and we will have complex eigen values. From equation (2.7) we know that these eigen values will be given by
[TABLE]
or if we denote the real part of these complex numbers as : and the imaginary part as we can rewrite (4.27) as
[TABLE]
Which type of dynamics do we expect here. If we drop for a while the imaginary part in (4.28), we get that , i.e. both eigen values will be the same and hence they will have the same sign. Which type of equilibria do we have for two real eigen values which have the same sign? We can have either a non-stable node if , or a stable node, if (see section 4.2). In the case of a non-stable node all the trajectories diverge from the equilibrium (fig.4.6b), while in the case of a stable node all the trajectories converge to it (fig.4.7b). It turns out that we will get similar behavior for complex eigen values: if we will get divergence of the trajectories, if we will get convergence of the trajectories to the equilibrium. However, the picture will be slightly different from fig.4.6b and fig.4.7b, as the imaginary part of the eigen values will add rotation to the trajectories.
The complete derivation of the formula for the general real solution of system (4.1) with complex eigen values is given in the appendix at the end of this chapter for extra reading. Here we will just demonstrate that adds rotation.
4.4.2 Center, spiral
Let illustrate that the imaginary part of the complex eigen value of the characteristic equation results in a rotation of the trajectories. For that let us consider a system:
[TABLE]
In this case the eigen values are given by: Det\left|\begin{array}[]{lr}-\lambda&2\\ -2&-\lambda\end{array}\right|=\lambda^{2}+4=0\;\;\lambda_{1,2}=\pm 2i. Thus we have purely imaginary eigen values. It turns out that we can draw a phase portrait of this system using the following trick. Let us multiply the first equation in (4.29) by , the second equation by and let us add them together. We get
[TABLE]
Now note that
[TABLE]
(just check this by applying the chain rule for ), and similarly
[TABLE]
hence eq.(4.30) can be rewritten as:
[TABLE]
We know that if the derivative of the function is zero () then the function is a constant, thus the above equation implies:
[TABLE]
Expression (4.31) gives a so-called first integral of our system: combination of variables which are preserved in time. It is not equivalent to the solution of our system, but using it we can draw the phase portrait of system (4.29).
Because in (4.31) is an arbitrary positive number, let us denote it as , where is just another arbitrary constant. Thus we will get the equation , which represents a graph of a circle with radius and with the center at the origin (fig.4.8a).
To find the dynamics on this circle let us find at a point (the bold point in fig.4.8a). From eq.(4.29 , i.e. the coordinate decreases at this point in the course of time. Thus we conclude that motions of trajectory along the circle will be in the clockwise direction. If we draw a series of such circles for different we will get the phase portrait as shown in fig.4.8b. The dynamics of the and variables can be found by considering motion along trajectories of the system: motion here is a rotation of a point along the circle. During this rotation the variable changes periodically (between the values and ), and thus we will see periodic oscillations. The same is valid for the variable . An example of such dynamics is represented in fig.4.8c. The equilibrium point in such system is called a center. The phase portrait in Fig.4.8 is a set of circles, which is a consequence of the symmetry of system (4.29). In a general case, if , we will also get a center point, similar to that in fig.4.8, but instead of circles we can get a series of embedded ellipses. The dynamics of the variables will always be oscillations.
Conclusion 5
If the eigen values of system (4.1) are , we have an equilibrium point called a center. The dynamics of variables are oscillations. The phase portrait is a set of embedded ellipses.
A computer generated phase portrait of the system is shown in fig.4.10a. The eigen values in this case are . The time-plot for the and variables is shown in fig.4.10(left).
Conclusion 6
The imaginary part of the eigen values results in the rotation of trajectories on a phase portrait.
Now let us consider the next two cases: .
As we discussed in section 4.3, determines convergence or divergence of the trajectories to the equilibrium. In the case we expect that the real part of the eigenvalue will give a behavior similar to the case of both negative real eigenvalues, hence we expect the convergence of trajectories to the equilibrium, as for a stable node. In addition to this, as we saw in the previous section, the imaginary part of the eigen values causes the rotation of the trajectory. If we add these two processes together we will get convergence to the equilibrium with rotation, hence trajectories will have the form of spirals. A computer generated phase portrait for this case is shown in fig.4.10b. The system is . The eigen values in this case are . The time-plot for the and variables is shown in fig.4.10(middle). The dynamics of the system are oscillations with gradually decreasing amplitude.
Conclusion 7
If the eigen values of system (4.1) are , we have an equilibrium point called a stable spiral, fig.4.10b.
The last case occurs if . This case is similar to the previous one. The only difference is that because , the real part gives motion equivalent to a non-stable node, or divergence of trajectories from the equilibrium. So, together with rotation from the imaginary part of we get the following phase portrait (fig.4.10c). This phase portrait is generated by computer for the system . The eigen values in this case are . The time-plot for the and variables is shown in fig.4.10(right). The dynamics of the system are oscillations with gradually increasing amplitude.
Conclusion 8
If the eigen values of system (4.1) are , we have an equilibrium point called a non-stable spiral, fig.4.10c.
We have found all possible types of equilibria which can occur in 2D systems: saddle, non-stable node, stable node, center, non-stable spiral and stable spiral. The next question which we will discuss here is the stability of these equilibria.
4.5 Stability of equilibrium
We will call an equilibrium point stable, if there is a neighborhood of this equilibrium, such that all trajectories which start in this neighborhood will converge to the equilibrium (Fig.4.11a). We will call the equilibrium point non-stable, if there is at least one diverging trajectory from the close neighborhood of this equilibrium (fig.4.11b).
If we analyze the stability of the 6 types of equilibria studied in the previous section we find the following:
Stable equilibria
Stable node real 2. 2.
Stable spiral
Non-stable equilibria
Non-stable node real 2. 2.
Non-stable spiral 3. 3.
Saddle point ; or real
In case of spirals and nodes the stability and non-stability is obvious. In case of a saddle point we have a converging trajectory, however the existence of the diverging trajectories (fig.4.4b) implies that this equilibrium point is non-stable. The last case, (center point), is non conclusive. Trajectories do not converge and do not diverge from the equilibrium. Usually this case is treated as neutrally stable. All these cases can be formulated in the following theorem.
Theorem 1
If all eigenvalues of the linear system (4.1) have negative real parts, then the equilibrium point is stable.
It is easy to see, that this theorem includes all listed stable equilibria. It obviously works for a stable spiral, but it also works for a stable node, because any real number can be considered as a complex number with imaginary part equal to zero. For example: can be represented as , and .
4.6 Exercises
Exercises for section 4.2
Find the general solution of the following systems of ordinary differential equations:
- (a)
\left(\begin{array}[]{c}{dx\over dt}\\ {dy\over dt}\end{array}\right)=\left(\begin{array}[]{lr}-2&1\\ 1&-2\end{array}\right)\left(\begin{array}[]{c}x\\ y\end{array}\right) 2. (b)
\left(\begin{array}[]{c}{dx\over dt}\\ {dy\over dt}\end{array}\right)=\left(\begin{array}[]{lr}3&-1\\ -2&4\end{array}\right)\left(\begin{array}[]{c}x\\ y\end{array}\right) 2. 2.
Find the solution for the following initial value problem:
[TABLE]
(Hint: See an example of solution in section 11.0.1). 3. 3.
Two different concentrations of a solution are separated by a membrane through which the solute can diffuse. The rate at which the solute diffuses is proportional to the difference in concentrations between two solutions. The differential equations governing the process are:
[TABLE]
where and are the two concentrations, and are the volumes of the respective compartments, and is a constant of proportionality. If , , and and if initially and , find and as functions of time.
Exercises for sections 4.3 and 4.4
Find eigen values and eigen vectors (for real eigen values only) of the following systems. Determine equilibrium type and sketch phase portraits. For real eigen values show non-stable, stable manifolds and several trajectories between the manifolds.
- (a)
\left\{\begin{array}[]{l}{dx\over dt}=x+4y\\ {dy\over dt}=2x+3y\end{array}\right. 2. (b)
\left\{\begin{array}[]{l}{dx\over dt}=5x-y\\ {dy\over dt}=3x+y\end{array}\right. 3. (c)
\left\{\begin{array}[]{l}{dx\over dt}=3x-5y\\ {dy\over dt}=x-y\end{array}\right. 4. (d)
\left\{\begin{array}[]{l}{dx\over dt}=-2x+y\\ {dy\over dt}=x-2y\end{array}\right. 5. (e)
\left\{\begin{array}[]{l}{dx\over dt}=-2y\\ {dy\over dt}=x-2y\end{array}\right. 6. (f)
\left\{\begin{array}[]{l}{dx\over dt}=-x-y\\ {dy\over dt}=2x+y\end{array}\right. 7. (g)
\left\{\begin{array}[]{l}{dx\over dt}=-2x-y\\ {dy\over dt}=3x+2y\end{array}\right. 5. 5.
Study the following linear system with a parameter :
[TABLE]
Find the types of equilibrium which are possible for different values of . Give the parameter region for each equilibrium and draw qualitative phase portraits. For which parameter values is the equilibrium stable? 6. 6.
For which values of parameters and does the following system has periodic oscillations (i.e. a center equilibrium point):
[TABLE] 7. 7.
Compartmental models play an important role in different parts of population biology, pharmacology and biochemistry. They describe the interaction between several processes, which may be interactions of populations, chemical reactions, etc. A two compartment model is schematically shown in fig.4.12. It represents two interacting species and . The concentration of the species can be changed either due to a transition to species with the rate given by , or can die with the rate . Similar transitions exist also for the species . The rates of these processes are specified in the figure.
- (a)
Derive a system of two differential equations for species and . 2. (b)
If find equilibrium type. Is it stable? 3. (c)
Determine the stability of the equilibrium in a general case when
4.7 Additional concepts (appendix)
4.7.1 General solution for complex eigen values
Notes about general solution of a system with complex eigenvalues
It turns out that if formula (4.6) is still valid, so the solution can be represented in the form:
[TABLE]
but as we get
[TABLE]
This expression gives the correct solutions of (4.5). However, the form of the solution is not good. First, because are complex, the eigenvectors \left(\begin{array}[]{c}v_{1x}\\ v_{1y}\end{array}\right) and \left(\begin{array}[]{c}v_{2x}\\ v_{2y}\end{array}\right) will also be complex. Next, we should consider as general complex constants. Therefore (4.33) is a quite complicated expression which gives and as complex valued functions of time . Mathematically it is correct and if we substitute (4.33) into the original equation (4.5) we get an equality. However, we need to draw a phase portrait of (4.5) on the -plane, i.e. we only need the real solutions of our system. Such real solutions are present in the general expression (4.33), i.e. they are a subset of all the possible solutions. However, extracting them from (4.33) is not a simple task. We will highlight the main idea behind this derivation below. To derive the general expressions we will need the following Euler formula.
Euler formula
The Euler formula gives a representation of in terms of trigonometric functions. It is quite unexpected:
[TABLE]
or in another representation:
[TABLE]
When you see this formula for the first time it looks quite crazy. We know that and come from the simple geometry of triangles, and is a special exponential function. Why are these functions connected together in such a simple way (4.35)?
To prove this formula, one should use Taylor series. However, here I will present another simpler derivation of (4.35) on the basis of differential equations.
At the beginning of this chapter, in order to find a solution of (4.5), we first considered a one dimensional differential equation , and we found its solution . Consider the following initial value problem for this equation:
[TABLE]
This initial value problem has the unique solution . So we say, that the solution of (4.36) is . But we can also say it vice versa: we can define the exponential function as the function which satisfies the initial value problem (4.36). For example, if we give to a person just this equation and a computer, he will be able to solve it and to draw the graph of , even without knowledge about exponential functions. The advantage of such a definition is that it can be easily extended to complex numbers. So, let us define as a function which satisfies the initial value problem (4.36) with
[TABLE]
In other words: must be the function , such that , and the derivative of this function must be equal to this function times , i.e. . Let us find an expression which satisfies these conditions. It turns out that it will be . Let us check it. The first condition is satisfied:
[TABLE]
To check the second condition we write:
[TABLE]
if we replace by we get:
[TABLE]
i.e. the second condition is also satisfied. So the function gives the solution of (4.37), hence it is the same as or and we get the Euler formula (4.35). The formula (4.34) is just the formula (4.35) in which instead of the letters are used. To find we write:
[TABLE]
General solution
Now let us find a solution of a system with imaginary eigenvalues. As we know the general solution is given by the formula (4.33) and because of the Euler formula we can rewrite it in the following way:
[TABLE]
where are arbitrary complex constants and \left(\begin{array}[]{c}v_{1x}\\ v_{1y}\end{array}\right) and \left(\begin{array}[]{c}v_{2x}\\ v_{2y}\end{array}\right) are complex eigen vectors. Now we should find a real part of this complicated expression and get a general real solution of our system in this case. A general solution of a system of two differential equations should depend on two arbitrary constants as the initial value of each of the variables can be arbitrary. We will use this fact to find the general solution. Our idea is instead of extracting all real solutions from (4.39) we will find just two real solutions. By multiplying them by two arbitrary constants we will get a general solution.
We will be able to find these two solutions from the first term of (4.39):
[TABLE]
Let us extract real and an imaginary parts of this term. If we use the formula (2.36) for the eigen value we find the eigen vector : It has the real and imaginary parts:
[TABLE]
The vector has the real part and imaginary part . So, the term (4.40) can be written as:
[TABLE]
If we denote:
[TABLE]
the term (4.40) can be rewritten as
[TABLE]
Let us prove that both and are the real solutions of (4.5). For that we will use the fact that (4.39) gives a solution of (4.5) and hence is a complex solution of (4.5) as it is a part of (4.39).
System (4.5) in a matrix form can be written as:
[TABLE]
As is a solution, it satisfies (4.43):
[TABLE]
Equating the real and imaginary parts yields
[TABLE]
Hence and are real solutions of (4.5). Finally, because and are real solutions of (4.5) the general solution is given by the formula:
[TABLE]
where and are arbitrary constants and and are given by (4.42).
Example Find the general solution of the following system.
[TABLE]
Solution. The eigen values are given by: Det\left|\begin{array}[]{lr}-\lambda&2\\ -2&-\lambda\end{array}\right|={\lambda}^{2}+4=0,\;\lambda_{1,2}=\pm\sqrt{-4}=\pm 2i. So the eigen vector corresponding to the eigen value is:
[TABLE]
So,
[TABLE]
Therefore the solution from the formulas (4.44) is:
[TABLE]
Or if we denote and we get:
[TABLE]
Let us check that (4.50) does gives a solution of (4.45). Substitution of (4.50) into equation (4.45) yields:
[TABLE]
So, (4.50) is a solution of (4.45).
Finally note that formula (4.31) which we got for the same equation (4.45) is of course valid for functions given in (4.50). For that you need to find with and given by (4.50). A direct computation will give us
[TABLE]
i.e. the same result as in (4.33). (Note that in order to get the final result you need to apply several times the well-known formula from trigonometry ).
Chapter 5 System of two non-linear differential equations
5.1 Introduction and first definitions
5.1.1 Phase portrait
After analyzing linear systems let us consider a general non-linear system which can be written in the following general form:
[TABLE]
Many biological systems are described by such systems. One of the classical examples of ecological models (the predator-prey model) can be derived as follows. Let us consider the prey population with a logistic growth given by eq.(3.21): , which interacts with the predator and let us assume that the effect of the predator on the prey population is given by the term . Then, if we assume that the growth of the predator population is proportional to the predator prey interaction and that the death rate of the predator is given by , we will get the following system of differential equations:
[TABLE]
Formally system (5.2) describes the predator-prey interactions with competition in the prey population. It has several parameters, which account for the specific properties of the populations. Let us study it for :
[TABLE]
If, as in section 4.1, we solve this system on a computer, we will get the following phase portrait (fig.5.1b). We see, that in the course of time all trajectories approach the stationary values of .
The main aim of our course is to develop the procedure of drawing a phase portrait of a general system of two non-linear differential equations without using a computer. We expect that as for 1D differential equation (section 3.2.1) and for 2D linear systems (section 4.1) the phase portrait should include two main elements: equilibria points and flows (trajectories) between them. Let us define first equilibria points of a general non-linear system (5.1).
5.1.2 Equilibria
In the 1D case and for 2D linear systems the equilibria were the points where our system is stationary: placed at equilibrium point the system will stay there forever. Therefore, for 1D equation equilibria were determined as the points where , i.e. where . For 2D linear system (section 4.1) we required that both variables and are stationary at equilibria points, i.e. both and . For a general non-linear system (5.1) these conditions yield the following definition of equilibria:
Definition 7
A point is called an equilibrium point of a system (5.1) if
[TABLE]
Equilibria in two dimensions are also stationary points, i.e. if system is placed at the equilibrium it will stay there forever. Thus this trajectory will contain just one point.
Example. Find the equilibria of the system (5.3):
Solution To find the equilibria we need to solve a system of algebraic equations (5.4) which in our case becomes:
[TABLE]
From the second equation we find , which can be either when or when . Substitution of to the first equation yields . This equation has two solutions and . Substitution of the other case to the first equation gives , or . Thus we have found three equilibria points: and .
We see in fig.5.1 that point is indeed an important attractor of our system which determines the final state of the populations. The other two points are not apparent in fig.5.1, however, as we will see later they also account for important changes of trajectories of our system.
Thus we have defined equilibria for 2D systems. Our next step is to understand what is the 2D analog of flows, which on the 1D phase portrait were represented by the ’’ or arrows.
5.1.3 Vector field
In 1D flows, visualizations of the direction of change of the variable were given via the sign of its derivative . In 2D, both variables can change and the rate of their change is given by the derivatives and . In 1D we were able to find the direction of flow at any point from the right hand side function of the equation . Similarly in 2D we can find and at any point from the right hand sides of system (5.1) ( functions and ). For example for system (5.3) at a point we find , and However, what do these two numbers show? They tell us that if the size of the prey population and the size of the predator population is , then the prey population decreases with the rate of and the predator population grows with the rate of . On the phase plane this will result in a shift of a point representing populations from point (1,1) (point in fig.5.2a) to some point which is to the left and upward from point . Let us make it more quantitative. We know that the rate of change of in our case is times larger than the rate of change of . This determines the direction of shift of point relative to point . The easiest way to represent it is to draw from point a horizontal arrow heading to the left with the length of and a vertical arrow heading upward with the length of . The direction of the overall shift will be given by the resultant vector of these two vectors fig.5.2b. The resultant vector will give us the direction tangent to the trajectory which goes through the given point. We can generalize this result as:
Conclusion 9
At any point of a phase space for system (5.1), we can define the vector with the components . Such vectors will be tangent to the trajectories of our system. We can find this vector field without a solution of our system, just from the right hand sides of our system.
Note, that the length of the vector in fig.5.2 is not important, as we are interested in the direction, only. If we apply the same procedure at many points and represent the directions by shorter vectors we will get the following vector field of the system (5.3) (fig.5.3).
We see that although the vector field describes qualitatively the direction of trajectories in the phase portrait of our system, it does not give us information on convergence/divergence of trajectories, thus we cannot determine the attractors of our system, which is crucial for our study. However, as we will see in chapter 6 we can elaborate the vector field based methods and in many cases will be able to obtain convergence/divergence information using the so-called graphical Jacobian approach.
But in order to derive this approach we need to understand how to find attractors of a general non-linear system (5.1) using an analytical approach. For that we need to establish the relation of the general non-linear system (5.1) and the general linear system (4.1) that we studied in chapter 4.
5.2 Linearization of a system: Jacobian
Consider a general system of two differential equations:
[TABLE]
In this section we will show that close to equilibrium point the phase portrait of this general non-linear system (5.6) can be found from the solution of the linear system (4.1) that we studied in chapter 4. As a consequence we will get six possible types of equilibria, that are: saddle, non-stable and stable node, non-stable and stable spiral, and a center.
Our main tool here will be a formula (5.7) for the approximation of function of two variables around point which we derived in section 2.4:
[TABLE]
where and are the values of the partial derivatives at the point , i.e. they are just numbers.
We will apply this formula to approximate the functions and of our system (5.6) and later solve the approximated system and find the phase portraits close to the equilibrium.
Let us start the derivation. Assume that system (5.6) has an equilibrium point at (). This means (see (5.4)) that:
[TABLE]
Let us approximate close to the equilibrium using the formula (5.7):
[TABLE]
As we assumed () is an equilibrium, i.e. and we get
[TABLE]
A similar approach for yields:
[TABLE]
If we replace the right hand sides of (5.6) by their approximations (5.9), (5.10), we get the following system:
[TABLE]
The system (5.11) is simpler than the original system (5.6), as the partial derivatives in (5.11) are constants (numbers, as they are evaluated at the equilibrium point ). So we can rewrite (5.11) as :
[TABLE]
where . We can simplify (5.12) even more. For that let us introduce new variables:
[TABLE]
were are new unknown functions of . If we substitute them into the right hand side of (5.12),we get:
[TABLE]
In order to substitute and into the left hand side of (5.14), note that , i.e. (here because is a constant). Similarly, . After replacing by and by in (5.14) we get:
[TABLE]
System (5.15) coincides with a general linear system (4.1) that we studied in chapter 4 and for which we found six possible types of solutions resulting in six possible phase portraits: saddle, non-stable and stable node, non-stable and stable spiral, and a center. However, how can we use the results of study of system (5.15) for the study of the original system (5.6)? In order to derive (5.15) we made two steps: (1) we used formula (5.7) for function approximation; (2) we changed variables to . Let us analyze each of these two steps. (1) As we discussed earlier, formula (5.7) gives a good approximation of the function only if is close to the point of approximation . So we can use the linear system (5.15) for approximating the solutions of the non-linear system (5.6) only close to the equilibrium point . (2) about change of variables. Equation (5.13) gives the variables via the variables . However, we can also solve the equations and find how will be expressed in terms of :
[TABLE]
Using this expression we can draw the trajectories for our original variables if we know the trajectories of the linearized system (5.15). Indeed, as the coordinate of the trajectory equals plus number , and the coordinate equals plus number , the only thing what we need to do is just to draw the trajectory and shift its -coordinate by and the -coordinate by units. As we know in a general linear system (5.15) the phase portrait is centered around a point , therefore all what we would need to do in order to draw the phase portrait of the non-linear system (5.6) close to its equilibrium is just to shift the phase portrait of the linear system (5.15) to a location of equilibrium in non-linear system (5.6). If, for example, the linearized system (5.15) will have a stable node type phase portrait, then a non-linear system (5.6) will have the same stable node point but around an equilibrium as shown in fig.5.4.
Conclusion 10
The phase portrait of liner system (5.15) close to the origin () is similar to the phase portrait of non-linear system (5.6) close the to equilibrium point . To draw the phase portrait of non-linear system (5.6) close to equilibrium, we need to shift a phase portrait of linear system (5.15) from the origin the equilibrium point .
To find the linearized system (5.15) we need to find the equilibrium point () and compute the following four numbers: the values of derivatives of right hand sides of our system at this equilibrium:
[TABLE]
So, system (5.15) can be written as:
[TABLE]
From coefficients of this system we can construct a matrix that is called the Jacobian
[TABLE]
Example Find qualitative phase portrait of system (5.3) close to a nontrivial equilibrium point (i.e. at the equilibrium where ).
[TABLE]
Solution In section 5.1.2 we found that this system has three equilibria , , and . The nontrivial equilibrium is . To find the Jacobian of our system we compute the partial derivatives (5.18) and evaluate them at the equilibrium point. In our case .
at point this derivative equals . Similarly: , hence the linearization of our system at point is
[TABLE]
or the Jacobian is:
[TABLE]
In order to find a phase portrait of this linear system we need to find eigen values of the Jacobian matrix from the following characteristic equation (2.31):
[TABLE]
Using ’abc’ formula we find that:
[TABLE]
or, and . Because both eigen values are real and negative we will have a stable node and a phase portrait qualitatively similar to that in fig.5.4b.
5.3 Determinant-trace method for finding the type of equilibrium
In this section we derive a simple method for finding signs of eigen values and the type of equilibrium of the linear system (4.1). The results will also be applicable for a general non-linear system (5.6), because, as discussed in the previous section, the equilibrium type of a non-linear system can be found from its linearization.
The eigen values of system (4.1) are the roots of the characteristic equation (2.31):
[TABLE]
or:
Det\left|\begin{array}[]{lr}a-\lambda&b\\ c&d-\lambda\end{array}\right|=(a-\lambda)(d-\lambda)-cb={\lambda}^{2}-\lambda(a+d)+ad-cb=0
Let us express the last equation in a slightly different form using the following definition:
Definition 8
The trace of the matrix \left(\begin{array}[]{lr}a&b\\ c&d\end{array}\right) is .
The determinant of the matrix is
So we can rewrite the characteristic equation using the definition of the trace and the determinant as follows:
[TABLE]
We see that although the original system depends on four parameters the characteristic equation depends only on two parameters and , thus if we know the determinant and the trace of our system we can find the eigen values and the type of the equilibrium of the system. Indeed, from (5.22) we can easily find that the eigen values are:
[TABLE]
Roots of the equation (5.23) are as the roots of any quadratic equation, connected in the following way to the coefficients of the equation:
[TABLE]
[TABLE]
To prove the properties (5.24) and (5.25), just note that if and are the roots of a quadratic (5.22), then it can be written as: . The direct computation yields:
[TABLE]
If we now compare the left and the right hand sides of the last equation we will get both properties (5.24) and (5.25).
Let us start classification.
If , then , so we have real roots. From (5.25) we conclude that their product is negative, i.e. roots have different signs, i.e. , or and we have a saddle point. 2. 2.
If , then can be negative as well as positive. This means the roots can be real, or complex. Let us consider the case of real roots first, i.e.
[TABLE]
If (5.26) holds, the roots are real. Next, let us use the property . In the case of , the product of the roots is positive, i.e. the roots have the same sign. They can be both positive, or both negative. The sign of the roots can be found from the trace of the matrix (). When , then , and and we have a non-stable node. When , and and we have a stable node. Let us formulate it as a separate case: 3. 3.
If , and the equilibrium is a stable node.
Let us put this information into a graph (fig.5.5). On this graph let us use as the -axis and as the -axis. Case 1 of a saddle point then corresponds to the lower half plane (region 1). The line (5.26), which separates real and complex roots, is the parabola given by , or . Real roots are below this line. Therefore, in region 2, where , we have case 2 of a non-stable node. In region 3 and we have case 3 of a stable node. 4. 4.
If , and , we have complex roots. In accordance with (5.24) and (4.27) they are:
[TABLE]
Hence, . From this we immediately see that if , we have a non-stable spiral (region 4). 5. 5.
If , and , we have a stable spiral (region 5). 6. 6.
The last case of a center point appears when , or when and . In our graph is it the upper part of the axis (region 6)
Let us apply this method for several examples:
Example. Find the equilibrium type for the following systems:
a)\left(\begin{array}[]{c}{dx\over dt}\\ {dy\over dt}\end{array}\right)=\;\left(\begin{array}[]{lr}1&2\\ 3&4\end{array}\right)\left(\begin{array}[]{c}x\\ y\end{array}\right);\quad b)\left(\begin{array}[]{c}{dx\over dt}\\ {dy\over dt}\end{array}\right)=\;\left(\begin{array}[]{lr}4&1\\ 1&2\end{array}\right)\left(\begin{array}[]{c}x\\ y\end{array}\right);\quad
c)\left(\begin{array}[]{c}{dx\over dt}\\ {dy\over dt}\end{array}\right)=\left(\begin{array}[]{lr}-2&3\\ -2&1\end{array}\right)\left(\begin{array}[]{c}x\\ y\end{array}\right);\quad d)\left(\begin{array}[]{c}{dx\over dt}\\ {dy\over dt}\end{array}\right)=\left(\begin{array}[]{lr}1&-1\\ 2&-1\end{array}\right)\left(\begin{array}[]{c}x\\ y\end{array}\right);\quad
Solution Our plan is to find for the corresponding matrices and make a conclusion.
a) The matrix is A=\left(\begin{array}[]{lr}1&2\\ 3&4\end{array}\right);\;trA=1+4=5;\;detA=1*4-2*3=-2. Therefore we have case 1, hence the system has a saddle point.
b) The corresponding matrix is A=\left(\begin{array}[]{lr}4&1\\ 1&2\end{array}\right);\;trA=4+2=6;\;detA=4*2-1*1=7. As we need to check that the roots are real. We have , hence the roots are real and we have case 2, hence the system has a non-stable node.
c) , we have complex roots and we are in region 5 and have a stable spiral
d) , so we have case 6 and we have a center point, or oscillation in our system.
5.4 Exercises
Exercises sec.5.1
(A) Find equilibria of the following systems \left\{\begin{array}[]{l}{dx\over dt}=f(x,y)\\ {dy\over dt}=g(x,y)\end{array}\right. (see definition in section 5.1.2).
(B) Find the following partial derivatives at each equilibrium point ().
(a) \left\{\begin{array}[]{l}{dx\over dt}=-4y\\ {dy\over dt}=4x-x^{2}-0.5y\end{array}\right. (b) \left\{\begin{array}[]{l}{dx\over dt}=9x+y^{2}\\ {dy\over dt}=x-y\end{array}\right.
(c) \left\{\begin{array}[]{l}{dx\over dt}=2x-xy\\ {dy\over dt}=-y+y^{2}x\end{array}\right. (d) \left\{\begin{array}[]{l}{dx\over dt}=(1-x-3y)x\\ {dy\over dt}=(1-2x-2y)y\end{array}\right.
(Hint: See an example of solution in section 11.0.2). 2. 2.
Find equilibria of the following Lotka Volterra model with competition in the prey population. Determine for which parameter values all equilibria are non-negative.
[TABLE] 3. 3.
Protein synthesis depends on DNA transcription () making mRNA molecules () and translation () of mRNA into proteins (). Some proteins inhibit the transcription of their own mRNA (). mRNA and proteins are degraded at rates b and d, respectively. This process gives the following of two differential equations. Find equilibria of this model.
[TABLE] 4. 4.
Mathematical epidemiology also makes use of simple ODE models. One of these models describes the number of susceptible individuals () and infected individuals (). Individuals are born at rate (), and die at rate (). Susceptible individuals can become infected when they come into contact with infected individuals (). Once infected, an individual has a certain death rate (); this may be different from the death rate of non-infected individuals. This process can therefore be modeled by the following equations: \left\{\begin{array}[]{l}dS/dt=B-\beta SI-\mu S\\ dI/dt=\beta SI-\alpha I\;\;\;\;\;\;\;B,\alpha,\beta,\mu>0\end{array}\right.
Find equilibria of this model. Determine for which parameter values all equilibria are non-negative.
Exercises sec. 5.2 and 5.3
Find the type of equilibria using the - method. Determine the stability of the equilibrium.
(a) \left\{\begin{array}[]{l}{dx\over dt}=3x+y\\ {dy\over dt}=-20x+6y\end{array}\right. (b) \left\{\begin{array}[]{l}{dx\over dt}=2x+y\\ {dy\over dt}=2x-10y\end{array}\right. (c) \left\{\begin{array}[]{l}{dx\over dt}=2x+y\\ {dy\over dt}=5x-2y\end{array}\right. (d) \left\{\begin{array}[]{l}{dx\over dt}=x+10y\\ {dy\over dt}=-10x-y\end{array}\right. 6. 6.
Consider the following model for the algae population: \left\{\begin{array}[]{l}{dx\over dt}=2x(1-y)\;\;x\geq 0;\\ {dy\over dt}=2-y-x^{2}\;\;y\geq 0.\end{array}\right.
- (a)
Find equilibria 2. (b)
Find the general expression for the Jacobian of this system 3. (c)
Determine type of each equilibrium using - method 4. (d)
Draw qualitative local phase portraits around each equilibrium point 7. 7.
Consider the following biological model:
[TABLE]
- •
For which values of parameter the system has only one equilibrium?
- •
Determine stability and type of this equilibrium in found parameter range (in which system (5.27) has only one equilibrium). 8. 8.
Study the Lotka-Volterra model for a predator-prey system: \left\{\begin{array}[]{l}{dN\over dt}=aN-bNP\\ {dP\over dt}=cNP-dP\end{array}\right.
here denotes the prey population, denotes the predator population and are parameters
- (a)
Find the nontrivial equilibrium of the system (i.e. an equilibrium where ). 2. (b)
Find the linearization of the system at this point (i.e. the Jacobian matrix) 3. (c)
Determine the type of the equilibrium 4. (d)
Sketch the phase portrait around this equilibrium. Which kind of dynamics do we expect here?
Additional exercises
Consider the system:
- (a)
Check that () is an equilibrium point of the system 2. (b)
Find the general expression for the Jacobian of this system 3. (c)
Find the Jacobian at the point () 4. (d)
Write the linearization of the system close to the equilibrium () 10. 10.
Find the equilibria of the following systems. Compute the Jacobian at the equilibria points. Determine type of each equilibrium using - method. Draw qualitative local phase portraits around each equilibrium point.
(a)\left\{\begin{array}[]{l}{dx\over dt}=y^{2}-3x+2\\ {dy\over dt}=x^{2}-y^{2}\end{array}\right. (b) \left\{\begin{array}[]{l}{dx\over dt}=y\\ {dy\over dt}=-x+x^{3}\end{array}\right. 11. 11.
Revisit solutions of problems 1. Use found equilibria and partial derivatives at these equilibria points to determine type of each equilibrium using - method. Draw qualitative local phase portraits around each equilibrium point. 12. 12.
Consider a modification of the Lotka-Volterra model, which includes competition in the prey population (): \left\{\begin{array}[]{l}{dN\over dt}=aN-eN^{2}-bNP\\ {dP\over dt}=cNP-dP\end{array}\right., where the parameters and the variables , .
- (a)
Find all equilibria of this system. 2. (b)
Compute the Jacobian at each equilibrium point. 3. (c)
At which parameter values do we have a non-trivial equilibrium (i.e. an equilibrium at which both and are positive). Find stability of this equilibrium. 13. 13.
Consider the following model for cardiac tissue:
[TABLE]
Here the variable accounts for the transmembrane potential, the variable accounts for the refractory period and are the parameters.
The shape of the action potential in cardiac tissue is an important characteristic of myocardium. If the recovery of the transmembrane potential shows oscillation as in fig.5.6b. it can cause dangerous cardiac arrhythmias. From a mathematical point of view the oscillations in fig5.6b occur when system (5.28) has an equilibrium point which is a stable spiral. Monotonous recovery occurs when this equilibrium is a stable node.
Determine for which parameter values we will have situation fig.5.6a and for which parameter values we will have situation fig.5.6b.
Chapter 6 Graphical methods to study systems of differential equations
In the first section of this chapter we will start from the vector field of a general non-linear system introduced in section 5.1.3, and find how we can approximate the vector field by the so-called null-cline method, without using a computer. Then, in section 6.2, we will show that the null-clines can be used not only for vector field approximation but also for determining the type of an equilibrium point without explicit computation of the Jacobian of the system.
6.1 Null-clines
We introduced the vector field in section 5.1.3, where we showed that for a general 2D system
[TABLE]
the direction of the trajectories on the phase portrait will be along the vector with the components . Thus, to draw the vector field of a particular system we need to evaluate the values of the functions in many points which usually requires using a computer. In this section we will develop a so-called method of null-clines, which will allow us to sketch a qualitative picture of the vector field analytically. The main idea here is similar to what we did in the 1D case, where we have represented the derivative by arrows of two types: ’’ for and ’’ for . In 2D the vector field has two components . Each of these components can be positive, or negative. Therefore, we can have the following four cases shown fig.6.1. (Note, that we use different line types there that will be important in the future).
The method of null-clines represents the vector field using these four cases. The main idea behind this method is the following. If we compare cases I and II we see that they differ by the sign of the derivative: for case I and for case II . Therefore these cases are separated by the boundary where (fig.6.2a). We know that for system (6.1) , therefore the boundary between cases I and II is given by the condition
[TABLE]
Geometrically equation (6.2) gives a graph of one or more lines in the -plane (see section 1.4). Note that at this line the horizontal component of the vector field is zero, therefore the vector is vertical.
Similarly, the transition from case I to case III occurs when (fig.6.2b). As we know for the system (6.1) , hence the separation line in this case is given by the equation:
[TABLE]
and the direction of the vector field at this line is horizontal.
Equation (6.3) will give us not only the boundary between cases I and III but also a boundary between cases II and IV, as the transition between cases II and IV also occurs when . In general equations (6.2),(6.3) give two (or more) lines on the -plane, which separate the plane into several regions with different directions of vectors (cases I-IV).
These lines are called null-clines.
Definition 9
The -null-cline (or null-cline) is the set of points satisfying the condition . The -null-cline (or null-cline) is the set of points satisfying the condition .
To use the method of null-clines it is useful to note that at the -null-cline the -component of the vector changes its sign, and at the -null-cline the -component of the vector changes its sign. To use this rule effectively we will always denote the vector field components using lines of different types: the horizontal component as a dashed line and the vertical component as a solid line. Let us use these ideas and formulate a plan for finding the vector field using null-clines.
Plan of null-cline analysis for system (6.1)
We assume that on the -plane the -axis is the horizontal axis and the -axis is the vertical axis.
Draw null-clines from the equation using a dashed line and null-clines from the equation using a solid line. 2. 2.
Choose a point in one of the regions in the plane and find the and the -components of the vector field. Use the dashed line for the component and the solid line for the component. For finding the directions use the following rule: if the component is directed as ’’, if it is directed as ’’; if the -component is directed as , it is directed as . 3. 3.
Find the vector field in the adjacent regions using the following rule:
- (a)
change the direction of the dashed component of the vector field if in order to get to the adjacent region you cross the dashed null-cline 2. (b)
change the direction of the solid component of the vector field if in order to get to the adjacent region you cross the solid null-cline 3. (c)
show the direction of the vector field on the null-clines.
Note, that instead of dashed and solid lines you can use lines of different colors. Then the last step of this plan would be: change the direction of the component of the same color as the color of the null-cline which we cross to get to the adjacent region.
Note, that although this plan works good in most of cases, there are some situations when components of the vector field do not change their sign at the corresponding null-cline. These are special so-called degenerate cases (exceptions). We will not consider them in our course.
Example Find the vector field of the following system using null-clines.
[TABLE]
Solution. We follow our plan as follows
The null-clines are given by the equation , i.e. . This equation has two solutions: (the vertical line which coincides with the -axis) and which is a straight line with the negative slope which goes through the point . The graphs are shown using dashed lines in fig.6.3. The null-clines are given by the equation , i.e. , which also has two solutions: (horizontal line which coincides with the -axis) and (vertical line through the point ). Graphs are shown by solid lines in fig.6.3. 2. 2.
We find that at point and , hence the direction of the dashed arrow is to the left ’’ and of the solid arrow is upward . 3. 3.
Now we complete the picture. For example, to get into the region to the left from point () we cross the solid line, thus we change the direction of the solid component here. Similarly for the other regions. We get the picture as in fig.6.3. Finally we show the vector field on the null-clines.
We see that the vector field in fig.6.3 is a good approximation of the flow in fig.5.3. We also see that attractor () is a special point in fig.6.3: the point of intersection of the and null-clines. This is not a coincidence. If we compare the conditions for finding the equilibrium point (5.4) and equations for null clines (6.2), (6.3), we see that the first equation for finding equilibria is also the equation for null-cline and the second equation for finding the equilibrium point , is the equation for the -null-cline. Thus, the solution of system (5.4), which gives the points satisfying both equations, gives the points which belong to both null-clines, i.e. the points of intersection of the null-clines. So we found that:
Conclusion 11
Equilibria are the points of intersection of the and -null-clines.
Note that this definition applies points of intersection of different null-clines only. For example intersection of null clines of the same type in fig.6.3 at points (0.5,0) and (0,2) do not give equilibria of system (5.3), while intersection of the different null-clines at points (0,0),(1,0), (0.5,1) give the equilibria of this system.
Our next step will be to study how can we apply the null-clines for finding the types of equilibrium.
6.2 Graphical Jacobian
The method which we present here allows, in a number of cases, to find the type of an equilibrium from the null-clines of the system, i.e., even without computation of partial derivatives of the Jacobian in the equilibrium. The main idea behind this method can be seen from the scheme in Fig.6.4. Let us consider an equilibrium point and two close points: one located to the right, with the coordinates , the other upward with the coordinates . Because we assume that is an equilibrium point, and (see (5.4)). We can approximate the partial derivative at using a formula similar to (1.9) as : , but because , we get . If we apply the same approach for all derivatives constituting the Jacobian at the equilibrium point we get:
[TABLE]
This approximation will be better if points and are closer to the equilibrium point, i.e., if is small. It turns out that in many cases the exact values of the derivative will not be important for us and we will be able to find the equilibrium type from just the sign of the components of the Jacobian. From (6.5) it is clear that the sign of the Jacobian components is the same as the sign of the functions at the appropriate points, e.g. the sign of is the same as the sign of , etc. Let us now recall, that the sign of the functions is represented on the vector field of our system. In fact, ’’ means that and occurs at the points where , the ’’ means that and occurs at the points where (see fig.6.2), etc. Thus from the vector field we can easily determine the sign of the components of the Jacobian matrix. For example, the negative direction of the -component in fig.6.4a will mean that and hence , the positive direction of the -component in fig.6.4a will mean that and thus . Similarly, in Fig.6.4b the vector field at point is vertical, thus and .
We will formulate this result as the following conclusion:
Conclusion 12
The sign of the and vector field components to the right from the equilibrium point give the sign of and components of the Jacobian. The sign of the and vector field components upward from the equilibrium point give the sign of and components of the Jacobian matrix J=\left(\begin{array}[]{lr}{\partial f\over\partial x}&{\partial f\over\partial y}\\ {\partial g\over\partial x}&{\partial g\over\partial y}\end{array}\right) at this equilibrium point (fig.6.4c).
Three notes on the application of this rule:
- •
The testing points must be exactly horizontal and *exactly vertical * with respect to the equilibrium point.
- •
Testing points should be as close as possible to the equilibrium and should never cross a null cline when going from the equilibrium. (Putting the testing point to the right from the dashed line in fig.6.4a will be wrong).
- •
If a null cline is exactly horizontal or exactly vertical then one of the corresponding derivatives will be zero. ( In fig.6.4b, ).
Examples. Let us consider several examples of the application of this graphical approach for null-clines presented in Fig.6.5 and Fig.6.6. On all these figures the testing points are marked by filled circles.
From fig.6.5a we find the following components of the Jacobian: J=\left(\begin{array}[]{lr}-\alpha&-\beta\\ +\gamma&-\delta\end{array}\right), where stand for unknown positive numbers. Thus we see that , and , thus (see fig.5.5) we have a stable equilibrium (stable node or stable spiral). On the basis of these data we cannot say whether this equilibrium is a node or a spiral as we cannot compute whether the discriminant of the Jacobian is positive or negative, as this depends on the exact values of the partial derivatives.
From fig.6.5b we find: J=\left(\begin{array}[]{lr}-\alpha&0\\ 0&-\delta\end{array}\right). Thus , and , thus we have a stable node.
From fig.6.5c we find: J=\left(\begin{array}[]{lr}+\alpha&0\\ 0&-\delta\end{array}\right). Thus , and we have a saddle.
From fig.6.6a we find: J=\left(\begin{array}[]{lr}+\alpha&+\beta\\ +\gamma&-\delta\end{array}\right), hence , and we again have a saddle.
From fig.6.6b we find: J=\left(\begin{array}[]{lr}+\alpha&+\beta\\ -\gamma&-\delta\end{array}\right). , . Because we do not know the values of the coefficients , just their signs, we do not know if the and are positive or negative, thus we cannot determine the equilibrium type in this case using the graphical Jacobian.
Finally from 6.6c we find: J=\left(\begin{array}[]{lr}+\alpha&+\beta\\ +\gamma&+\delta\end{array}\right). Thus , . Again we cannot determine the sign of the , but positive implies that the equilibrium will be unstable. We do not know its type: non-stable node, non-stable spiral and a saddle are all possible.
We see that the method of graphical Jacobian is a useful tool for finding the equilibrium type from the vector field, but sometimes it is not sufficient and we need to know not only the sign but also the values of the Jacobian coefficients.
6.3 Exercises
Sketch the phase portraits of the following systems using null-clines. Try to find the type of equilibrium using the graphical Jacobian (section 6.2). If you are unable to find the equilibrium type using the graphical Jacobian, find it using the method of sec.5.3. Sketch a qualitative phase portrait (without computation of eigen values and eigen vectors).
- (a)
\left\{\begin{array}[]{l}{dx\over dt}=3x+y\\ {dy\over dt}=-x-y\end{array}\right. 2. (b)
\left\{\begin{array}[]{l}{dx\over dt}=x+2y\\ {dy\over dt}=-2x-2y\end{array}\right. 3. (c)
\left\{\begin{array}[]{l}{dx\over dt}=y\\ {dy\over dt}=3x\end{array}\right. 4. (d)
\left\{\begin{array}[]{l}{dx\over dt}=x-4y\\ {dy\over dt}=x+y\end{array}\right. 2. 2.
The figure below shows the null-clines of three systems of two non-linear differential equations
- (a)
Complete the figures by showing the direction of the vector field in all regions on the plane and on the null-clines. 2. (b)
Mark the equilibria and find their types using the graphical Jacobian (if possible). 3. 3.
Algae model using graphical Jacobian (see also problem 6 from section 5.4):
[TABLE]
- (a)
Sketch the vector field for the Algae system using null-clines. 2. (b)
Find equilibria. 3. (c)
Find the type of each equilibrium using the graphical Jacobian and draw qualitative local phase portraits. 4. 4.
Lotka Volterra model using graphical Jacobian (see also problem 8 from section 5.4):
[TABLE]
- (a)
Sketch the vector field of the system using null-clines. 2. (b)
Find equilibria. 3. (c)
Find the type of each equilibrium using the graphical Jacobian and draw qualitative local phase portraits. 4. (d)
Does the vector field change if we change the values of parameters ?
Additional exercises
Find the vector fields of these systems using null-clines. Find equilibria. Determine the type of each equilibrium using the graphical Jacobian and draw qualitative local phase portraits.
- (a)
\left\{\begin{array}[]{l}{dx\over dt}=x-3y\\ {dy\over dt}=x+y\end{array}\right. 2. (b)
\left\{\begin{array}[]{l}{dx\over dt}=y\\ {dy\over dt}=-x\end{array}\right. 3. (c)
\left\{\begin{array}[]{l}{dx\over dt}=9x+y^{2}\\ {dy\over dt}=x-y\end{array}\right. 4. (d)
\left\{\begin{array}[]{l}{dx\over dt}=y^{2}-x^{2}\\ {dy\over dt}=y-1\end{array}\right.
Chapter 7 Plan of qualitative analysis and examples
7.1 Plan
Let us now formulate a plan to qualitatively study a system of two differential equations with two variables and consider several examples.
We study the system:
[TABLE]
Our main aim is to plot the phase portrait of this system and then predict its dynamics. Based on methods which we have developed we will do it in two steps:
- I
Null-cline analysis 2. II
Jacobian analysis
Where the Jacobian analysis can be either performed using the determinant-trace method from section 5.3, or using the graphical Jacobian from sec.6.2. Note, that the determinant-trace method always give us a definitive answer, while the graphical Jacobian method sometimes fails. In more details:
Null-cline analysis
We assume that on the -plane the -axis is the horizontal axis and the -axis is the vertical axis.
Draw the null-clines from the equation using dashed lines and the null-clines from the equation using solid lines. 2. 2.
Choose a point in one of the regions on the plane and find the vector field for the -component. Denote the component as a dashed ’’ if and as a dashed ’’ if . 3. 3.
Find the vector field for the -component at the same point. Denote the component as a solid if and as a solid if 4. 4.
Find the vector field in the adjacent regions using the following rule:
- (a)
change the direction of the dashed component of the vector field if to get to the adjacent region you cross the dashed null-cline. 2. (b)
change the direction of the solid component of the vector field if to get to the adjacent region you cross the solid null-cline. 3. (c)
show the direction of the vector field on the null-clines close to the equilibrium
Jacobian analysis using the determinant-trace method
Find equilibria from equations:
[TABLE] 2. 2.
For each equilibrium , find the Jacobian at that equilibrium point
[TABLE]
Note: Do not forget to substitute into the Jacobian. 3. 3.
Determine the type of each equilibrium . For this compute
[TABLE]
[TABLE]
[TABLE]
To find the type of equilibrium use fig.5.5, or the following list:
- (a)
If ; the point is a saddle point 2. (b)
If ; the point is a non-stable node 3. (c)
If ; the point is a stable node 4. (d)
If ; the point is a non-stable spiral 5. (e)
If ; the point is a stable spiral 6. (f)
If ; the point is a center 4. 4.
Draw local phase portraits using both knowledge on the equilibrium type and the vector fields obtained using null-cline analysis 5. 5.
Connect local phase portraits to get the global picture and show attractors and their basins of attraction.
Jacobian analysis using the graphical Jacobian
For each equilibrium point (point of intersection of different null-clines) choose two points, one of which is located to the right and the other upward from the equilibrium. Find the components of the Jacobian using the following rule: the sign of the component of the vector field to the right of the equilibrium point gives the sign of and the sign of the component of the vector field to the right of the equilibrium point gives the sign of . The sign of the and vector field components upward of the equilibrium point give the sign of and . 2. 2.
Note that if one of the components is zero, then the corresponding derivative of the Jacobian will be zero. The latter can happen if at a given equilibrium point one or both null-clines are exactly horizontal or exactly vertical. 3. 3.
Try to compute the sign of the determinant and of the trace of the Jacobian and try to identify the type of equilibrium from fig.5.5. 4. 4.
If the type of equilibrium cannot be identified, use the analytical determinant-trace method formulated above. 5. 5.
Draw local and then global phase portraits of the system and make predictions about its dynamics.
7.2 Examples
Example 1 Study system (5.3), using null-clines and the determinant-trace method.
[TABLE]
Solution We have made studies of different aspects of this model throughout the reader, so let us collect the information.
I.Null-cline analysis.
The null-clines of this system are given in fig.6.3. We repeat it here in fig7.1a:
II.Jacobian analysis:
Equilibria. In example eq.(5.5) we found that this system has three equilibria points: and . 2. 2.
Jacobian. We compute the Jacobian as: , thus: J=\left(\begin{array}[]{lr}3-6x-1.5y&-1.5x\\ 0.5y&0.5x-0.25\end{array}\right).
Let us find equilibria types from the Jacobian.
At point the Jacobian is: J_{1}=\left(\begin{array}[]{lr}3&0\\ 0&-0.25\end{array}\right), , thus this is a saddle point.
At the point the Jacobian is: J_{2}=\left(\begin{array}[]{lr}-3&-1.5\\ 0&0.25\end{array}\right), , thus this is a saddle point.
At the point the Jacobian is: J_{3}=\left(\begin{array}[]{lr}-1.5&-0.75\\ 0.5&0\end{array}\right), , thus we need to find and thus this is a stable node. 3. 3.
Local phase portraits are presented in fig.7.1b. From the null-cline analysis we find the approximate locations of the manifolds of the saddle points and the direction of the trajectories around the stable node. 4. 4.
Global picture. There are no general rules to draw the global picture. However in the case of fig.7.1b we can expect that the non-stable manifold of the saddle point at will end at the other saddle point , the non-stable manifold of the saddle point as well as most of other trajectories should go to the stable node at as this is the only attractor here. Thus we get the phase portrait presented in fig.7.1c. We see that we will have a stable global attractor with the basin of attraction the whole region (except two axes ).
We see that the phase portrait which we got as a result of our study is qualitatively the same as the phase portrait of this system obtained using a computer (see fig.5.1b).
Let us try to apply for this problem the method of the graphical Jacobian
Example 2 Study the same system (7.7), using the graphical Jacobian method.
Solution From fig.7.1a we find in the following components of the Jacobian:
Point , J=\left(\begin{array}[]{lr}+\alpha&0\\ 0&-\delta\end{array}\right). Thus , thus we have a saddle. 2. 2.
Point , J=\left(\begin{array}[]{lr}-\alpha&-\beta\\ 0&+\delta\end{array}\right). Thus , thus we have a saddle. 3. 3.
Point , J=\left(\begin{array}[]{lr}-\alpha&-\beta\\ +\gamma&0\end{array}\right). Thus , , thus we have a stable equilibrium (stable node or stable spiral). 4. 4.
Drawing of the local and global phase portraits is exactly the same as with the previous method
Thus we see that in this case we were able to solve the problem completely using the graphical Jacobian method. The only difference with the determinant-trace method is that we were not able to say that the stable equilibrium is a node. However, this has almost no consequences for the phase portrait.
7.3 Exercises
Complete the vector field approximations for the null-clines shown in fig7.2a,b,c. Mark equilibria and determine their type using the graphical Jacobian method. Draw local and then global phase portraits. Try to describe the dynamics of the system, by saying what happens with the variables and in the course of time.
Draw the phase portrait of the following systems of differential equations. Explain their dynamics.
- (a)
\left\{\begin{array}[]{l}{dx\over dt}=x(1-x-y)\\ {dy\over dt}=y(1-2x)\end{array}\right.x\geq 0;\;y\geq 0 2. (b)
\left\{\begin{array}[]{l}{dN\over dt}=2N-NP-N^{2}\\ {dP\over dt}=3P-2NP-P^{2}\end{array}\right.N\geq 0;\;P\geq 0 3. 3.
In theoretical biology the following model has been used to study mRNA () protein () interaction:
[TABLE]
Study this system using the graphical Jacobian approach, i.e. draw null-clines, mark equilibria as points of intersection of the null-clines, determine stability of these equilibria and sketch a qualitative phase portrait. 4. 4.
The following equations describe the dynamics of predator () and prey () populations:
[TABLE]
here are parameters
This system of differential equations always has an equilibrium corresponding to an extinct population of the predator and non-zero population of prey ().
- (a)
Find this () equilibrium. 2. (b)
Find the Jacobian at this equilibrium point. 3. (c)
For which values of the parameters the predator population can be driven to extinction (i.e. to that equilibrium ).
(N.B. Do not use the ’graphical Jacobian approach for this problem!).
Additional exercises
Draw the phase portrait of the following systems of differential equations. Explain their dynamics.
- (a)
\left\{\begin{array}[]{l}{dx\over dt}=2y\\ {dy\over dt}=x-x^{2}-0.5y\end{array}\right. 2. (b)
\left\{\begin{array}[]{l}{dx\over dt}=x+y^{2}\\ {dy\over dt}=x+y\end{array}\right. 3. (c)
\left\{\begin{array}[]{l}{dx\over dt}=x(25-x^{2}-y^{2})\\ {dy\over dt}=y(x-3)\end{array}\right.x\geq 0;\;y\geq 0 4. (d)
\left\{\begin{array}[]{l}{dx\over dt}=xy\\ {dy\over dt}=4-y-x^{2}\end{array}\right.\quad-\infty<x<\infty;\quad-\infty<y<\infty 6. 6.
Lotka Volterra model with competition in the prey population:
[TABLE]
for . 7. 7.
A swimming pool is infested with algae whose population is . The owner attempts to control the infestation with an algicidal chemical, poured into the pool at a constant rate. In the absence of algae, the chemical decays naturally; when algae are present it is metabolized by them and kills them. The equations of the rates of change of and the concentration of the chemical in the pool, , are
[TABLE]
where are positive constants. Discuss the meaning of each term in these equations.
- (a)
Put and show, that the system has two non-negative equilibria if , find them, draw the phase portrait of this system and explain the dynamics. 2. (b)
What happens if ? 3. (c)
- Find the conditions for control of the infestation for arbitrary positive values of .
Chapter 8 Limit cycle
8.1 Stable and non-stable limit cycles
In previous chapters we found several possible types of equilibria: saddle, node, spiral and center. Some of these equilibria can be attractors which determine the ultimate dynamics of a system. However, it turns out that there exists another important attractor for systems of two differential equations. It is called a limit cycle. The geometrical image of a limit cycle on a phase portrait is a closed curve. The usual phase portrait of a system with a limit cycle is shown in fig.8.3a. Here the limit cycle is shown as a bold ellipse. The figure also shows two more trajectories: one inside the limit cycle and one outside the limit cycle. Inside the limit cycle you can see a spiral. This is not unusual. There is a theorem which states that for systems of two equations there is always at least one equilibrium point inside the limit cycle. In most of the cases such an equilibrium point is a spiral. Outside the limit cycle we can see a trajectory which approaches the limit cycle and winds onto it.
The limit cycle which is shown in fig.8.3a is called a stable limit cycle. This is because if we start on a trajectory close to this limit cycle, this trajectory will approach this limit cycle in the course of time.
There also exists another type of limit cycle called a non-stable limit cycle (fig.8.3b). If we start on a trajectory close to a non-stable limit cycle, it will diverge from this limit cycle. In order to distinguish the stable and non-stable limit cycles we will draw the non-stable limit cycle using a dashed line.
The main questions regarding the limit cycle are: what will be the dynamics of systems with limit cycles and how do limit cycles occur in systems of two differential equations? We will also consider an example of a biological system with a limit cycle.
8.2 Dynamics of a system with a limit cycle.
What will be the dynamics of systems with a limit cycle? In section 4.4.2 we studied an equilibrium point called a “center”. The phase portrait of that point (fig.4.8) was a set of embedded ellipses. The corresponding dynamics are oscillations of the and variables. Therefore, the dynamics which corresponds to the trajectory which starts on the limit cycle will also be oscillations. Fig.8.3a, fig.8.3a shows an example of dynamics of the variable for the trajectory which starts on the limit cycle (at the point C in fig.8.3).
The dynamics of trajectories originating around the limit cycles will depend on its type. If the limit cycle is stable, then a trajectory which starts inside the limit cycle will approach it in the course of time and we obtain oscillations with initially increasing amplitude (fig.8.3b). If the trajectory starts outside the limit cycle, then we obtain oscillations with an initial decrease of amplitude until the trajectory reaches the limit cycle (fig.8.3c).
For a non-stable limit cycle the dynamics are different. Only the trajectory which starts on a limit cycle will have oscillatory dynamics (fig.8.3a). Other trajectories will have different behavior. The trajectory which originates inside the limit cycle will approach the stable equilibrium inside the limit cycle and we obtain oscillations with gradually decreasing amplitude (fig.8.3b). The trajectory which starts outside the limit cycle blows up to infinity or to another equilibrium (fig.8.3c). The oscillations in fig.8.3a with a non-stable limit cycle are highly improbable in real systems. This is because for such oscillations the trajectory must start exactly on the limit cycle and even small disturbances will switch our system either to the behavior of fig.8.3b or fig.8.3c. Therefore, for real systems the non-stable limit cycle just determines the basin of attraction of the stable equilibrium which is located inside this limit cycle.
8.3 How do limit cycles occur?
In many cases the limit cycle occurs as a result of the changing of a parameter in a system of differential equations. The most usual process of formation of a limit cycle is the following. Assume we have a system of two differential equations with a parameter :
[TABLE]
Assume that at some parameter value system (8.1) has a global attractor which is a stable spiral (fig.8.4a). This means that all the trajectories which originate close to or even far from this equilibrium approach it in the course of time. Because our equilibrium point is a stable spiral the Jacobian of the system at this point will have two complex eigen values . However, because our system will now depend on the parameter the eigen values will also depend on this parameter:
[TABLE]
Because we have assumed that at system (8.1) has an equilibrium which is a stable spiral, then . When we change the parameter the value of will change and at some it can become a positive number . This means that the stable spiral in fig.8.4a will become an unstable spiral. However, as we found in chapter 4 (fig.4.3), the eigenvalues only give us the dynamics close to the equilibrium point of our system. Therefore, it can happen that although close to the equilibrium point our spiral becomes unstable, the global behavior far from the equilibrium remains the same, i.e. far from the equilibrium we still have a converging flow (fig.8.4b). Hence, in our phase portrait we have two types of flow: the diverging flow from the unstable spiral around the equilibrium and the converging flow from the periphery. Simple geometrical consideration shows that these two flows must be separated from each other. The line of separation will be the limit cycle in our system. Therefore, we will obtain a phase portrait as in fig.8.3a.
Note, that such a mechanism of limit cycle formation frequently occurs in systems of two equations. It is called the Hopf bifurcation. From our analysis it follows that the Hopf bifurcation occurs when changes its sign, i.e. at the parameter value where:
[TABLE]
8.4 Example of a system with a limit cycle
Let us consider the Holling-Tanner model for predator-prey interactions.
[TABLE]
here denotes the prey and the predator population, the term describes the growth of the prey in absence of predator, accounts for the predator-prey interaction. At system (8.4) has one equilibrium point (for ) which is a stable spiral. The null-clines, phase portrait and dynamics of this system are shown in fig.8.5.
If we increase the value of the parameter which accounts for the carrying capacity () the type of equilibrium changes and it becomes an unstable spiral. As we discussed in section 8.3 we expect the formation of a stable limit cycle. We can clearly see it in fig.8.6b. The trajectory, which starts from the same initial conditions as the trajectory from fig.8.5b will now approach some closed curve which is a limit cycle and the dynamics of the system will be oscillatory (fig.8.6c).
If we start a trajectory inside the limit cycle then, as we predicted in fig.8.3a and fig.8.3b, the trajectory will approach the limit cycle and the dynamics will be oscillation with increasing amplitude fig.8.7a,b. The complete phase portrait of this system at is shown in fig.8.7c.
You will learn more about the biological consequences of this type of dynamics in the course “Theoretical Biology” by R. de Boer.
Conclusion 13
A limit cycle is a closed trajectory on the phase portrait of a system of two differential equations.
If trajectories around the limit cycle converge onto it, then the limit cycle is called a stable limit cycle. If trajectories around the limit cycle diverge away from it, then the limit cycle is called a non-stable limit cycle.
The dynamics of a system with a stable limit cycle is oscillatory. The dynamics of a system with an non-stable limit cycle is either converging to the equilibrium which is located within the limit cycle or diverging, possibly to infinity.
Limit cycles can appear as a result of a Hopf bifurcation, i.e. the process where the real part of the complex eigenvalues change their sign.
8.5 Exercises
Assume that a system of two differential equations has two equilibria which are a saddle point and a non-stable spiral. The phase portrait of this system is partially shown in fig.8.8a What is missing here?
- (a)
Complete the phase portrait of this system 2. (b)
Qualitatively sketch the dynamics of the variable ( dependence of the variable on time) for the initial conditions which are shown in fig.8.8 by the point A in fig.8.8a. 2. 2.
Complete the phase portrait of this system shown in fig.8.8b,c. Qualitatively sketch the dynamics of the variable ( dependence of the variable on time) for the initial conditions which are shown by points A,B,C in fig.8.8b,c. (Note, that in fig.8.8b,c a stable limit cycle is shown by the solid line and a non-stable limit cycle is shown by the dashed line)
As we discussed in (8.3) the necessary condition for the appearance of a limit cycle via a Hopf bifurcation is that the eigenvalues of the Jacobian matrix of a system at the equilibrium point are complex and the real part of the eigenvalues are zero . Prove that this condition can be rewritten in the following way using the and of the Jacobian matrix:
[TABLE] 4. 4.
One of the classical models for oscillatory phenomena in biochemical systems is a model called the Brusselator. In dimensionless form this model can be written as the following system of two differential equations:
[TABLE]
here and are concentrations of two biochemical species and and are parameters. Study system (8.6) for .
- (a)
Find a non-trivial equilibrium. 2. (b)
Determine the stability of this equilibrium as a function of the parameter . 3. (c)
Find the value of when system (8.6) undergoes a Hopf bifurcation. (Note: it can be helpful to use equations (8.5)). 4. (d)
For which values of do you expect oscillations in system (8.6)? 5. 5.
Study the following predator-prey model
[TABLE]
- (a)
Draw null-clines of this system for and two values of , and . (Hint: the maximum of the parabola is reached at the middle between its roots (i.e. at )). 2. (b)
Determine the stability of non-trivial equilibrium in both cases using the graphical Jacobian. 3. (c)
For which value of do you expect oscillatory behavior? 4. (d)
Extend your analysis for arbitrary positive values of the parameters () provided . Find for which values of the non-trivial equilibrium is stable. When do we expect oscillations? (Note, that critical values of should be a function of the other parameters of the system).
Chapter 9 Historical notes
I developed this course in 1995-2010 for biology students from Utrecht University. The main idea was to develop a course which will allow students with minimal mathematical background to understands basic methods used in mathematical biology. Over the years the course was many times adjusted and modified. The current text contains the latest version of the course from 2010.
Chapter 10 Dictionary
[TABLE]
[TABLE]
Chapter 11 Hints:
11.0.1 Solution of the initial value problem for a linear system
Let us illustrate how we can find a solution for the initial value problem, on example of system (4.11).
Problem: Find the solution for the following initial value problem:
[TABLE]
Solution: We already found a general solution of (11.1) as formula (4.20):
[TABLE]
Using it we can find a particular solution, corresponding to any given initial conditions.
We proceed as follows. At time (11.2) gives:
[TABLE]
This expression will satisfy the initial conditions if: C_{1}\left(\begin{array}[]{c}-4\\ 2\end{array}\right)+C_{2}\left(\begin{array}[]{c}-4\\ -2\end{array}\right)=\left(\begin{array}[]{c}4\\ 6\end{array}\right). This gives the following system of equations for unknowns and : \left\{\begin{array}[]{l}-4C_{1}-4C_{2}=4\\ 2C_{1}-2C_{2}=6\end{array}\right.
We now solve this system using the method described in section 1.1.4. From the second equation we find:
, or . After substitution of this expression to the first equation we find:
, i.e. , or and .
We now find from our substitution as .
Thus give the solution of our system satisfying given initial conditions. Let us rewrite it as:
[TABLE]
thus the particular solution is given by , and .
11.0.2 Equilibria/derivatives
Problem:(A) Find equilibria of the following systems \left\{\begin{array}[]{l}{dx\over dt}=f(x,y)\\ {dy\over dt}=g(x,y)\end{array}\right. (see definition in section 5.1.2).
(B) Find the following partial derivatives at each equilibrium point ().
(a) \left\{\begin{array}[]{l}{dx\over dt}=4x-2xy\\ {dy\over dt}=2xy-4y\end{array}\right.
Solution: (A) We rewrite the first equation as . It has solutions or .
Let us substitute them to the second equation:
For , we get , thus . Therefore we found one equilibrium .
Now substitute . We get , thus . We found the second equilibrium .
(B) , thus .
, thus .
At the equilibiria points:
Equilibrium : .
Equilibrium : .
Chapter 12 Answers for selected exercises and Formulas lists
Exercises Chapter 1
-
-
(a)
2. (b)
2. 2.
- (a)
2. (b)
3. 3.
- (a)
; 2. (b)
; 3. (c)
; , or . 4. (d)
, or , thus . 5. (e)
, or 4. 4.
- (a)
From 1st eq. , substitution to 2nd eq. gives , thus . 2. (b)
from 1st eq. to 2nd eq. gives 3. (c)
From 2nd eq. , this or . Substitution to the 1st eq. gives or and . Now substitute , we get , thus all solutions are given: . 4. (d)
From 2nd eq. , thus or . Substitution to the 1st eq. gives:
. 5. (e)
From 2nd eq. , or . Substitution to the 1st eq. gives: : . 5. 5.
- (a)
2. (b)
and 3. (c)
4. (d)
6. 6.
- (a)
-intercept , zeros , horizontal asymptote: 2. (b)
-intercept , zeros and . 7. 7.
- (a)
fig.a 2. (b)
fig.b 3. (c)
horizontal asymptote , intercept . (fig.c). Dependence on . If increases, the horizontal asymptote will be approached at a slower rate. (fig.c)
-
-
(a)
(see below,fig.a) 2. (b)
(see below,fig.b)
- (c)
(see above, fig.c) 4. (d)
(see above, fig.d) 5. (e)
and (see below, fig.a) 6. (f)
and (see below, fig.b)
- (g)
and (see above, fig.c). 8. (h)
(see above, fig.d). 9. 9.
- (a)
2. (b)
10. 10.
no answer provided for this exercise 11. 11.
- (a)
, thus 2. (b)
, thus 3. (c)
and 4. (d)
and 5. (e)
6. (f)
7. (g)
8. (h)
9. (i)
12. 12.
no answer provided for this exercise 13. 13.
- (a)
, 2. (b)
, 3. (c)
14. 14.
- (a)
-intercept ,zero ,horizontal asymptote , as , vertical asymptote is at , i.e. , and . 2. (b)
thus -intercept , one zero: 15. 15.
- (a)
no answer provided for this exercise 2. (b)
, as we have ’+’ at the parabola is opened upward. (Graph fig.a) 3. (c)
horizontal asymptote , vertical asymptote , no zeros (fig.b)
- (d)
(fig.a below). Parameter shifts the horizontal asymptote. 5. (e)
no answer provided for this exercise 6. (f)
.
no answer provided for this exercise
Exercises Chapter 2
-
-
(a)
. 2. (b)
. 2. 2.
- (a)
\left(\begin{array}[]{lr}2&-4\\ 1&1\end{array}\right)\left(\begin{array}[]{c}x\\ y\end{array}\right)=\left(\begin{array}[]{c}3\\ 1\end{array}\right), . 2. (b)
\left(\begin{array}[]{lr}a&b\\ c&d\end{array}\right)\left(\begin{array}[]{c}x\\ y\end{array}\right)=\left(\begin{array}[]{c}0\\ -b\end{array}\right), . 3. 3.
- (a)
, \mathbf{v_{1}}=k\left(\begin{array}[]{c}-1\\ -1\end{array}\right); , \mathbf{v_{2}}=k\left(\begin{array}[]{c}-1\\ 1\end{array}\right), where is an arbitrary number. 2. (b)
, \mathbf{v_{1}}=k\left(\begin{array}[]{c}-4\\ 2\end{array}\right); , \mathbf{v_{2}}=k\left(\begin{array}[]{c}-4\\ -2\end{array}\right), where is an arbitrary number. 3. (c)
, \mathbf{v_{1}}=k\left(\begin{array}[]{c}-5\\ -2-i\end{array}\right); , \mathbf{v_{2}}=k\left(\begin{array}[]{c}-5\\ -2+i\end{array}\right), where is an arbitrary number. 4. 4.
- (a)
at it is 2; at it is 4; 2. (b)
, at it is -18; 3. (c)
, at it is ; at it is [math]; , at it is [math]; at it is . 4. (d)
, ; 5. (e)
, ; 6. (f)
, ; 7. (g)
and 8. (h)
5. 5.
- (a)
2. (b)
3. (c)
4. (d)
6. 6.
- (a)
2. (b)
7. 7.
- •
A=\left(\begin{array}[]{lr}1&2\\ 2&1\end{array}\right), ; D_{x}=\left(\begin{array}[]{lr}5&2\\ 4&1\end{array}\right), D_{y}=\left(\begin{array}[]{lr}1&5\\ 2&4\end{array}\right), , thus .
- •
By usual method: thus after subs into 2nd eq. we get , or and hence . 8. 8.
- (a)
, \mathbf{v_{1}}=k\left(\begin{array}[]{c}-6\\ -3\end{array}\right); , \mathbf{v_{2}}=k\left(\begin{array}[]{c}-6\\ 4\end{array}\right), where is an arbitrary number. 2. (b)
, \mathbf{v_{1}}=k\left(\begin{array}[]{c}-1\\ -1\end{array}\right); , \mathbf{v_{2}}=k\left(\begin{array}[]{c}-1\\ 7\end{array}\right), where is an arbitrary number. 9. 9.
Exercises Chapter 3
at , The double size at : . 2. 2.
. 3. 3.
- (a)
Phase portrait in fig.a. (below), (zeros of parabola and ), attractor , basin . 2. (b)
Phase portrait in fig.a. (below), attractor , basin . 3. (c)
Phase portrait in fig.c. (below), attractor , basin and , basin, .
- (d)
fig.d. (above), attractor , basin and , basin, . 5. (e)
fig.a (below) attractor basin ; fig.b (below), attractor , basin and attractor , basin .
- (f)
graph is shown in fig.c (above), attractor , basin and attractor , basin . 4. 4.
General consideration. Typical graphs are shown in the figure. shifts the graph upward. If the total shift is less then the minimum of the graph, nothing will change. If the shift is more than the minimum (fig.b), the will go to the right equilibrium.
For questions (a) (b) final is , for question (c) final is . (d) 5. 5.
The maximal yield if equals maximum of (see the last section of this chapter), which gives 6. 6.
(a) the general solution , ;
(b) ;
(c) ; 7. 7.
the steady state value is and the characteristic time . 8. 8.
stable equilibrium at , yield is given by , which has maximal value . 9. 9.
the last strategy is better as population is more stable: small decrease in population size is OK for the last strategy, but for the other case small decrease will result in the population extinction.
Exercises chapter 4
-
-
(a)
\left(\begin{array}[]{c}x\\ y\end{array}\right)=C_{1}\left(\begin{array}[]{c}-1\\ -1\end{array}\right)e^{-t}+C_{2}\left(\begin{array}[]{c}-1\\ 1\end{array}\right)e^{-3*t} 2. (b)
\left(\begin{array}[]{c}x\\ y\end{array}\right)=C_{1}\left(\begin{array}[]{c}1\\ 1\end{array}\right)e^{2t}+C_{2}\left(\begin{array}[]{c}1\\ -2\end{array}\right)e^{5t} 2. 2.
\left(\begin{array}[]{c}x\\ y\end{array}\right)=\frac{3}{2}\left(\begin{array}[]{c}2\\ -2\end{array}\right)e^{3t} 3. 3.
system is \left\{\begin{array}[]{l}\frac{dC_{1}}{dt}=-0.01C_{1}+0.01C_{2}\\ \frac{dC_{2}}{dt}=0.04C_{1}-0.04C_{2}\end{array}\right., the solution: \left(\begin{array}[]{c}C_{1}\\ C_{2}\end{array}\right)=-240\left(\begin{array}[]{c}-0.01\\ -0.01\end{array}\right)-60\left(\begin{array}[]{c}-0.01\\ 0.04\end{array}\right)e^{-0.05t}, or , and . 4. 4.
- (a)
det\left|\begin{array}[]{lr}1-\lambda&4\\ 2&3-\lambda\end{array}\right|=(1-\lambda)(3-\lambda)-8=\lambda^{2}-4\lambda-5=0, ,\mathbf{v_{1}}=k\left(\begin{array}[]{c}-4\\ -4\end{array}\right); ,\mathbf{v_{2}}=k\left(\begin{array}[]{c}-4\\ 2\end{array}\right), thus this is a saddle point (the corresponding phase portrait form the figure (a) below). 2. (b)
characteristic eq. is , ,\mathbf{v_{1}}=k\left(\begin{array}[]{c}1\\ 3\end{array}\right); ,\mathbf{v_{2}}=k\left(\begin{array}[]{c}1\\ 1\end{array}\right), i.e. a non-stable node (phase portrait see fig.b). 3. (c)
characteristic eq. is , , ,(no eigen vectors required). A non-stable spiral (qualitative phase portrait see fig.c).
- (d)
characteristic eq. is , ,\mathbf{v_{1}}=k\left(\begin{array}[]{c}-1\\ 1\end{array}\right), ,\mathbf{v21}=k\left(\begin{array}[]{c}-1\\ -1\end{array}\right) i.e. a stable node ( phase portrait see fig.a below). 5. (e)
characteristic eq. is , , ,(no eigen vectors required). A stable spiral (qualitative phase portrait see fig.b). 6. (f)
characteristic eq. is , , ,(no eigen vectors required). A center (qualitative phase portrait see fig.c). 7. (g)
characteristic eq. is , ,\mathbf{v_{1}}=k\left(\begin{array}[]{c}1\\ -3\end{array}\right); ,\mathbf{v_{1}}=k\left(\begin{array}[]{c}1\\ -1\end{array}\right), i.e. a saddle ( phase portrait see fig.d).
Characteristic equation is given by
det\left|\begin{array}[]{lr}-2-\lambda&-a\\ 3&-1-\lambda\end{array}\right|=(-2-\lambda)(-1-\lambda)+3a=\lambda^{2}+3\lambda+2+3a=0, . If , i.e. we have complex roots. Because the real part for complex roots is negative we have a stable spiral. If the roots are real (), then is always negative. The first root will be positive is , i.e. if the expression under the root we have a value more than ( ), this gives us , or . Thus if , then and we have a saddle point. In the interval we have , thus both real roots negative, and we have a stable node. Conclusion, if we increase from , we will first have a saddle point (unstable eq.) till , then the saddle will become stable node (stable eq.), until . If becomes more than we will have a stable spiral (stable eq.). Qualitative phase portrait can be plotted as in the previous exercise. 6. 6.
In linear system oscillation occur if the equilibrium type is a center. Characteristic equation is given by . The center occurs if the roots complex with a zero real part (). These conditions give , thus oscillations occur, if and . 7. 7.
(a) system is \left\{\begin{array}[]{l}\frac{dx}{dt}=-(a+c)x+by\\ \frac{dy}{dt}=ax-(b+e)y\end{array}\right.,
(b) system is \left\{\begin{array}[]{l}\frac{dx}{dt}=-5x+2y\\ \frac{dy}{dt}=0.5x-5y\end{array}\right., eigen values , , stable node, stable equilibrium. (c) formula for eigen values in general case are:
. We see that , thus expression under is less than , thus we can have either stable node or stable spiral, both types of equilibria are stable.
Exercises chapter 5
(a) Equilibria: . are:
at :;
at : . 2. (b)
Equilibria: . Derivatives are:
at :;
at : . 3. (c)
Equilibria: Derivatives are:
at : ; at : . 4. (d)
Equilibria: Derivatives are: at :;
at : ; at :;
at : . 5. 2.
From 2nd or , which after substitution to 1st equation gives 3 equilibria , ,. All non-negative if: 6. 3.
From 2nd . Substitution to 1st gives one non-negative equilibrium: and thus 7. 4.
From 2nd , or . Substitution to 1st gives equilibria: () and (. The first equilibrium is always positive, the second one is positive if . 8. 5.
- (a)
non-stable spiral, non-stable. 2. (b)
saddle, non-stable. 3. (c)
saddle, non-stable. 4. (d)
center, neutrally stable. 9. 6.
(a) 2 equilibria , which are stable node and saddle; (d) see fig.a below.
only one equilibrium . At this equilibrium, ,, thus equilibrium is stable. 11. 8.
(a) , linearization , where J=\left(\begin{array}[]{lr}0&-\frac{bd}{c}\\ \frac{ac}{b}&0\end{array}\right); (c) center point; (d) see fig.b above. 12. 9.
no answer provided for this exercise 13. 10.
no answer provided for this exercise 14. 11.
no answer provided for this exercise 15. 12.
no answer provided for this exercise 16. 13.
Equilibria of system. From 2nd , to 1st eq., , i.e. equilibrium is . Jacobian at the equilibrium is: , at it is , , , thus J=\left(\begin{array}[]{lr}-a&-1\\ \varepsilon&0\end{array}\right) and we get , thus equilibrium is always stable, and it is a node if , i.e. , and a stable spiral if .
Exercises chapter 6
(a) Null-clines are given below in fig.a, thus the graphical Jacobian on basis of ’black’ points is: J=\left(\begin{array}[]{lr}\alpha&\beta\\ -\gamma&-\delta\end{array}\right). This gives for , and we do not know what is the sign. Thus graphical Jacobian does not work here. The real Jacobian here is: J=\left(\begin{array}[]{lr}3&1\\ -1&-1\end{array}\right), that gives , thus we have a saddle point. Phase portrait is in fig.b.
- (b)
Null-clines are given above in fig.c, Graphical Jacobian does not work here. The real Jacobian gives stable spiral. Phase portrait is in fig.d. 3. (c)
Null-clines are given below in fig.a, the graphical Jacobian gives saddle point. Phase portrait in fig.b. 4. (d)
Null-clines are given below in fig.c, the graphical Jacobian gives a non-stable equilibrium (node, spiral). The real Jacobian gives a non-stable spiral. Phase portrait fig.d.
For fig.a below we have two equilibria (marked). For the left one the graphical Jacobian gives a stable equilibrium (stable node or stable spiral). The right equilibrium gives a saddle.
For fig.b we have one equilibrium (marked). The graphical Jacobian gives a stable equilibrium (stable node or stable spiral).
For fig.c we have one equilibrium (marked). The graphical Jacobian gives a stable node. 6. 3.
(a) See fig.a below. (b,c) equilibria: ,graphical Jacobian gives a stable node; : saddle. For phase portrait see solution problem 6 chapter 5.
(a) See fig.b above. (b,c) equilibria: ,graphical Jacobian: saddle; 2nd equilibrium , graphical Jacobian: center. Phase portrait fig.c above. (d) no, for positive parameters. 8. 5.
(a), (b) no answer provided for this exercise 9. (c)
vector field see fig.a; Equilibria ,graphical Jacobian: J=\left(\begin{array}[]{lr}\alpha&\beta\\ \gamma&-\delta\end{array}\right), ,saddle; 2nd equilibrium ,graphical Jacobian: J=\left(\begin{array}[]{lr}\alpha&-\beta\\ \gamma&-\delta\end{array}\right), , sign unknown, , equilibrium type unknown; phase portrait fig.b.
- (d)
no answer provided for this exercise
Exercises chapter 7
(a) graphical Jacobian: at 1 stable node. At 2 saddle. Null-clines fig.a above, phase portrait fig.d. 2. (b)
At 1 saddle. At 2 not known, stable node/spiral. At 3 saddle. Null-clines fig.b above, phase portrait fig.e. 3. (c)
At 1 not known, stable node/spiral. At 2 saddle. At 3 not known, stable node/spiral. Null-clines fig.c above, phase portrait fig.f. 4. 2.
(a) Null-clines : , : , fig.a below.
Equilibria . Graphical Jacobian: unstable node. : stable node. saddle. Phase portrait fig.b. One attractor . Basin of attraction shaded.
- (b)
Null-clines : , : , fig.c above
Equilibria . Graphical Jacobian: , unstable node; , stable node; , stable node: cannot determine equilibrium type using graphical Jacobian. Need ’real’ Jacobian, which gives saddle. Phase portrait fig.d. Two attractors and , basins of attraction shaded. 6. 3.
Null-clines : , : , fig.a below
One equilibrium. Graphical Jacobian gives stable equilibrium node/spiral. Phase portrait fig.b. One attractor 1. Basin of attraction shaded.
(a)Equilibrium ; (b) J=\left(\begin{array}[]{lr}r-\frac{2rN}{K}-bP&-bN\\ bP&bN-2b\end{array}\right)=\left(\begin{array}[]{lr}-r&-bK\\ 0&bK-2b\end{array}\right); (c) Equilibrium is stable if . 8. 5.
(a)-(d) no answer provided for this exercise 9. 6.
Fig c,d above
Exercises Chapter 8
see figure below
see figure below