APPENDIX 

SOME  MATHEMATICAL TOOLS FOR POPULATION BIOLOGISTS

Introduction

I assume a basic knowledge of elementary algebra and calculus. This appendix discusses some other mathematical tools which are frequently used in the development of ecological and evolutionary models, and shows how they can be implemented in Mathematica. Some worked  examples are shown. You should do these examples by activating the Mathematica input, and make up more examples of your own for any of this material which is unfamiliar to you.

A   Taylor series expansions

If f(x) is a function of x, it is often useful to linearize it about a value x* which is of particular interest. We can write
            f(x) ~ f(x*) + (x – x*) f '(x*)                                                        							               			[1]
where f '(x*) is the derivative of f(x) evaluated at x = x*. This approximation is valid for values of x sufficiently close to x*. It relies on the function being locally linear; the derivative is the slope of the curve at x*, so that the righthand side of [1] is the best linear approximation to the curve.
This linear approximation is given in Mathematica by the function Series[f[x], {x,x*, 1}], the '1' indicating that the approximation includes only the first (linear) term. For example, we find the linear approximation to ln x near x = 2 and assess its accuracy:

Series[Log[x],{x,2.,1}]

f[x_]=Normal[Series[Log[x],{x,2.,1}]]

Table[{x,Log[x],f[x]},{x,1.5,2.5,.1}]  //TableForm

The last term in Series shows that the approximation is accurate to order (x–2)2; it can be deleted to give an ordinary expression by using Normal.
          The linear approximation [1] depends on the linearity of the function near x*, and will break down as soon as x moves far enough from x* for any curvature to become appreciable. One might hope that a better approximation could be obtained by taking into account higher powers of (x–x*), leading to a power series expansion of the form
          f(x) ~ a0 + a1 (x–x*) + a2 (x–x*)2 +  a3 (x–x*)3 + ...
For the power series to have the same derivatives as f(x) at x=x* requires that ar = f (r)(x*)/r!, where f (r)(x*) denotes the rth derivative of f(x). Thus the power series is

f(x) ~ f(x*) + (x–x*)f '(x*) + (x–x*)2  f ''(x*) + (x–x*)3 f '''(x*) + ...                  [2]
                                                            2!                         3!
This is known as a Taylor series expansion. By taking more terms into account, one expects to get a better and better approximation. The Taylor series up to terms of order (x–x*)n can be generated with Series[f[x],{x,x*,n}].

Find the Taylor series expansions to ln x near x = 2 including (a) quadratic and (b) cubic terms and compare their accuracy with that of the linear approximation.

Three important Taylor series expansions about x* = 0 are those for the exponential, sine and cosine functions. The divisors are the appropriate factorials:

Series[Exp[x],{x,0,5}]

Series[Sin[x],{x,0,5}]

Series[Cos[x],{x,0,5}]

This technique can be extended to a function of several variables, f(x), where x = (x1, x2, ..., xn), expanded about the point x = x*. The best linear approximation is
          f(x) ~ f(x*) + S (xi - xi*) df/dxi                                                                          [3]
where df/dxi is the partial derivative of f with respect to xi evaluated at x = x*. For example, we find the linear approximation to ln x1/x2 near x1 = 2, x2 = 1:

f[x1_,x2_]=Normal[Series[Log[x1/x2],
				{x1,2.,1},{x2,1.,1}]]

f[1.9,1.1]

Log[1.9/1.1]

B  Complex numbers

This section is addressed to readers who are unfamiliar with the use of complex numbers, which arise frequently in stability analysis in determining whether oscillatory behavior is likely to occur. The quadratic equation
																				ax2 + bx + c = 0
has two solutions
																				x = [–b + (b2 – 4ac).5]/2a and x = [–b – (b2 – 4ac).5]/2a
When b2 < 4ac, this presents the problem of finding the square root of a negative number. A similar problem can occur in solving higher order polynomial equations, and in other contexts. It is convenient to resolve this problem by inventing an "imaginary" number i whose square is –1, in the same vein as earlier mathematicians invented negative numbers. In Mathematica i  is writtten I in accordance with the convention that inbuilt quantities begin with capitals. Thus

Solve[x^2 - 2x + 5 ==0]

A number like 1 – 2i , with a real part (1) and an imaginary part (–2) are called complex numbers.
										Complex numbers provide a consistent system when manipulated by the ordinary laws of algebra, provided that i2 is replaced by –1 whenever it occurs. For example,
										(1 + 2i) + (2 – i) = 3 + (2–1)i = 3 + i
										(1 + 2i)(2 – i) = 2 – i + 4i -2i2 = 4 + 3i

z1 = 1 + 2 I;
z2 = 2 - I;

z1 + z2

z1  z2

The ratio z3 = z1/z2 is the number such that z2 z3 = z1. In the above example, if z3 = x + iy
										(2 – i)(x + iy) = (2x + y) +(–x + 2y)i = 1 + 2i
										2x + y = 1; –x + 2y = 2
										x = 0, y = 1
										z3 = i

z3=z1/z2

z2 z3

A complex number has two parts, real and imaginary, which can be represented in a two-dimensional diagram (the Argand diagram) with the real part along the horizontal axis and the imaginary part along the vertical axis. (See Figure B1 in the book.) It is sometimes convenient to represent a complex number by its polar coordinates in the Argand diagram, that is to say its modulus or absolute value  (its distance from the origin) and its argument (the angle between the line joining the point to the origin and the real axis):

r1=Abs[z1] //N

theta1=Arg[z1] //N

r2=Abs[z2]  //N

theta2=Arg[z2]  //N

From the diagram in the book, for any complex number z = x + iy with r = Abs(z), q = Arg(z), then 
																			x = r cos q
																			y = r sin q
																			z = x + iy = r( cos q + i sin q)

r1 Cos[theta1]

r1 Sin[theta1]

It is possible to define functions of a complex variable by analogy with the corresponding definition for a real variable. Of particular importance is the exponential function, which can be defined by the series expansion
          	exp(z) = 1 + z + z2/2! + z3/3! + ...

Exp[z1]//N

Use the following function to evaluate the series expansion for exp(z1) for different values of n.

approxexp[z_,n_]:=Sum[z^i/i!,{i,0,n}]//N

If z1 and z2 are complex numbers, it can be verified from the series expansion definition that
          exp(z1 + z2) = exp(z1) exp(z2)
a well known property for real variables. For example,

Exp[z1+z2]  //N

Exp[z1] Exp[z2] //N

In particular, if z = x + iy
          exp(z) = exp(x + iy) = exp(x) exp(iy)
The exponential of an imaginary number is
          exp iy = 1 + iy + i2y2/2! + i3y3/3! + ...
                     = (1 – y2/2! + y4/4!  ...) + i(y  – y3/3! + y5/5! ...)
                     = cos y + i sin y
 Hence 
 										z = r(cos q + i sin q) = r exp iq

z1

r1 Exp[I theta1]

It follows that
          z1 z2 = r1 exp(i q1) r2 exp(i q2) = r1 r2 exp i(q1 + q2)
Thus the rule for multiplying complex numbers is simple in polar coordinates: multiply the absolute values and add the arguments.

Abs[z1 z2]

r1 r2

Arg[z1 z2]  //N

theta1 + theta2

C   Elementary vector operations

A vector is an ordered list of numbers. For example, agefem below is the vector giving the numbers of breeding female Great Tits of ages 1 through 8 shown in Table 4.2, and agemale is the corresponding vector for males:

agefem={92,37,22,12,6,3,0,0};
agemale={83,67,31,31,16,7,3,2};

Geometrically, a vector of length n (i.e. with n elements) can be thought of as a point in n-dimensional space; it may also be thought of as a directed line from the origin to the point, the line having a magnitude (the distance from the origin) and a direction. These two ways of representing the vector {1.5, 1} in two dimensions are shown in Figure C1 in the book.

To demonstrate elementary arithmetic operations on vectors, I define two symbolic vectors:

v1={a1,a2,a3};
v2={b1,b2,b3};

Addition of two vectors of the same length is defined by adding their corresponding elements:

v1+v2

For example, the age distribution of Great Tits of both sexes is

agefem + agemale

Multiplication by a scalar (an ordinary number) is defined by multiplying each element by that number:

2 v1

The inner product of two vectors of the same length is a scalar obtained by multiplying corresponding elements and then summing these products. It is often called the dot product, and is represented in Mathematica by a dot (.):

v1.v2

Geometrically, the inner product of a vector with itself is the square of its magnitude:

v1.v1

If the magnitudes of v1 and v2 are r1 and r2 and the angle between them is q, it can be shown that
          v1.v2 = r1 r2 cos q
This can be used to find cos q and hence q, as we now do for the age distributions of female and male Great Tits:

agefem.agemale/Sqrt[agefem.agefem*agemale.agemale]//N

ArcCos[%]

(This is in radians; if you want the answer in degrees divide it by Degree and take the numerical value.) If two vectors have the same direction, so that one is a multiple of the other, q = 0 and cos q = 1; if the two vectors are at right angles (orthogonal), q = p/2 and cos q = 0.
	         Another use of the inner product is illustrated by the following calculation of the mean age of the female Great Tits:

agefem.Range[8]/Sum[agefem[[i]],{i,8}]  //N

The outer product of two vectors v1 and v2 of lengths m and n is a matrix of order m x n whose ijth element is the product of the ith element of v1 and the jth element of v2. It will be written v1 x v2 and is obtained in Mathematica as follows:

Outer[Times,v1,v2]

D   Elementary matrix operations

Matrix algebra was invented by Cayley in the last century to provide a compact way of writing simultaneous linear equations. A matrix is a rectangular array of numbers, which is represented in Mathematica as a list of lists. For example, let us make a list of the two lists agefem and agemale defined above:

age={agefem, agemale}

This can be displayed in conventional matrix format:

MatrixForm[age]

Thus age is a 2 x 8 matrix (with 2 rows and 8 columns), the rows representing females and males and the columns representing age classes.
           We now set up some matrices to use as examples in discussing their properties.

A1={{1,3,4},{9,-1,6}};
A2={{3,-1,2},{10,3,4}};
Aa={{a11,a12},{a21,a22}};
Ab={{b11,b12},{b21,b22}};

Addition of two matrices  of the same order is defined by adding their corresponding elements.

Aa     //MatrixForm
Ab     //MatrixForm
Aa+Ab  //MatrixForm

A1     //MatrixForm
A2     //MatrixForm
A1+A2  //MatrixForm

The transpose of a matrix A, usually denoted AT, has the rows and columns transposed:

A1

Transpose[A1]

%  //MatrixForm

Multiplication of a matrix by a scalar is defined by multiplying each element by the scalar.

1.5 A1 //MatrixForm

Matrix multiplication is an extension of the vector dot product. Suppose A and B are matrices of orders m x n and n x p respectively. Their product A.B is defined as the matrix with dimensions m x p, whose ijth element is the dot product of the ith row of A and the jth column of B.

Aa.Ab    //MatrixForm

Transpose[A1].A2   //MatrixForm

Notice that the matrix product A.B is only defined when the number of columns of A is equal to the numbers of rows of B. If A and B are square matrices with the same dimensions, both A.B and B.A are defined, but they are not usually the same; that is to say, matrix multiplication is not commutative:

{{1,3},{2,-1}}.{{4,6},{-2,8}}

{{4,6},{-2,8}}.{{1,3},{2,-1}}

[NOTE: A vector in Mathematica is simply a list, and there is no distinction between row and column vectors. You can use "dot" for either left or right multiplication of vectors by matrices and it carries out whatever operation is possible. If A is a vector of order m x n and v is a vector of length n , A.v is a vector of length m obtainedby multiplying the rows of A into v; if v is of length m, v.A is a vector of length n obtained by multiplying v into the columns of A:

A1.{x1,x2,x3}

{x1,x2}.A1

In matrix algebra it is usual to distinguish between row vectors and column vectors. A row vector is really a matrix with one row, and likewise a column vector is a matrix with one column, and if they are to be used in Mathematica they must be defined in this way:

v1={1,2,3};
v2={4,5,6};
rowv={{1,2,3}};
colv={{4},{5},{6}};

Note that

colv.rowv

is the same as the outer product

Outer[Times,v2,v1]

However

rowv.colv

is a matrix with one row and one column whereas

v1.v2

is a scalar.]

The motivation for defining matrix multiplication like this is that it leads to useful applications.  Consider the problem of solving 3 linear equations in 3 unknowns:

3x1 + 4x2 –  x3 = 4
															4x1 +  x2 + 12x3 = 12
																9x1 – 3x2 + x3 = 6

This can be written in matrix terminology more compactly as

A.x = c

where A is the matrix of known coefficients, x is the vector of unknowns and c is the vector of known constants on the right hand side:

A = {{3,4,-1},{4,1,12},{9,-3,1}};   
c={4,12,6};

Note the value of  A.x:

A.{x1,x2,x3}  //ColumnForm

It is not only more convenient to use this compact terminology, but also leads to new methodology such as solution of the equations by using the inverse matrix:

																																								x = A–1.c

Inverse[A].c  //N

This is the solution for x as will be explained in the next section.

E   Determinants, linear independence and matrix inversion

Suppose A is a square matrix. Its determinant is a scalar which allows one to determine whether there are any linear relationships among its rows or columns. If such relationships exist the determinant is zero and the matrix is called singular; otherwise the determinant is non-zero and the matrix is called non-singular.

A3={{1,2,4},{3,1,2},{7,4,8}};
A4={{1.,2,3},{2.,1,2},{3.,4,5}};

Det[A3]

Det[A4]

A3 is singular because the third row is equal to twice the second row plus the third row. The rows (and columns) of A4 on the other hand are linearly independent; no such linear relationship exists between them.
	         The definition of a determinant can be seen by calculating it for symbolic matrices of different dimensions. Array[a,{m,n}] generates a matrix of dimensions {m, n} with a[i, j] in the ijth position.

Det[Array[a,{2,2}]]

Det[Array[a,{3,3}]]

Each term contains one element  from each row and from each column. The column subscripts in each term are a permutation of the integers {1, 2, ..., n} and the sum (with a + or – sign) is taken over all possible permutations, of which there are n!. A permutation is classified as odd or even depending on whether it can be obtained from the natural ordering by an odd or an even number of reciprocal interchanges (this can be determined with Signature), a + sign being given to even permutations and a – sign to odd ones.
          The following properties follow almost immediately from this definition. If all the elements of a row (or column) are multiplied by a constant, the determinant is multiplied by the same constant. If two rows (or columns) are interchanged the sign of the determinant changes (because even permutations become odd and vice versa) but its absolute value is unchanged. If two rows (or columns) are identical, the determinant is zero (because interchanging them leaves the matrix the same). Hence for any square matrix, adding one row (or column) to another leaves the determinant unchanged. It follows that if the rows are linearly dependent the determinant is zero. (Suppose that there is a linear dependence involving the first row, for example. Then by adding multiples of the other rows to it, it is possible, without changing the determinant, to make this row zero; a matrix with a zero row has zero determinant.) The converse is also true: if the rows are independent then the determinant is non-zero. Thus there is a one-one correspondence between the dependence or independence of the rows (and columns) and whether the determinant is zero or non-zero.
          If A is a non-singular matrix with dimensions {n, n} it has a unique inverse matrix A–1 such that A.A–1 = A–1.A = In, where In is the identity matrix, a square matrix with dimensions {n, n} with all its diagonal elements equal to 1 and all off-diagonal elements zero. Note how the identity matrix and the inverse are obtained in Mathematica:

IdentityMatrix[4]  //MatrixForm

B=Inverse[A4]

A4.B  //MatrixForm

B.A4

A singular matrix cannot have an inverse. Suppose that A.B = In,  write ri for the ith row of A and cj for the jth column of B, and suppose also that A is singular with r3 = r4 + 2r5. Then r3.c3 = 1 and (r4 + 2r5).c3 = r4.c3 + 2r5.c3 = 0. This contradicts the assumption that r3 = r4 + 2r5.

Solution of linear equations

We return to the problem of solving n linear equations in n unknowns:

where A is a matrix of known coefficients aij, x is the vector of unknowns and c is a vector of known constants. 
          Consider first the non-homogeneous equation with c = /  0. If A is non-singular it has a unique inverse; pre-multiplying both sides by this inverse we find the unique solution

If A is singular, there are linear relationships among its rows. If these relationships do not hold among the c's, the equations are inconsistent and have no solution; if the relationships among the rows of A also happen to hold among the c's, a range of solutions exists. For example, the equations

have no solution whereas the equations

are compatible with any values of x1 and x2 on the line x2 = 1 – x1.
          Consider finally the homogeneous equation

If A is non-singular, premultiplication by its inverse shows that the only solution is the trivial solution x = 0. If A is singular, then any linear relationship among the rows of A automatically holds on the right-hand side so that a range of non-trivial solutions exists. For example, the equations

are satisfied by any values on the line x2 = –x1.

F   Eigenvalues and eigenvectors

If A is a square matrix of dimensions n x n and x is a vector of length n, then A.x is also a  vector of length n. Thus A can be regarded as defining a linear transformation on points in n-dimensional space. Consider a simple example in two-dimensional space:

A = {{1,2},{3,2}};

A.{0,1}

A.{1,1}

A has transformed the point {0,1} into the point {2,2} and the point {1,1} into {3,5}.
          A common application of this idea occurs in the set of linear difference equations:
          n(t+1) = A.n(t).
In this context n(t) might be the vector of numbers of animals of different ages in year t, and A  defines how this year's numbers are transformed into next year's (the Leslie matrix, Chapter 3).
          We shall now consider an important tool in studying linear transformations. If l is a scalar and x a non-zero vector such that
          A.x = lx                                                                                  [1]
then l is called an eigenvalue of A and x is the corresponding eigenvector. This equation can be written
          (A - lIn).x = 0.                                                                                                   [2]
This homogeneous equation can only have a non-zero solution if
          det(A - lIn) = 0.                                                                                                   [3]
The determinant defines a polynomial of degree n in l and the equation will in general have n solutions, though there may be fewer if there are repeated roots.  In practice there are usually n distinct eigenvalues, and for simplicity we shall from now on assume that this is the case. These eigenvalues can be found thus:

Eigenvalues[A]

Thus there are two eigenvalues, –1 and 4, of the 2 x 2 matrix A.
          To each eigenvalue there corresponds an eigenvector. They can be found with Eigenvectors, but there is no guarantee that they will be listed in the same order as the eigenvalues, and it is better to use Eigensystem  which lists first the eigenvalues and then the corresponding eigenvectors. The following assignment sets vals to the list of eigenvalues, and vecs to the list of eigenvectors:

{vals,vecs}=Eigensystem[A];

vals

vecs

The eigenvector {–1, 1} corresponds to the eigenvalue –1, and the eigenvector {2/3, 1) to the eigenvalue 4. In fact {2, 3} or {4, 6} or any multiple thereof is also an eigenvector corresponding to the eigenvalue 4; an eigenvector is arbitrary up to a multiplicative constant as is clear from the equation defining it. We now verify that 4 is an eigenvalue of A with the associated eigenvector {2, 3} or any multiple thereof:

A.{2,3}

A.{20,30}

An eigenvector can be thought of as a direction in n-dimensional space which is conserved under the linear transformation A while the magnitude  of the vector is multiplied by l. In the set of difference equations
          n(t+1) = A.n(t),
if the population vector n lies in the direction of a particular eigenvector this year, it will lie in the same direction next year but will be multiplied in magnitude by l; that is to say, all the elements of n will be multiplied by the same constant l. Linear difference equations will be considered in more detail in the section 6.
          So far we have been considering right eigenvalues and vectors. One can define in the same way left eigenvalues and vectors satisfying
          x.A = lx.                                                                                                            [4]
The left eigenvalues satisfy the same determinantal equation as the right eigenvalues and are therefore the same, but the corresponding eigenvectors are different. Mathematica only calculates right eigenvectors, but the left eigenvectors can be found as the right eigenvectors of the transpose of the matrix since (x.A)T = AT.x. Eigenvectors are arbitrary up to a multiplicative constant, and it is convenient to standardize each of the left eigenvectors so that its product with the corresponding right eigenvector is unity. As an example we find the standardized left eigenvectors of A:

{lvals, lvecs}=Eigensystem[Transpose[A]];

lvals

Do[lvecs[[i]]=lvecs[[i]]/(lvecs[[i]].vecs[[i]]),
				{i,2}]

lvecs

A check should be made that the eigenvalues are listed in the same order in vals and lvals; if not, the left eigenvectors should be rearranged to appear in the same order as the right eigenvectors.

An important fact is that a left and a right eigenvector corresponding to different eigenvalues are orthogonal. Suppose that e1 is a right eigenvector corresponding to l1 and that e2* is a left eigenvector corresponding to l2. Then
										e2*.A.e1 = e2*.(l1e1) = (l2e2*).e1
It follows that e2*.e1 = 0 if l1 differs from  l2.  We have standardized the left eigenvectors so that ei*.ei = 1.If we define R = {e1, e2, ... , en} as the matrix of right eigenvectors and L = {e1*, e2*, ... , en*} as the matrix of standardized  left eigenvectors, it follows that L is the inverse of RT:
          L.RT = In                                                                                                             [5]
Since a left inverse is also a right inverse
          RT.L =  In                                                                                                            [6]
We now verify this for the matrix A:

lvecs.Transpose[vecs]

Transpose[vecs].lvecs

Finally, we consider the effect of a small perturbation of the elements of a matrix A on an eigenvalue l with right and left eigenvectors e and e*. A perturbation from A to A+dA will change l to l+dl and e to e+de such that
          	A.e = le
         (A+dA).(e+de) = (l+dl)(e+de)
Ignoring second order terms
										dA.e + A.de = lde +	dle
Multiplying both sides by e* and remembering that e*.A = le* and that e*.e = 1:
          e*.dA.e = dl                                                                                                       [7]
If only one element in A is perturbed, say aij  to aij + daij, then
          dl = ei* ej daij                                                                                                      [8]
where ei* and ej are the ith and jth elements of e* and e respectively. This result is useful in sensitivity analysis discussed in Chapters 4 and 5.

G   Spectral decomposition of a matrix

It is often useful to represent a matrix in a particular way known as its spectral decomposition. We first observe that  
         A.RT = RT.L                                                                                                     [1]
where R is the matrix of right eigenvectors defined in Section F and L is a diagonal matrix with the eigenvalues on the diagonal and zero everywhere else. This follows from the definitions of eigenvalues and eigenvectors together with the fact that post-multiplying a matrix by a diagonal matrix multiplies the ith column of that matrix by  the ith diagonal element, so that the ith column of the right hand side is the ith eigenvalue multiplied by the ith eigenvector. Verifying this for A defined above:

A.Transpose[vecs]

Lambda = DiagonalMatrix[vals]

Transpose[vecs].Lambda

Post-multiplying [1] by L and using [F6] we find that
          A = RT.L.L                                                                                                          [2]
 This equation can be expressed as

A =  S li Ti                                                                                                                                                                 [3]                                                                                                               
where                                                                                                                               
           Ti = ei x ei*,

the outer product of the right and left eigenvectors. This is called the spectral decomposition  of A. To verify that [3] follows from [2] we set up two arbitrary 3 x 3 matrices R and L and a diagonal matrix with 1 in the leading position and zero everywhere else which will isolate the coefficient of l1:

R={{r11,r12,r13},{r21,r22,r23},{r31,r32,r33}}

L={{l11,l12,l13},{l21,l22,l23},{l31,l32,l33}}

Diag=DiagonalMatrix[{1,0,0}]

Transpose[R].Diag.L //MatrixForm

Outer[Times,R[[1]],L[[1]]] //MatrixForm

The significance of the spectral decomposition [3] is that 
          Ti.Ti = Ti
          Ti.Tj = 0, i /=  j                                                                                                        [4]
(since Ti.Tj =  ei x ei*. ej x ej*, and ei*. ej = 1 if i = j and 0 otherwise.) Thus the kth power of A is
           Ak = (S li Ti )k =   S lik Ti                                                                                [5]
This is useful for determining the behavior of large matrix powers. If l1 dominates the other eigenvalues in the sense that it is larger than any of them in absolute value, then [5] will be dominated by this eigenvalue when k is large so that
            Ak ~   l1k T1                                                                                                       [6]

Do[T[i]=Outer[Times,vecs[[i]],lvecs[[i]]],{i,2}]

T[1]

T[1].T[1]

T[1].T[2]

MatrixPower[A,10]

vals[[1]]^10 T[1] + vals[[2]]^10 T[2]

vals[[2]]^10 T[2]  //N

This idea can be extended to other functions of matrices. The exponential of a matrix, exp(A), can be defined by the power series appropriate for the exponential of a real variable (see Section A): 
          exp(A) = I + A + A2/2! + A3/3! + A4/4! + ...                                                     [7]
Using [5], we find that 
          exp(A) =   S (1 + li +  li2/2! + li3/3! + ...)Ti  =    S exp(li) Ti                         [8]                                           
We shall be interested in exp(At) where t is time:
          exp(At) =   S exp(lit) Ti  ~  exp(l1t) T1                                                               [9]

MatrixExp[3.5 A]  //MatrixForm

Exp[3.5 vals[[1]]] T[1] + Exp[3.5 vals[[2]]] T[2] //
										MatrixForm

Exp[3.5 vals[[2]]] T[2] 	//MatrixForm

Another application is to find the square root of a matrix A, defined as a matrix B such that B2 = A. This can be calculated as
          B = S li0.5 Ti

H   Linear differential equations and local stability analysis in        continuous time

Linear differential equations

The linear differential equation
          dx/dt = a x                                                                                                            [1]
expresses the model that the rate of change of x is a constant multiple a of its current value. Starting with the value x(0) at time zero, the solution is
										x(t)= x(0) exp (at)                                                                                                 [2]

Verify that this is the solution in two ways, (a) by using DSolve, (b) by showing that it satisfies both the initial condition and the differential equation.

Consider now the set of n linear differential equations in n variables:
          dx1/dt = a11 x1 + a12 x2 + ... + a1n xn
          dx2/dt = a21 x1 + a22 x2 + ... + a2n xn        
          .                                                                                                                             [3]
          .
          dxn/dt = an1 x1 + an2 x2 + ... + ann xn
which can be expressed in matrix notation as
          dx/dt= A.x                                                                                [4]
where differentiation of a vector means differentiating each of its elements. The solution subject to the initial values x(0) is 
           x(t) = exp(At).x(0)                                                                                            [5]
This satisfies the initial conditions and it also satisfies the differential equations since d exp(At)/dt = A exp(At). Using the spectral decomposition in Equation[G9] and writing

ci = ei*.x(0)                                                                                                        [6]

we find that

x(t) = S ci exp(lit) ei                                                                                            [7]

(since S exp(lit) Ti.x(0)  = S exp(lit)(ei x ei*) .x(0) = S exp(lit) ei x (ei* .x(0)) 
	                                                                                       = S ci exp(lit) ei ).

Example: Solve
          dx1/dt = x1 + 2x2
          dx2/dt = 3x1 + 2x2
with x1(0) = x2(0) = 1.

x1[0]={1,1};
A1={{1,2},{3,2}};

{vals1,vecs1}=Eigensystem[A1]

{lvals1,lvecs1}=Eigensystem[Transpose[A1]]

Do[lvecs1[[i]]=lvecs1[[i]]/(vecs1[[i]].lvecs1[[i]]),
					{i,2}]

c1=lvecs1.x1[0]

x1[t]=Sum[c1[[i]] Exp[vals1[[i]] t] vecs1[[i]],{i,2}]

Local stability analysis of nonlinear differential equations

Consider the differential equation
          dx/dt = f(x).
At equilibrium there is no change so that f(x) = 0; a value satisfying this condition is called an equilibrium (or stationary) value. Such a value is locally stable if the system tends to return to it after a small perturbation; otherwise it is locally unstable. Local stability can be established by linearizing about the equilibrium value in a Taylor series expansion.
          Suppose that x* is an equilibrium value so that f(x*) = 0, and consider what happens when x is near this value, so that x = x* + p, where p  is a small perturbation. Expanding f(x) in a Taylor series expansion about x*,
          f(x) ~ f(x*) + f '(x*) p.
Also
         dx/dt = d(x* + p)/dt = dp/dt
since x* is a constant. Thus to a linear approximation
          dp/dt =  f '(x*) p
whose solution is
          p(t) = p(0) exp(f '(x*) t).
The behavior of the system depends on the sign of the derivative f '(x*). If it is negative p(t) tends to zero so that the system returns to x* after a small perturbation; the equilibrium is locally stable. If it is positive, p(t) increases exponentially and the system moves away from x*, at least initially, after a small perturbation; the equilibrium is locally unstable.
          Consider for example the logistic model
          dn/dt = r(1 – n/K)n
We find the equilibria and evaluate the derivative at each of them:

f[n_]:=r (1 - n/k) n

equil=Solve[f[n]==0,n]

D[f[n],n]  /.equil

The equilibria are at n = K and n = 0; the first is stable and the second unstable.
          This analysis can be extended to a set of n differential equations in n variables:
          dx1/dt = f1(x1, x2, ..., xn)
          dx2/dt = f2(x1, x2, ..., xn)
          .
          .
          dxn/dt = fn(x1, x2, ..., xn)
Write x for the vector of variables. An equilibrium satisfies
          fi(x) = 0, i = 1, 2, ..., n.
Suppose that x* is an equilibrium value and consider a small perturbation from it so that x = x* + p, where p is a vector of perturbations. Then
          fi(x) ~ Sj (dfi(x*)/dxj) pj
Hence the set of differential equations governing the perturbations is approximately
          dp/dt = Jp
where J is the Jacobian matrix whose ijth element is dfi(x*)/dxj. The asymptotic behavior of the solution to  this equation is determined by the eigenvalue of J with the largest real part, say l = a + bi. If a is negative (that is to say, if all the eigenvalues have negative real parts), the perturbations will diminish and the equilibrium is locally stable; otherwise it is locally unstable. If l is real (b = 0), the perturbations will converge to (or diverge from ) 0 exponentially; if  l is complex, the perturbations will tend to oscillate about the equilibrium (with decreasing or increasing amplitude depending on the sign of a).
          Consider for example the Lotka-Volterra prey-predator model:
          dn1/dt = (r1 – a1n2) n1
          dn2/dt = (–r2 +  a2n1) n2
The equilibrium with both species present is at (n1, n2) = (r2/a2,  r1/a1). We find the eigenvalues of the Jacobian at this point:

f[1][n1_,n2_]:=(r1-a1 n2) n1
f[2][n1_,n2_]:=(-r2+a2 n1) n2

Jac=Table[D[f[i][n[1],n[2]],n[j]],{i,2},{j,2}]  /.
					{n[1]->r2/a2,n[2]->r1/a1}

Eigenvalues[Jac]

J   Linear difference equations and local stability analysis in discrete time

Linear difference equations

The linear difference equation
          x(t+1) = a x(t)                                                                                                      [1]
expresses the model that the value of x  next time is a constant multiple a of its value this time. Starting with the value x(0) at time zero, we find
          x(1) = a x(0)
          x(2) = a x(1) = a2 x(0)                                        
and so on. The solution for any time t is clearly
          x(t) = at x(0)                                                                                                         [2]
We can use NestList to demonstrate this:

f[x_]:=a x
NestList[f,x[0],10]

Solve  the linear difference equation with a = .5, 1.5, -.5, –1.5 with n(0) = 100 in each case for 10 generations.

Consider now the set of n linear difference equations in n variables:
          x1(t+1) = a11 x1(t) + a12 x2(t) + ... a1n xn(t)
          x2(t+1) = a21 x1(t) + a22 x2(t) + ... a2n xn(t)                                                   [3]
          .
          .
          xn(t+1) = an1 x1(t) + an2 x2(t) + ... ann xn(t)
which can be expressed in matrix notation as
          x(t+1) = A.x(t)                                                                                                    [4]
where x(t) is the vector {x1(t), x2(t), ..., xn(t)} and A is the matrix of coefficients aij. We want to solve this set of equations starting from the initial value x(0) at time zero. By analogy with the univariate case the solution is
          x(t) = At.x(0)                                                                                                      [5]
Using Equation [G5] and writing
          ci = ei*.x(0)                                                                                                         [6]
we find

x(t) =   S ci lit ei                                                                                               [7]

(since as before S lit Ti.x(0)  = S lit(ei x ei*) .x(0) = S lit ei x (ei* .x(0)) = S ci lit ei ).        
          We shall now consider as an example the model
										x1(t+1) = x1(t) + 2 x2(t)
										x2(t+1) = 3 x1(t) + 2 x2(t)                                                                                   [8]
with the initial values x1(0) = x2(0) = 1. These equations can be solved by iterating fron one generation to the next recursively without using any of the above formulae:

x[t_]:=x[t]=A.x[t-1]
x[0]={1,1};
A={{1,2},{3,2}};

Table[x[t],{t,0,10}]   //MatrixForm

They can also be solved by using Eqn [5]:

Table[MatrixPower[A,t].x[0],{t,0,10}]    //MatrixForm

Finally, they can be solved using Eqn [7]:

{vals,vecs}=Eigensystem[A]

{lvals,lvecs}=Eigensystem[Transpose[A]]

Do[lvecs[[i]]=lvecs[[i]]/(vecs[[i]].lvecs[[i]]),
					{i,2}]

c=lvecs.x[0]

lvecs

xform[t_]:=N[Sum[c[[i]] vals[[i]]^t vecs[[i]],{i,2}]]

MatrixForm[xform[t]]

Verify that the formula given by xform generates the same solutions as found above

The advantage of this approach is that it shows immediately what the asymptotic behavior of the system will be. In this example, the term in 4t will eventually dominate the term in (–1)t, and the system will asymptotically lie in the direction of the dominant right eigenvector {2/3, 1}, increasing fourfold each year (since the dominant eigenvalue is 4). The constant c multiplying the dominant term (in this case c2 = 1.2) is the inner product of the dominant left eigenvector and the vector of initial values. Thus the elements of this eigenvector reflect the relative importance of the different initial values (which might, for example, represent age classes) in determining c which decides the absolute population sizes. In the present case this eigenvector is {0.6, 0.6}, which means that the two initial values are equally important in determining the ultimate population sizes.
          We now consider an example in which the eigenvalues are complex:
          x1(t+1) = x1(t)+ 2x2(t)
										x2(t+1) = –2x1(t) +  x2(t)                                                                                   [9]
with initial values x1(0) = x2(0) = 1.

A={{1,2},{-2,1}};
x[0]={1,1};

{vals,vecs}=Eigensystem[A]

{lvals,lvecs}=Eigensystem[Transpose[A]]

Do[lvecs[[i]]=lvecs[[i]]/(vecs[[i]].lvecs[[i]]),
					{i,2}]

c=lvecs.x[0]

Since the eigenvalues are complex, it is simpler to express them in polar coordinates, r and q, so that 
          lt = rt (cos q t + i sin q t)

r=Abs[vals]  //N

theta=Arg[vals]  //N

func[i_,t_]:=r[[i]]^t (Cos[theta[[i]] t] +
					I Sin[theta[[i]] t])

x[t_]:=Sum[c[[i]] func[i,t] vecs[[i]],{i,2}];

x[t]=Chop[Simplify[x[t]]]

Verify that this solution agrees with that found by direct iteration of [4]. Plot the solution. Quantify the period of the oscillations. Note that imaginary numbers appear in the intermediate calculations but not in the final answer.

This method can be used for sets of equations with any number of variables. The asymptotic behavior of the system depends on the eigenvalues of the matrix A and can be determined as follows. (We suppose that these eigenvalues are distinct.) Rank the eigenvalues by their absolute value. There are two common situations. (i) There may be a real eigenvalue l greater in absolute value than any of the others. The system will increase or decrease exponentially according as |l| >1 or |l| <1; it will alternate in sign or not according as l is negative or positive. (ii) There may be a pair of complex conjugate eigenvalues greater in absolute value than any of the others. The system will oscillate with increasing or decreasing amplitude according as  |l| >1 or |l| <1.

Local stability analysis of nonlinear difference equations

Consider the difference equation
            x(t+1) = f(x(t)).
 At equilibrium x(t+1) = x(t), so that f(x) = x. Suppose that x*  is an equilibrium, and write x(t) = x* + p(t). Then
           x(t+1) = f(x(t)) = f(x* + p(t)) ~ f(x*) + f '(x*) p(t) = x* + f '(x*) p(t)
whence to a first order approximation
          p(t+1)  = f '(x*) p(t).
The solution of this equation is
          p(t) = [f '(x*)]t p(0)
which converges to zero or diverges from it (giving stability or instability) according as the absolute value of f '(x*) is less than or greater than unity.
          Consider now the set of n difference equations in n unknowns
          x1(t+1) = f1(x1(t), x2(t), ..., xn(t))
          x2(t+1) = f2(x1(t), x2(t), ..., xn(t))
          .
          .
          xn(t+1 )= fn(x1(t), x2(t), ..., xn(t))
An equilibrium satisfies 
           fi(x) = x, i = 1, 2, ..., n.
Thus near an equilibrium
          xi(t+1) = fi(x(t)) = fi(x* + p(t)) ~ f(x*) + Sj (dfi(x*)/dxj) pj. = x* + Sj (dfi(x*)/dxj) pj
so that to a first approximation
          p(t+1) = J p(t)
where J is the Jacobian matrix of derivatives. The ultimate behavior of the solution to this equation depends on the eigenvalue of J with the largest absolute value. The solution converges to zero (giving local stability) or diverges from it (giving instability) according as this absolute value is less than or greater than unity. If this eigenvalue is real, the perturbations will decrease or increase geometrically (with alternating sign if it is negative). If this eigenvalue is complex with argument q (which implies that it has a complex conjugate with the same absolute value but with argument –q), the perturbations will oscillate with period 2p/q (being damped or explosive according as the absolute value is less than or greater than unity).

K  The Poisson distribution

The Poisson distribution was first derived by S.D.Poisson in 1837 and has many applications in biology. It is a discrete probability distribution with the formula
																				P(x) = e-mmx/x!,            x = 0, 1, 2, ...
P(x) represents the probability that some random event occurs exactly x times. It depends on the parameter m which is the mean of the distribution, representing the average number of occurences of the event; a characteristic of the Poisson distribution is that its mean is also its variance.

Write a function to generate the Poisson probabilities from x = 0 through 50 given the parameter m. Implement it with m = 0.5, 1, 2, 5, 20. Display the results, and in each case find the sum of the probabilities and the mean and variance, and verify that they agree with their theoretical values.

The Poisson distribution is a limiting case of the binomial distribution which gives the probability of getting x "successes" and n – x "failures" in n trials if each trial results in one of two outcomes with a constant probability p of success at each trial, which is independent of the outcomes of previous trials. For example, the number of boys in families of fixed size is to a good approximation binomial because the probability of a male birth is very nearly the same for all families and independent of birth order. The formula for this distribution is
          										P(x) = (n!/x!(n – x)!) px(1 – p)n–x,          x = 0, 1, 2, ... , n
The reason is that the probability of getting x successes and n – x failures in a particular order is px(1 – p)n–x and there are n!/x!(n – x)! ways in which this can happen. 
										The Poisson distribution is the limiting case of this distribution when the probability p becomes very small but the number of trials n becomes very large. Write m = np for the mean of the binomial distribution (this is clearly its mean since p is the probability of success in each trial, whence np is the Expected number of successes in n trials), and let n increase and p decrease with m held constant. Then with m constant and n->infinity,
																				n!/x!(n – x)!  -> nx/x!
																				px = mx/nx
																				(1 – p)n–x = (1 – m/n)n–x -> e-m
Hence
																				P(x) -> 		e-mmx/x!,
the Poisson distribution with parameter m.
										A classic example of the Poisson distribution is the numbers of deaths from horse kicks in the Prussian army between 1875 and 1894. Data on 14 army corps were collected by von Bortkewiecz for each of these years, giving 280 observations in all on the number of deaths from this cause in a corps in a year. The data are shown below:
										number of deaths        0      1      2      3      4      5+      Total
										Frequency                 144   91   32     11      2      0         280
One would expect them to follow a Poisson distribution since there is a very large number of men at risk each of whom has a very small chance of being killed by a horse in a year; otherwise the conditions for the binomial distribution hold. The Poisson distribution can be fitted to these data by equating the observed mean number of deaths with the theoretical mean:

data={144,91,32,11,2};
mean = Sum[(i-1)*data[[i]],{i,5}]/280  //N

expected=Table[280 Exp[-mean] mean^x/x!,{x,0,10}]

The Poisson distribution can be derived by an independent argument as the distribution of the number of random events occurring in a fixed time. Suppose that in a short time interval dt there is a chance m dt that an event occurs, which does not depend on the number of events which have already occurred. (This is what we mean by saying that events occur at random.) Write P(n, t) for the probability that exactly n events have occurred by time t, starting with zero events at time zero. Now n events can occur by time (t + dt) in two ways: (1) n events occur by time t and none between t and t + dt, (2) n – 1 events occur by time t and another event occurs between t and t + dt. (The second possibility only applies when n >= 1.) Thus
          P(n, t + dt) = P(n, t) (1 – m dt) + P(n – 1, t) m dt
whence
          dP(n, t)/dt = m (P(n – 1, t) – P(n, t))                 for n>=1
          dP(0, t)/dt = –m P(0, t),
with initial conditions
          P(n, 0) = 0                                                     for n>=1
          P(0, 0) = 1.
It is straightforward to verify that the solution to these differential equations is the Poisson distribution with parameter mt:

p[n_,t_]:=Exp[-mu*t] (l t)^n/n!

p[n,0]

Limit[p[0,t],t->0]

D[p[0,t],t] + mu*p[0,t]

Simplify[D[p[n,t],t] - mu*(p[n-1,t]-p[n,t])]

%/.n/n!->1/(n-1)!

The properties of the Poisson distribution are most easily investigated by calculating its probability generating function (pgf), G(s). The pgf of a discrete random variable X with probability function
          Prob(X = x) = P(x)
is defined as
          G(s) = E(sX) = Sx sxP(x).
For the Poisson distribution
										G(s) =  Sx e–m(sm)x/x! = e–mesm = e(s–1)m
          If we differentiate the pgf with respect to s and then set s = 1, we obtain the expected value of the random variable (the mean of the distribution):
          G '(s) at s=1 = Sx xsx–1P(x) at s=1 = Sx xP(x) = E(X).
Taking the second derivative at s = 1 gives
										G ''(s) at s=1 = E[X(X – 1)] = E(X2) – E(X)
from which the variance can be calculated as
          Var(X) = E(X2) - E2(X) = E[X(X – 1)] + E(X) - E2(X).
For the Poisson distribution we can calculate:

G=Exp[(s-1)*mu];

mean=D[G,s] /.s->1

variance=D[G,{s,2}] + mean - mean^2 /.s->1

This verifies that the mean and the variance of the distribution are both equal to l.
										The rth derivative  of the pgf evaluated at s = 0 is equal to r! P(r). For example consider a variable which can take values 0, 1, ..., 5:

G[s_]:=P[0]+s*P[1]+s^2*P[2]+s^3*P[3]+s^4*P[4]+s^5*P[5]

Table[D[G[s],{s,r}]/r!,{r,0,5}]   /.s->0

This enables the probability distribution to be calculated from the pgf, which is useful if the pgf has been obtained indirectly without direct knowledge of the probability distribution.
          Another useful property is that the pgf of the sum of two independent random variables is equal to the product of their pgf's. Suppose X and Y are independent random variables with pgf's GX(s) and GY(s); the pgf of their sum, X + Y, is
           GX+Y(s) = E(sX+Y) = E(sXsY) (under independence) = GX(s) GY(s)
If X and Y are independent Poisson variates with parameters m1 and  m2 respectively, the pgf of their sum is
          e(s–1)m1 e(s–1)m2 = e(s–1)(m1+ m2)
which is the pgf of a Poisson distribution with parameter m1+ m2. Hence the sum of two independent Poisson variates is itself a Poisson variate, with mean equal to the sum of the means.
          The sum of n identically and independently distributed random variables Xi with pgf GX(s), say Y = X1 + X2 + ... Xn, has pgf GY(s) = (GX(s))n. Suppose however that the number of Xi's is not fixed but is itself a random variable N with probability distribution PN(n) and pgf GN(s). The pgf of Y = X1 + X2 + ... XN is
          GY(s) = PN(0) + PN(1) GX(s) + PN(2) [GX(s)]2 + PN(3) [GX(s)]3 + ...
                    = GN(GX(s))
For example, suppose an individual has a Poisson number of offspring with mean m and then dies, and that each of these offspring has a Poisson number of offspring with the same mean. Then the number of grandoffspring of the original individual has pgf G(G(s)), where G(s) is the pgf of the Poissson distributioin; this can be used to find the probability distribution of the number of grandoffspring. Suppose for example that m = 1; we compare the distribution of the number of offspring (Poisson with mean 1) and the number of grandoffspring:

G[s_]:=Exp[(s-1)]

Table[Exp[-1]/r!,{r,0,10}]  //N

Table[D[G[G[s]],{s,r}]/r!,{r,0,10}]   /.s->0  //N

This process is called a Poisson branching process. If G(s) is the Poisson pgf, the pgf of the descendants of a single individual after t generations is
										Gt(s) = G(G(G(...G(s)...)))
where there are t G's on the right hand side; this can be calculated from the recurrence relationship
											Gt(s) = G(	Gt–1(s)) 
In particular if q(t) is the probability that the line has gone extinct by t years, which is the probability that there are no descendants, then q(t) = Gt(0), which satisfies the recurrence
										q(t) = G(q(t – 1)).
If G(s) is the pgf of a Poisson distribution with mean 1, this becomes
										q(t) = eq(t–1)–1
with the initial condition that q(0) = 0.

q[0]=0.0;
q[t_]:=q[t]=Exp[q[t-1]-1]

Table[{t,q[t]},{t,10,100,10}]  //TableForm