(*^ ::[ Information = "This is a Mathematica Notebook file. It contains ASCII text, and can be transferred by email, ftp, or other text-file transfer utility. It should be read or edited using a copy of Mathematica or MathReader. If you received this as email, use your mail application or copy/paste to save everything from the line containing (*^ down to the line containing ^*) into a plain text file. On some systems you may have to give the file a name ending with ".ma" to allow Mathematica to recognize it as a Notebook. The line below identifies what version of Mathematica created this file, but it can be opened using any other version as well."; FrontEndVersion = "Macintosh Mathematica Notebook Front End Version 2.2"; MacintoshStandardFontEncoding; fontset = title, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e8, 24, "Times"; fontset = subtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6, 18, "Times"; fontset = subsubtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, e6, 14, "Times"; fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, a20, 18, "Times"; fontset = subsection, inactive, pageBreak, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, a15, 14, "Times"; fontset = subsubsection, inactive, noPageBreak, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, a12, 12, "Times"; fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Times"; fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, L-5, 12, "Courier"; fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5, 12, "Courier"; fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R65535, L-5, 12, "Courier"; fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5, 12, "Courier"; fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, B65535, L-5, 12, "Courier"; fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287, 12, "Courier"; fontset = name, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic, 10, "Geneva"; fontset = header, inactive, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = leftheader, inactive, L2, 12, "Times"; fontset = footer, inactive, noKeepOnOnePage, preserveAspect, center, M7, 12, "Times"; fontset = leftfooter, inactive, L2, 12, "Times"; fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Times"; fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, cellOutline, M16, italic, 12, "Times"; fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; paletteColors = 128; showRuler; automaticGrouping; currentKernel; ] :[font = subtitle; inactive; preserveAspect; startGroup] APPENDIX SOME MATHEMATICAL TOOLS FOR POPULATION BIOLOGISTS :[font = subsection; inactive; Cclosed; dontPageBreak; noPageBreak; preserveAspect; startGroup] Introduction :[font = text; inactive; noPageBreak; preserveAspect; endGroup] 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. ;[s] 5:0,0;252,1;263,0;345,1;356,0;479,-1; 2:3,13,9,Times,0,12,0,0,0;2,14,10,Times,2,12,0,0,0; :[font = subsection; inactive; Cclosed; pageBreak; preserveAspect; startGroup] A Taylor series expansions :[font = text; inactive; preserveAspect; leftWrapOffset = 3; leftNameWrapOffset = 4; rightWrapOffset = 469] 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: ;[s] 48:0,0;13,1;18,0;35,1;36,0;87,1;88,2;89,1;90,0;148,1;158,2;159,1;169,2;170,1;173,3;174,1;177,2;178,1;260,0;263,1;264,0;270,1;271,3;272,1;275,2;276,1;278,0;299,1;304,0;317,1;322,2;323,0;367,1;369,0;391,1;392,2;393,0;487,1;488,2;489,0;610,1;621,0;638,4;655,5;656,4;661,0;670,4;671,0;835,-1; 6:15,13,9,Times,0,12,0,0,0;19,14,10,Times,2,12,0,0,0;8,18,11,Times,34,10,0,0,0;2,7,5,Times,2,6,0,0,0;3,13,9,Times,1,12,0,0,0;1,18,12,Times,33,10,0,0,0; :[font = input; preserveAspect] Series[Log[x],{x,2.,1}] :[font = input; preserveAspect] f[x_]=Normal[Series[Log[x],{x,2.,1}]] :[font = input; preserveAspect] Table[{x,Log[x],f[x]},{x,1.5,2.5,.1}] //TableForm :[font = text; inactive; preserveAspect] 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 ;[s] 65:0,0;17,1;24,0;79,5;80,0;140,1;146,0;233,2;234,6;235,0;268,2;270,0;292,2;293,6;294,2;295,0;437,2;441,6;442,2;445,0;503,2;511,7;512,2;516,7;517,3;518,2;522,6;523,2;528,7;529,2;534,6;535,2;536,5;537,2;542,7;543,2;548,6;549,2;550,5;551,4;552,2;558,0;611,2;616,0;619,2;622,6;623,0;638,2;639,8;640,2;645,6;648,2;650,6;651,2;654,0;663,2;665,6;668,2;670,6;671,2;673,0;685,2;686,0;703,2;707,0;734,-1; 9:15,13,9,Times,0,12,0,0,0;2,13,9,Times,1,12,0,0,0;27,14,10,Times,2,12,0,0,0;1,22,14,Times,66,12,0,0,0;1,22,14,Times,34,12,0,0,0;3,18,12,Times,32,10,0,0,0;11,18,11,Times,34,10,0,0,0;4,18,12,Times,64,10,0,0,0;1,18,11,Times,66,10,0,0,0; :[font = text; inactive; preserveAspect; leftNameWrapOffset = 414] 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}]. ;[s] 47:0,0;10,1;20,6;21,1;29,6;30,1;32,10;33,1;36,6;37,1;41,3;45,6;46,3;47,7;48,8;49,1;51,10;52,1;56,6;57,1;61,3;65,6;66,3;67,7;68,1;70,10;71,1;76,6;77,0;84,2;102,0;105,2;164,5;172,2;193,5;195,2;196,0;368,1;372,6;373,1;374,6;375,1;376,0;398,4;414,9;415,4;419,0;421,-1; 11:6,13,9,Times,0,12,0,0,0;13,14,10,Times,2,12,0,0,0;4,22,14,Times,34,12,0,0,0;4,14,10,Times,6,12,0,0,0;2,13,9,Times,1,12,0,0,0;2,21,13,Times,32,12,0,0,0;9,18,11,Times,34,10,0,0,0;2,18,12,Times,32,10,0,0,0;1,12,8,Times,2,10,0,0,0;1,18,12,Times,33,10,0,0,0;3,7,5,Times,2,6,0,0,0; :[font = special1; inactive; preserveAspect] 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. :[font = text; inactive; preserveAspect] Three important Taylor series expansions about x* = 0 are those for the exponential, sine and cosine functions. The divisors are the appropriate factorials: ;[s] 3:0,0;58,1;59,0;167,-1; 2:2,13,9,Times,0,12,0,0,0;1,18,12,Times,32,10,0,0,0; :[font = input; preserveAspect] Series[Exp[x],{x,0,5}] :[font = input; preserveAspect] Series[Sin[x],{x,0,5}] :[font = input; preserveAspect] Series[Cos[x],{x,0,5}] :[font = text; inactive; preserveAspect] 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: ;[s] 64:0,0;77,1;79,2;80,3;81,0;89,2;91,1;95,6;96,1;99,6;100,1;108,7;109,1;110,0;137,2;142,10;143,0;188,1;190,2;191,1;197,2;198,11;199,1;203,9;206,1;207,7;208,1;212,7;213,11;214,1;216,9;217,1;219,9;220,1;221,7;224,4;296,0;299,4;300,0;306,9;307,1;309,9;310,1;311,7;312,4;313,0;342,1;344,0;360,1;361,7;362,8;363,0;376,2;381,11;382,0;437,6;438,0;440,6;441,0;448,6;449,0;456,6;457,5;458,0;463,-1; 12:14,13,9,Times,0,12,0,0,0;17,14,10,Times,2,12,0,0,0;6,14,10,Times,3,12,0,0,0;1,13,9,Times,1,12,0,0,0;3,22,14,Times,66,12,0,0,0;1,21,13,Times,64,12,0,0,0;6,18,12,Times,64,10,0,0,0;6,18,11,Times,66,10,0,0,0;1,12,8,Times,2,10,0,0,0;5,14,10,Symbol,0,12,0,0,0;1,18,11,Times,35,10,0,0,0;3,18,11,Times,34,10,0,0,0; :[font = input; preserveAspect] f[x1_,x2_]=Normal[Series[Log[x1/x2], {x1,2.,1},{x2,1.,1}]] :[font = input; preserveAspect] f[1.9,1.1] :[font = input; preserveAspect; endGroup] Log[1.9/1.1] :[font = subsection; inactive; Cclosed; pageBreak; preserveAspect; startGroup] B Complex numbers :[font = text; inactive; preserveAspect] 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 ;[s] 45:0,0;251,1;253,3;255,1;268,0;306,1;310,0;311,1;318,3;320,1;322,0;323,1;326,3;328,0;329,1;330,0;331,1;332,0;337,1;341,0;342,1;349,3;350,1;353,0;354,1;357,3;359,0;360,1;361,0;362,1;364,0;369,1;370,3;372,1;374,0;375,1;377,0;627,1;629,0;722,1;733,0;734,1;736,0;749,2;751,0;835,-1; 4:18,13,9,Times,0,12,0,0,0;20,14,10,Times,2,12,0,0,0;1,13,9,Times,1,12,0,0,0;6,18,12,Times,32,10,0,0,0; :[font = input; preserveAspect] Solve[x^2 - 2x + 5 ==0] :[font = text; inactive; preserveAspect] 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 ;[s] 31:0,0;19,1;21,0;82,1;90,0;217,1;218,4;219,2;220,0;287,1;288,0;297,1;298,0;311,1;312,0;318,1;320,0;337,1;338,0;344,1;345,0;352,1;354,0;358,1;359,0;362,1;363,4;364,3;365,0;372,1;373,0;374,-1; 5:14,13,9,Times,0,12,0,0,0;13,14,10,Times,2,12,0,0,0;1,22,14,Times,34,12,0,0,0;1,21,13,Times,32,12,0,0,0;2,18,12,Times,32,10,0,0,0; :[font = input; preserveAspect] z1 = 1 + 2 I; z2 = 2 - I; :[font = input; preserveAspect] z1 + z2 :[font = input; preserveAspect] z1 z2 :[font = text; inactive; preserveAspect] 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 ;[s] 60:0,0;9,1;11,2;12,0;15,1;16,2;17,0;18,1;19,2;20,0;45,1;46,2;47,1;49,2;50,0;53,1;54,2;55,0;81,1;83,2;84,0;87,1;88,0;91,1;93,0;108,1;110,0;112,1;113,0;115,1;118,0;124,1;125,0;128,1;129,0;134,1;135,0;139,1;140,0;141,1;142,0;150,1;151,0;163,1;164,0;167,1;168,0;175,1;176,0;180,1;181,0;196,1;197,0;203,1;204,0;219,1;220,2;221,0;223,1;225,0;226,-1; 3:26,13,9,Times,0,12,0,0,0;26,14,10,Times,2,12,0,0,0;8,18,12,Times,64,10,0,0,0; :[font = input; preserveAspect] z3=z1/z2 :[font = input; preserveAspect] z2 z3 :[font = text; inactive; preserveAspect] 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): ;[s] 7:0,0;383,1;391,0;394,1;409,0;449,1;458,0;538,-1; 2:4,13,9,Times,0,12,0,0,0;3,14,10,Times,2,12,0,0,0; :[font = input; preserveAspect] r1=Abs[z1] //N :[font = input; preserveAspect] theta1=Arg[z1] //N :[font = input; preserveAspect] r2=Abs[z2] //N :[font = input; preserveAspect] theta2=Arg[z2] //N :[font = text; inactive; preserveAspect] 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) ;[s] 43:0,0;53,2;54,0;57,2;58,0;61,2;63,0;69,2;70,0;77,2;78,0;81,1;83,0;89,2;90,0;118,2;119,0;122,2;123,0;127,2;128,1;129,0;149,2;150,0;152,2;154,0;159,1;179,3;180,2;182,0;183,2;185,0;187,2;190,0;192,2;194,0;199,2;200,1;202,0;204,2;206,0;210,1;211,0;213,-1; 4:20,13,9,Times,0,12,0,0,0;5,14,10,Symbol,0,12,0,0,0;17,14,10,Times,2,12,0,0,0;1,14,10,Symbol,2,12,0,0,0; :[font = input; preserveAspect] r1 Cos[theta1] :[font = input; preserveAspect] r1 Sin[theta1] :[font = text; inactive; preserveAspect] 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! + ... ;[s] 11:0,0;247,1;248,0;256,1;257,0;260,1;261,2;262,0;268,1;269,2;270,0;280,-1; 3:5,13,9,Times,0,12,0,0,0;4,14,10,Times,2,12,0,0,0;2,18,12,Times,32,10,0,0,0; :[font = input; preserveAspect] Exp[z1]//N :[font = special1; inactive; preserveAspect] Use the following function to evaluate the series expansion for exp(z1) for different values of n. :[font = input; preserveAspect] approxexp[z_,n_]:=Sum[z^i/i!,{i,0,n}]//N :[font = text; inactive; preserveAspect] 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, ;[s] 19:0,0;3,1;4,2;5,0;9,1;11,2;12,0;109,1;110,2;111,0;113,1;115,2;116,0;124,1;125,2;126,0;132,1;133,2;134,0;191,-1; 3:7,13,9,Times,0,12,0,0,0;6,14,10,Times,2,12,0,0,0;6,18,12,Times,64,10,0,0,0; :[font = input; preserveAspect] Exp[z1+z2] //N :[font = input; preserveAspect] Exp[z1] Exp[z2] //N :[font = text; inactive; preserveAspect] 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 ;[s] 61:0,0;18,2;19,0;22,2;24,0;26,2;28,0;43,2;44,0;52,2;53,0;55,2;58,0;66,2;67,0;73,2;75,0;133,2;135,0;142,2;144,0;147,2;148,3;149,2;150,3;151,0;157,2;158,3;159,2;160,3;161,0;199,2;200,3;201,0;207,2;208,3;209,0;221,2;222,0;223,2;224,0;228,2;229,3;230,0;236,2;237,3;238,0;304,2;306,0;308,2;310,0;315,1;317,0;318,2;320,0;325,1;326,0;330,2;331,0;336,2;337,1;339,-1; 4:25,13,9,Times,0,12,0,0,0;3,14,10,Symbol,0,12,0,0,0;25,14,10,Times,2,12,0,0,0;8,18,12,Times,32,10,0,0,0; :[font = input; preserveAspect] z1 :[font = input; preserveAspect] r1 Exp[I theta1] :[font = text; inactive; preserveAspect] 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. ;[s] 35:0,0;26,1;27,2;28,0;29,1;30,2;31,0;34,1;35,2;36,0;41,1;43,3;44,2;45,0;47,1;48,2;49,0;53,1;56,3;57,2;58,0;62,1;63,2;64,0;65,1;66,2;67,0;71,1;73,0;74,3;75,2;76,0;79,3;80,2;81,0;213,-1; 4:12,13,9,Times,0,12,0,0,0;9,14,10,Times,2,12,0,0,0;10,18,12,Times,64,10,0,0,0;4,14,10,Symbol,0,12,0,0,0; :[font = input; preserveAspect] Abs[z1 z2] :[font = input; preserveAspect] r1 r2 :[font = input; preserveAspect] Arg[z1 z2] //N :[font = input; preserveAspect; endGroup] theta1 + theta2 :[font = subsection; inactive; Cclosed; pageBreak; preserveAspect; startGroup] C Elementary vector operations :[font = text; inactive; preserveAspect] 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: ;[s] 5:0,0;64,1;71,0;184,1;192,0;232,-1; 2:3,13,9,Times,0,12,0,0,0;2,13,9,Times,1,12,0,0,0; :[font = input; preserveAspect] agefem={92,37,22,12,6,3,0,0}; agemale={83,67,31,31,16,7,3,2}; :[font = text; inactive; preserveAspect] 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. ;[s] 11:0,0;44,1;45,0;57,1;58,0;101,1;102,0;215,1;225,0;262,1;271,0;378,-1; 2:6,13,9,Times,0,12,0,0,0;5,14,10,Times,2,12,0,0,0; :[font = text; inactive; preserveAspect] To demonstrate elementary arithmetic operations on vectors, I define two symbolic vectors: :[font = input; preserveAspect] v1={a1,a2,a3}; v2={b1,b2,b3}; :[font = text; inactive; preserveAspect] Addition of two vectors of the same length is defined by adding their corresponding elements: :[font = input; preserveAspect] v1+v2 :[font = text; inactive; preserveAspect] For example, the age distribution of Great Tits of both sexes is :[font = input; preserveAspect] agefem + agemale :[font = text; inactive; preserveAspect] Multiplication by a scalar (an ordinary number) is defined by multiplying each element by that number: :[font = input; preserveAspect] 2 v1 :[font = text; inactive; preserveAspect] 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 (.): ;[s] 5:0,0;212,1;223,0;234,2;235,0;238,-1; 3:3,13,9,Times,0,12,0,0,0;1,14,10,Times,2,12,0,0,0;1,13,9,Times,1,12,0,0,0; :[font = input; preserveAspect] v1.v2 :[font = text; inactive; preserveAspect] Geometrically, the inner product of a vector with itself is the square of its magnitude: :[font = input; preserveAspect] v1.v1 :[font = text; inactive; preserveAspect] 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: ;[s] 34:0,0;21,1;22,3;23,1;24,0;28,1;29,3;30,1;31,0;34,4;36,3;37,0;41,4;43,3;44,0;75,2;76,0;109,1;110,3;111,1;113,3;114,1;115,0;116,4;118,3;119,4;121,3;122,0;127,2;129,0;158,2;160,0;169,2;171,0;243,-1; 5:11,13,9,Times,0,12,0,0,0;7,13,9,Times,1,12,0,0,0;4,14,10,Symbol,0,12,0,0,0;8,18,12,Times,64,10,0,0,0;4,14,10,Times,2,12,0,0,0; :[font = input; preserveAspect] agefem.agemale/Sqrt[agefem.agefem*agemale.agemale]//N :[font = input; preserveAspect] ArcCos[%] :[font = text; inactive; preserveAspect] (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: ;[s] 13:0,0;68,1;75,0;186,2;188,0;200,2;202,0;260,2;262,0;264,2;267,0;276,2;278,0;412,-1; 3:7,13,9,Times,0,12,0,0,0;1,13,9,Times,1,12,0,0,0;5,14,10,Symbol,0,12,0,0,0; :[font = input; preserveAspect] agefem.Range[8]/Sum[agefem[[i]],{i,8}] //N :[font = text; inactive; preserveAspect] 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: ;[s] 39:0,0;43,1;44,3;45,0;50,1;51,3;52,1;53,0;64,2;66,0;70,2;72,0;93,2;95,0;97,2;99,0;105,2;107,0;140,2;141,0;155,1;156,3;157,1;158,0;166,2;167,0;181,1;182,3;183,0;204,1;205,3;206,1;207,0;209,1;210,3;211,1;212,0;231,2;242,0;255,-1; 4:15,13,9,Times,0,12,0,0,0;10,13,9,Times,1,12,0,0,0;8,14,10,Times,2,12,0,0,0;6,18,12,Times,64,10,0,0,0; :[font = input; preserveAspect; endGroup] Outer[Times,v1,v2] :[font = subsection; inactive; Cclosed; pageBreak; preserveAspect; startGroup] D Elementary matrix operations :[font = text; inactive; preserveAspect] 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: ;[s] 7:0,0;203,1;214,0;284,2;291,0;295,2;303,0;318,-1; 3:4,13,9,Times,0,12,0,0,0;1,14,10,Times,2,12,0,0,0;2,13,9,Times,1,12,0,0,0; :[font = input; preserveAspect] age={agefem, agemale} :[font = text; inactive; preserveAspect] This can be displayed in conventional matrix format: :[font = input; preserveAspect] MatrixForm[age] :[font = text; inactive; preserveAspect] 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. ;[s] 3:0,0;5,1;9,0;229,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,1,12,0,0,0; :[font = input; preserveAspect] A1={{1,3,4},{9,-1,6}}; A2={{3,-1,2},{10,3,4}}; Aa={{a11,a12},{a21,a22}}; Ab={{b11,b12},{b21,b22}}; :[font = text; inactive; preserveAspect] Addition of two matrices of the same order is defined by adding their corresponding elements. :[font = input; preserveAspect] Aa //MatrixForm Ab //MatrixForm Aa+Ab //MatrixForm :[font = input; preserveAspect] A1 //MatrixForm A2 //MatrixForm A1+A2 //MatrixForm :[font = text; inactive; preserveAspect] The transpose of a matrix A, usually denoted AT, has the rows and columns transposed: ;[s] 6:0,0;36,1;37,0;55,1;56,2;57,0;96,-1; 3:3,13,9,Times,0,12,0,0,0;2,13,9,Times,1,12,0,0,0;1,18,12,Times,32,10,0,0,0; :[font = input; preserveAspect] A1 :[font = input; preserveAspect] Transpose[A1] :[font = input; preserveAspect] % //MatrixForm :[font = text; inactive; preserveAspect] Multiplication of a matrix by a scalar is defined by multiplying each element by the scalar. :[font = input; preserveAspect] 1.5 A1 //MatrixForm :[font = text; inactive; preserveAspect] 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. ;[s] 29:0,0;83,1;85,0;89,1;91,0;114,2;115,0;118,2;119,0;124,2;125,0;128,2;129,0;158,1;162,0;203,2;204,0;207,2;208,0;216,2;218,0;255,2;256,0;266,1;268,0;276,2;277,0;290,1;291,0;294,-1; 3:15,13,9,Times,0,12,0,0,0;5,13,9,Times,1,12,0,0,0;9,14,10,Times,2,12,0,0,0; :[font = input; preserveAspect] Aa.Ab //MatrixForm :[font = input; preserveAspect] Transpose[A1].A2 //MatrixForm :[font = text; inactive; preserveAspect] 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: ;[s] 15:0,0;41,1;45,0;91,1;93,0;128,1;129,0;134,1;136,0;140,1;142,0;193,1;197,0;201,1;205,0;311,-1; 2:8,13,9,Times,0,12,0,0,0;7,13,9,Times,1,12,0,0,0; :[font = input; preserveAspect] {{1,3},{2,-1}}.{{4,6},{-2,8}} :[font = input; preserveAspect] {{4,6},{-2,8}}.{{1,3},{2,-1}} :[font = input; preserveAspect] :[font = text; inactive; preserveAspect; fontSize = 11] [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: ;[s] 33:0,0;19,1;30,0;244,2;246,0;267,1;269,0;271,1;273,0;277,2;279,0;301,1;303,0;305,2;309,0;331,1;333,0;368,2;370,0;375,2;376,0;381,2;383,0;396,1;397,0;399,2;403,0;425,1;427,0;451,2;453,0;473,2;474,0;477,-1; 3:17,12,8,Times,0,11,0,0,0;7,12,8,Times,2,11,0,0,0;9,12,8,Times,1,11,0,0,0; :[font = input; preserveAspect] A1.{x1,x2,x3} :[font = input; preserveAspect] {x1,x2}.A1 :[font = text; inactive; preserveAspect; fontSize = 11] 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: ;[s] 7:0,0;63,1;75,0;79,1;93,0;229,1;240,0;275,-1; 2:4,12,8,Times,0,11,0,0,0;3,12,8,Times,2,11,0,0,0; :[font = input; preserveAspect] v1={1,2,3}; v2={4,5,6}; rowv={{1,2,3}}; colv={{4},{5},{6}}; :[font = text; inactive; preserveAspect; fontSize = 11] Note that :[font = input; preserveAspect] colv.rowv :[font = text; inactive; preserveAspect; fontSize = 11] is the same as the outer product :[font = input; preserveAspect] Outer[Times,v2,v1] :[font = text; inactive; preserveAspect; fontSize = 11] However :[font = input; preserveAspect] rowv.colv :[font = text; inactive; preserveAspect; fontSize = 11] is a matrix with one row and one column whereas :[font = input; preserveAspect] v1.v2 :[font = text; inactive; preserveAspect; fontSize = 11] is a scalar.] :[font = text; inactive; preserveAspect] 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: :[font = text; inactive; preserveAspect] 3x1 + 4x2 Ð x3 = 4 4x1 + x2 + 12x3 = 12 9x1 Ð 3x2 + x3 = 6 ;[s] 28:0,0;16,1;17,2;18,0;22,1;23,2;24,0;28,1;29,2;30,0;51,1;52,2;53,0;57,1;58,2;59,0;64,1;65,2;66,0;89,1;90,2;91,0;95,1;96,2;97,0;100,1;101,2;102,0;107,-1; 3:10,13,9,Times,0,12,0,0,0;9,14,10,Times,2,12,0,0,0;9,18,12,Times,64,10,0,0,0; :[font = text; inactive; preserveAspect] This can be written in matrix terminology more compactly as :[font = text; inactive; preserveAspect; plain; bold] A.x = c :[font = text; inactive; preserveAspect] 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: ;[s] 7:0,0;6,1;7,0;44,1;46,0;77,1;78,0;136,-1; 2:4,13,9,Times,0,12,0,0,0;3,13,9,Times,1,12,0,0,0; :[font = input; preserveAspect] A = {{3,4,-1},{4,1,12},{9,-3,1}}; c={4,12,6}; :[font = text; inactive; preserveAspect] Note the value of A.x: ;[s] 3:0,0;19,1;22,0;24,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,1,12,0,0,0; :[font = input; preserveAspect] A.{x1,x2,x3} //ColumnForm :[font = text; inactive; preserveAspect] 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 ;[s] 4:0,0;203,1;208,2;210,1;213,-1; 3:1,13,9,Times,0,12,0,0,0;2,13,9,Times,1,12,0,0,0;1,18,12,Times,32,10,0,0,0; :[font = input; preserveAspect] Inverse[A].c //N :[font = text; inactive; preserveAspect; endGroup] This is the solution for x as will be explained in the next section. ;[s] 4:0,1;41,0;66,1;68,0;110,-1; 2:2,13,9,Times,0,12,0,0,0;2,13,9,Times,1,12,0,0,0; :[font = subsection; inactive; Cclosed; pageBreak; preserveAspect; startGroup] E Determinants, linear independence and matrix inversion :[font = text; inactive; preserveAspect] 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. ;[s] 9:0,0;18,1;20,0;44,2;55,0;245,2;253,0;318,2;330,0;332,-1; 3:5,13,9,Times,0,12,0,0,0;1,13,9,Times,1,12,0,0,0;3,14,10,Times,2,12,0,0,0; :[font = input; preserveAspect] A3={{1,2,4},{3,1,2},{7,4,8}}; A4={{1.,2,3},{2.,1,2},{3.,4,5}}; :[font = input; preserveAspect] Det[A3] :[font = input; preserveAspect] Det[A4] :[font = text; inactive; preserveAspect] 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. ;[s] 22:0,1;1,3;2,1;3,0;116,1;117,3;118,1;119,0;331,1;346,0;380,2;381,0;383,2;384,0;391,2;392,0;393,2;394,0;395,2;397,0;406,2;408,0;421,-1; 4:9,13,9,Times,0,12,0,0,0;5,13,9,Times,1,12,0,0,0;6,14,10,Times,2,12,0,0,0;2,18,12,Times,64,10,0,0,0; :[font = input; preserveAspect] Det[Array[a,{2,2}]] :[font = input; preserveAspect] Det[Array[a,{3,3}]] :[font = text; inactive; preserveAspect] 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: ;[s] 38:0,0;148,2;150,0;245,2;246,0;442,1;451,0;848,2;858,0;1702,1;1704,0;1746,2;1747,0;1749,2;1750,0;1783,1;1784,3;1786,4;1787,0;1797,1;1800,3;1802,0;1805,1;1806,3;1808,1;1810,0;1813,1;1814,5;1815,0;1823,1;1824,5;1825,0;1883,2;1884,0;1885,2;1887,0;2028,2;2039,0;2051,-1; 6:16,13,9,Times,0,12,0,0,0;8,13,9,Times,1,12,0,0,0;8,14,10,Times,2,12,0,0,0;3,18,12,Times,32,10,0,0,0;1,12,9,Times,0,10,0,0,0;2,22,14,Times,66,12,0,0,0; :[font = input; preserveAspect] IdentityMatrix[4] //MatrixForm :[font = input; preserveAspect] B=Inverse[A4] :[font = input; preserveAspect] A4.B //MatrixForm :[font = input; preserveAspect] B.A4 :[font = text; inactive; preserveAspect] 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. ;[s] 63:0,0;65,1;72,2;73,0;82,1;83,3;84,0;93,4;94,0;104,1;106,0;110,1;111,3;112,0;121,4;122,0;135,1;136,0;160,1;162,0;179,1;180,5;181,0;184,1;185,5;186,0;190,1;191,5;192,0;199,1;200,5;201,1;203,5;204,0;214,1;215,5;216,0;220,1;221,5;222,0;223,1;225,5;226,0;229,1;230,5;231,1;232,0;233,5;234,0;238,1;239,5;240,1;242,5;243,0;286,1;287,5;288,0;291,1;292,5;293,0;297,1;298,5;299,0;301,-1; 6:22,13,9,Times,0,12,0,0,0;21,13,9,Times,1,12,0,0,0;1,22,14,Times,66,12,0,0,0;2,18,11,Times,66,10,0,0,0;2,14,10,Times,2,12,0,0,0;15,18,12,Times,64,10,0,0,0; :[font = subsubsection; inactive; Cclosed; noPageBreak; preserveAspect; startGroup] Solution of linear equations :[font = text; inactive; preserveAspect] We return to the problem of solving n linear equations in n unknowns: :[font = output; input; inactive; preserveAspect] A.x = c ;[s] 3:0,0;12,1;14,0;16,-1; 2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,0,12,0,0,0; :[font = text; inactive; preserveAspect] 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 ;[s] 17:0,0;6,1;8,0;42,3;43,4;45,0;47,1;49,0;79,1;81,0;173,1;175,2;176,0;180,1;181,0;186,1;188,0;300,-1; 5:8,13,9,Times,0,12,0,0,0;6,13,9,Times,1,12,0,0,0;1,13,9,Times,129,12,0,0,0;1,14,10,Times,2,12,0,0,0;1,18,11,Times,66,10,0,0,0; :[font = output; input; inactive; preserveAspect] x = A-1.c ;[s] 5:0,0;10,1;12,0;13,2;15,0;18,-1; 3:3,12,10,Courier,1,12,0,0,0;1,12,10,Courier,0,12,0,0,0;1,17,12,Courier,32,10,0,0,0; :[font = text; inactive; preserveAspect] 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 ;[s] 9:0,0;13,1;15,0;120,2;121,0;217,1;219,0;248,2;250,0;310,-1; 3:5,13,9,Times,0,12,0,0,0;2,13,9,Times,1,12,0,0,0;2,14,10,Times,2,12,0,0,0; :[font = output; input; inactive; preserveAspect] x1 + x2 = 1 2x1 + 2x2 = 3 ;[s] 12:0,1;1,2;2,0;5,1;6,2;7,0;13,1;14,2;15,0;19,1;20,2;21,0;26,-1; 3:4,12,10,Courier,0,12,0,0,0;4,12,10,Courier,2,12,0,0,0;4,17,12,Courier,64,10,0,0,0; :[font = text; inactive; preserveAspect] have no solution whereas the equations :[font = output; input; inactive; preserveAspect] x1 + x2 = 1 2x1 + 2x2 = 2 ;[s] 12:0,1;1,2;2,0;5,1;6,2;7,0;13,1;14,2;15,0;19,1;20,2;21,0;26,-1; 3:4,12,10,Courier,0,12,0,0,0;4,12,10,Courier,2,12,0,0,0;4,17,12,Courier,64,10,0,0,0; :[font = text; inactive; preserveAspect] are compatible with any values of x1 and x2 on the line x2 = 1 Ð x1. Consider finally the homogeneous equation ;[s] 13:0,0;34,1;35,2;36,0;40,1;42,2;43,0;55,1;57,2;58,0;65,1;66,2;67,0;121,-1; 3:5,13,9,Times,0,12,0,0,0;4,14,10,Times,2,12,0,0,0;4,18,12,Times,64,10,0,0,0; :[font = output; input; inactive; preserveAspect] A.x = 0 ;[s] 4:0,0;7,1;11,0;13,1;15,-1; 2:2,12,10,Courier,0,12,0,0,0;2,12,10,Courier,1,12,0,0,0; :[font = text; inactive; preserveAspect] 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 ;[s] 11:0,0;3,1;5,0;108,1;110,0;112,1;113,0;118,1;120,0;180,1;182,0;301,-1; 2:6,13,9,Times,0,12,0,0,0;5,13,9,Times,1,12,0,0,0; :[font = output; input; inactive; preserveAspect] x1 + x2 = 0 2x1 + 2x2 = 0 ;[s] 12:0,1;1,2;2,0;5,1;6,2;7,0;13,1;14,2;15,0;19,1;20,2;21,0;26,-1; 3:4,12,10,Courier,0,12,0,0,0;4,12,10,Courier,2,12,0,0,0;4,17,12,Courier,64,10,0,0,0; :[font = text; inactive; preserveAspect; endGroup; endGroup] are satisfied by any values on the line x2 = Ðx1. ;[s] 7:0,0;40,1;41,2;42,0;46,1;47,2;48,0;50,-1; 3:3,13,9,Times,0,12,0,0,0;2,14,10,Times,2,12,0,0,0;2,18,12,Times,64,10,0,0,0; :[font = subsection; inactive; Cclosed; pageBreak; preserveAspect; startGroup] F Eigenvalues and eigenvectors :[font = text; inactive; preserveAspect] 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: ;[s] 21:0,0;12,1;14,0;32,1;33,0;47,2;48,0;51,2;52,0;57,1;59,0;81,2;82,0;89,1;93,0;121,2;122,0;129,1;131,0;196,2;197,0;270,-1; 3:11,13,9,Times,0,12,0,0,0;5,13,9,Times,1,12,0,0,0;5,14,10,Times,2,12,0,0,0; :[font = input; preserveAspect] A = {{1,2},{3,2}}; :[font = input; preserveAspect] A.{0,1} :[font = input; preserveAspect] A.{1,1} :[font = text; inactive; preserveAspect] 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: ;[s] 63:0,0;2,1;191,0;192,1;193,3;194,1;200,0;203,1;204,3;205,1;224,0;225,1;226,3;227,1;297,3;298,1;304,0;306,1;493,2;495,1;511,0;513,1;551,0;555,1;557,2;558,0;641,1;644,0;645,1;650,2;652,1;665,3;676,1;679,0;681,1;685,0;687,1;708,3;719,1;761,0;763,1;765,2;766,0;767,4;768,1;769,0;772,1;774,0;775,1;956,0;958,1;960,2;961,0;962,4;963,1;1119,3;1120,1;1124,2;1126,1;1164,3;1165,1;1263,3;1264,1;1394,-1; 5:17,13,9,Times,1,12,0,0,0;29,13,9,Times,0,12,0,0,0;6,14,10,Symbol,0,12,0,0,0;9,14,10,Times,2,12,0,0,0;2,18,11,Times,66,10,0,0,0; :[font = input; preserveAspect] Eigenvalues[A] :[font = text; inactive; preserveAspect] 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: ;[s] 11:0,0;62,1;63,0;151,1;163,0;278,1;290,0;396,1;401,0;433,1;438,0;467,-1; 2:6,13,9,Times,0,12,0,0,0;5,13,9,Times,1,12,0,0,0; :[font = input; preserveAspect] {vals,vecs}=Eigensystem[A]; :[font = input; preserveAspect] vals :[font = input; preserveAspect] vecs :[font = text; inactive; preserveAspect] 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: ;[s] 3:0,0;359,1;361,0;425,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,1,12,0,0,0; :[font = input; preserveAspect] A.{2,3} :[font = input; preserveAspect] A.{20,30} :[font = text; inactive; preserveAspect] 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: ;[s] 48:0,0;48,1;58,0;61,1;62,0;132,2;134,0;144,1;155,0;186,3;187,0;234,2;235,0;236,1;237,0;243,2;246,0;247,1;248,0;276,2;278,0;422,3;423,0;461,2;463,0;503,3;504,0;628,1;634,0;690,1;695,0;740,2;744,0;746,3;747,2;748,0;1025,1;1036,0;1176,2;1179,0;1180,4;1181,0;1184,2;1185,4;1186,2;1188,0;1447,2;1448,0;1450,-1; 5:23,13,9,Times,0,12,0,0,0;8,14,10,Times,2,12,0,0,0;11,13,9,Times,1,12,0,0,0;4,14,10,Symbol,0,12,0,0,0;2,18,12,Times,32,10,0,0,0; :[font = input; preserveAspect] {lvals, lvecs}=Eigensystem[Transpose[A]]; :[font = input; preserveAspect] lvals :[font = input; preserveAspect] Do[lvecs[[i]]=lvecs[[i]]/(lvecs[[i]].vecs[[i]]), {i,2}] :[font = input; preserveAspect] lvecs :[font = text; inactive; preserveAspect] 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. ;[s] 5:0,0;76,1;81,0;85,1;90,0;199,-1; 2:3,13,9,Times,0,12,0,0,0;2,13,9,Times,1,12,0,0,0; :[font = text; inactive; preserveAspect] 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: ;[s] 120:0,0;134,1;135,4;136,1;137,0;177,3;178,5;179,3;180,0;189,1;190,4;191,6;192,1;193,0;232,3;233,5;234,3;236,0;240,3;251,1;252,4;253,6;254,1;258,4;259,1;260,0;262,1;263,4;264,6;265,1;266,0;267,3;268,5;269,1;270,4;271,0;276,3;277,4;278,1;279,4;280,6;281,0;282,1;284,4;285,1;286,0;302,1;303,4;304,6;305,1;307,4;308,1;309,0;316,3;317,5;318,3;319,0;331,3;332,0;333,3;334,4;335,3;338,0;389,1;390,7;391,6;392,1;394,7;395,2;396,0;413,1;415,0;418,1;419,4;420,0;422,1;423,4;424,0;432,1;433,7;434,0;476,1;478,0;481,1;482,4;483,6;484,0;486,1;487,4;488,6;489,0;497,1;498,7;499,6;500,0;501,1;502,0;568,1;570,0;588,1;589,6;590,0;602,1;605,6;606,1;607,0;609,1;610,7;611,2;723,0;779,1;780,6;781,1;784,0;787,1;788,7;789,2;900,0;935,1;936,0;938,-1; 8:31,13,9,Times,0,12,0,0,0;36,13,9,Times,1,12,0,0,0;3,21,13,Times,64,12,0,0,0;12,14,10,Symbol,0,12,0,0,0;16,18,12,Times,64,10,0,0,0;4,18,12,Symbol,64,10,0,0,0;12,18,12,Times,32,10,0,0,0;6,18,11,Times,66,10,0,0,0; :[font = input; preserveAspect] lvecs.Transpose[vecs] :[font = input; preserveAspect] Transpose[vecs].lvecs :[font = text; inactive; preserveAspect; endGroup] 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. ;[s] 133:0,0;95,1;97,0;114,2;116,0;135,2;136,0;149,1;151,0;155,1;156,6;157,0;179,1;181,0;184,1;185,0;187,1;189,0;201,2;203,0;206,2;207,0;209,2;210,0;215,1;216,0;220,1;221,0;223,1;224,0;246,1;250,0;252,2;253,1;264,0;265,1;266,0;268,1;269,0;272,1;273,0;275,1;276,0;281,2;282,0;284,2;287,1;288,0;290,1;291,0;332,1;336,0;338,1;340,0;341,1;343,0;345,2;346,0;347,1;349,0;350,1;351,0;352,2;353,1;355,0;381,1;382,6;383,0;405,1;406,6;407,1;410,0;412,2;413,1;414,6;415,1;416,0;425,1;426,6;427,1;430,0;445,1;446,6;447,1;448,0;449,1;452,0;456,2;457,0;563,2;564,0;587,1;589,0;607,3;608,4;611,0;615,3;616,4;619,0;622,3;623,4;625,0;643,2;646,7;647,3;648,8;649,6;650,0;651,3;652,8;653,0;655,3;656,8;658,9;659,0;770,3;771,8;772,6;773,0;777,3;779,8;780,5;781,0;789,3;790,0;796,3;798,0;813,1;814,6;815,1;816,0;820,1;822,0;927,-1; 10:52,13,9,Times,0,12,0,0,0;37,13,9,Times,1,12,0,0,0;14,14,10,Symbol,0,12,0,0,0;10,14,10,Times,2,12,0,0,0;3,22,14,Times,66,12,0,0,0;1,21,13,Times,64,12,0,0,0;9,18,12,Times,32,10,0,0,0;1,14,10,Symbol,2,12,0,0,0;5,18,11,Times,66,10,0,0,0;1,12,8,Times,2,10,0,0,0; :[font = subsection; inactive; Cclosed; pageBreak; preserveAspect; startGroup] G Spectral decomposition of a matrix :[font = text; inactive; preserveAspect] 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: ;[s] 32:0,1;6,0;80,3;102,0;137,1;140,4;141,1;142,0;144,1;145,4;146,1;147,2;240,1;249,0;252,2;253,0;259,1;260,0;322,2;324,0;563,3;565,0;598,3;599,0;632,3;633,0;672,3;674,0;706,3;707,0;742,1;744,0;759,-1; 5:13,13,9,Times,0,12,0,0,0;8,13,9,Times,1,12,0,0,0;3,14,10,Symbol,1,12,0,0,0;6,14,10,Times,2,12,0,0,0;2,18,12,Times,32,10,0,0,0; :[font = input; preserveAspect] A.Transpose[vecs] :[font = input; preserveAspect] Lambda = DiagonalMatrix[vals] :[font = input; preserveAspect] Transpose[vecs].Lambda :[font = text; inactive; preserveAspect] Post-multiplying [1] by L and using [F6] we find that A = RT.L.L [2] This equation can be expressed as ;[s] 11:0,0;34,1;36,0;74,1;76,0;78,1;79,3;80,1;81,2;83,1;84,0;229,-1; 4:4,13,9,Times,0,12,0,0,0;5,13,9,Times,1,12,0,0,0;1,14,10,Symbol,1,12,0,0,0;1,18,12,Times,32,10,0,0,0; :[font = text; inactive; preserveAspect] A = S li Ti [3] where Ti = ei x ei*, ;[s] 25:0,0;10,1;12,0;14,2;18,4;19,3;20,1;21,4;183,0;186,4;295,3;298,0;327,3;415,0;442,1;443,4;444,5;445,0;447,1;448,4;449,0;452,1;453,6;454,7;455,0;457,-1; 8:8,13,9,Times,0,12,0,0,0;5,13,9,Times,1,12,0,0,0;1,14,10,Symbol,0,12,0,0,0;3,21,13,Times,64,12,0,0,0;5,18,11,Times,66,10,0,0,0;1,12,8,Times,2,10,0,0,0;1,18,12,Times,64,10,0,0,0;1,18,12,Times,32,10,0,0,0; :[font = text; inactive; preserveAspect] 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: ;[s] 12:0,0;73,1;96,0;100,2;101,0;178,2;180,0;184,2;186,0;302,3;303,4;304,0;306,-1; 5:6,13,9,Times,0,12,0,0,0;1,14,10,Times,2,12,0,0,0;3,13,9,Times,1,12,0,0,0;1,14,10,Symbol,0,12,0,0,0;1,18,12,Times,64,10,0,0,0; :[font = input; preserveAspect] R={{r11,r12,r13},{r21,r22,r23},{r31,r32,r33}} :[font = input; preserveAspect] L={{l11,l12,l13},{l21,l22,l23},{l31,l32,l33}} :[font = input; preserveAspect] Diag=DiagonalMatrix[{1,0,0}] :[font = input; preserveAspect] Transpose[R].Diag.L //MatrixForm :[font = input; preserveAspect] Outer[Times,R[[1]],L[[1]]] //MatrixForm :[font = text; inactive; preserveAspect] 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] ;[s] 99:0,0;80,1;81,5;82,1;84,5;85,0;88,1;89,5;90,2;101,1;102,5;103,1;105,5;106,6;107,0;109,1;110,0;112,3;113,7;114,8;116,7;117,6;119,0;234,1;235,5;236,1;238,5;239,2;240,0;243,1;244,5;245,0;248,1;249,5;250,9;251,0;253,1;254,5;255,0;258,1;259,5;260,9;261,0;267,1;268,5;269,9;270,0;272,1;273,5;274,2;275,0;282,3;283,0;286,3;287,0;315,3;316,0;328,1;330,0;343,1;345,10;346,0;350,4;353,5;354,2;355,1;356,5;357,2;358,0;359,10;360,0;363,2;364,4;368,5;369,10;370,2;371,1;372,5;373,2;440,0;456,2;457,0;528,4;529,11;530,0;681,3;683,0;712,1;713,10;714,1;715,0;717,4;720,11;721,10;722,2;723,1;724,11;727,2;828,0;835,-1; 12:25,13,9,Times,0,12,0,0,0;21,13,9,Times,1,12,0,0,0;11,21,13,Times,64,12,0,0,0;5,14,10,Times,2,12,0,0,0;4,14,10,Symbol,0,12,0,0,0;17,18,11,Times,66,10,0,0,0;2,12,8,Times,2,10,0,0,0;2,12,9,Times,0,10,0,0,0;1,12,9,Times,128,10,0,0,0;3,18,12,Times,32,10,0,0,0;5,18,11,Times,34,10,0,0,0;3,18,12,Times,64,10,0,0,0; :[font = input; preserveAspect] Do[T[i]=Outer[Times,vecs[[i]],lvecs[[i]]],{i,2}] :[font = input; preserveAspect] T[1] :[font = input; preserveAspect] T[1].T[1] :[font = input; preserveAspect] T[1].T[2] :[font = input; preserveAspect] MatrixPower[A,10] :[font = input; preserveAspect] vals[[1]]^10 T[1] + vals[[2]]^10 T[2] :[font = input; preserveAspect] vals[[2]]^10 T[2] //N :[font = text; inactive; preserveAspect] 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] ;[s] 81:0,0;101,1;102,0;223,1;224,0;228,1;230,0;232,1;234,0;236,1;237,6;238,0;244,1;245,6;246,0;252,1;253,6;254,0;359,1;360,0;361,1;362,0;365,2;374,5;375,3;376,0;378,3;379,2;380,5;381,6;382,0;388,2;389,5;390,6;391,0;401,1;402,5;404,3;405,0;407,3;409,0;410,2;412,0;416,2;417,5;418,0;419,3;420,1;421,5;423,7;424,0;524,1;525,4;526,0;534,4;536,0;559,1;560,4;561,0;567,2;569,0;573,2;574,5;575,4;576,0;577,3;578,1;579,5;581,3;582,0;584,2;585,0;589,2;590,8;591,4;592,0;593,3;594,1;595,8;598,3;659,0;663,-1; 9:27,13,9,Times,0,12,0,0,0;15,13,9,Times,1,12,0,0,0;9,14,10,Symbol,0,12,0,0,0;9,21,13,Times,64,12,0,0,0;5,14,10,Times,2,12,0,0,0;8,18,11,Times,66,10,0,0,0;5,18,12,Times,32,10,0,0,0;1,12,8,Times,2,10,0,0,0;2,18,12,Times,64,10,0,0,0; :[font = input; preserveAspect] MatrixExp[3.5 A] //MatrixForm :[font = input; preserveAspect] Exp[3.5 vals[[1]]] T[1] + Exp[3.5 vals[[2]]] T[2] // MatrixForm :[font = input; preserveAspect] Exp[3.5 vals[[2]]] T[2] //MatrixForm :[font = text; inactive; preserveAspect; endGroup] 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 ;[s] 19:0,0;69,1;70,0;92,1;94,0;104,1;105,3;106,4;107,0;109,1;110,0;148,1;150,0;152,2;155,5;156,3;159,0;160,1;161,5;163,-1; 6:7,13,9,Times,0,12,0,0,0;6,13,9,Times,1,12,0,0,0;1,14,10,Symbol,0,12,0,0,0;2,18,12,Times,32,10,0,0,0;1,12,9,Times,1,10,0,0,0;2,18,11,Times,66,10,0,0,0; :[font = subsection; inactive; Cclosed; pageBreak; preserveAspect; startGroup] H Linear differential equations and local stability analysis in continuous time :[font = subsubsection; inactive; Cclosed; noPageBreak; preserveAspect; startGroup] Linear differential equations :[font = text; inactive; preserveAspect] 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] ;[s] 21:0,0;54,1;56,0;57,1;172,0;175,1;176,0;223,1;225,0;248,1;250,0;296,1;298,0;299,1;301,0;341,1;349,0;350,1;352,0;356,1;457,0;461,-1; 2:11,13,9,Times,0,12,0,0,0;10,14,10,Times,2,12,0,0,0; :[font = special1; inactive; preserveAspect] 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. ;[s] 3:0,0;59,1;65,0;158,-1; 2:2,14,10,Times,2,12,0,0,0;1,14,10,Times,3,12,0,0,0; :[font = text; inactive; preserveAspect] 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 ;[s] 112:0,0;34,1;36,0;69,1;71,0;93,1;94,5;95,1;96,0;97,1;102,5;104,2;105,1;106,5;107,1;111,5;113,1;115,5;117,1;126,5;127,6;128,2;129,1;130,6;131,1;142,0;143,1;144,5;145,1;146,0;147,1;152,5;154,2;155,1;156,5;157,2;158,1;161,5;163,2;164,1;165,5;166,2;167,1;176,5;177,6;178,1;180,6;181,7;184,1;326,0;329,1;352,0;353,1;354,6;355,1;356,0;357,1;362,6;363,5;364,2;365,1;366,5;367,2;368,1;371,6;372,5;373,2;374,1;375,5;376,1;386,6;388,1;390,6;391,1;392,0;448,3;449,0;451,1;454,3;455,4;537,0;540,1;541,0;662,4;663,1;664,0;665,1;667,0;681,1;682,4;683,1;689,0;693,3;694,1;696,3;697,4;698,1;699,0;700,1;793,0;796,2;797,0;896,3;897,1;898,0;901,1;902,0;905,3;907,0;911,3;912,1;913,0;977,-1; 8:23,13,9,Times,0,12,0,0,0;42,14,10,Times,2,12,0,0,0;10,22,14,Times,66,12,0,0,0;7,13,9,Times,1,12,0,0,0;4,14,10,Times,3,12,0,0,0;16,18,12,Times,64,10,0,0,0;9,18,11,Times,66,10,0,0,0;1,12,8,Times,2,10,0,0,0; :[font = text; inactive; preserveAspect] ci = ei*.x(0) [6] ;[s] 13:0,0;11,1;12,3;13,1;16,5;17,3;18,4;19,1;20,2;21,1;22,0;23,1;24,0;132,-1; 6:3,13,9,Times,0,12,0,0,0;5,14,10,Times,2,12,0,0,0;1,14,10,Times,3,12,0,0,0;2,18,11,Times,66,10,0,0,0;1,18,11,Times,34,10,0,0,0;1,13,9,Times,1,12,0,0,0; :[font = text; inactive; preserveAspect] we find that :[font = text; inactive; preserveAspect] x(t) = S ci exp(lit) ei [7] ;[s] 17:0,0;10,1;11,0;12,3;13,0;17,2;19,3;20,4;21,0;26,2;27,4;28,3;29,0;31,1;32,4;34,5;35,0;129,-1; 6:6,13,9,Times,0,12,0,0,0;2,13,9,Times,1,12,0,0,0;2,14,10,Symbol,0,12,0,0,0;3,14,10,Times,2,12,0,0,0;3,18,11,Times,66,10,0,0,0;1,12,8,Times,2,10,0,0,0; :[font = text; inactive; preserveAspect] (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 ). ;[s] 60:0,0;7,1;9,0;13,1;14,6;15,5;16,0;17,2;18,3;19,6;20,3;22,0;25,2;26,0;29,1;31,0;35,1;36,6;37,5;38,0;40,3;41,6;42,7;43,0;45,3;46,6;47,8;48,0;50,3;52,0;58,1;60,0;64,1;65,6;66,5;67,0;68,4;69,3;70,6;71,0;75,3;76,6;77,8;78,0;79,3;81,0;176,5;177,1;179,5;180,6;181,1;182,0;186,1;187,6;188,5;189,0;190,4;191,3;192,6;194,0;200,-1; 9:18,13,9,Times,0,12,0,0,0;9,14,10,Symbol,0,12,0,0,0;2,21,13,Times,64,12,0,0,0;9,13,9,Times,1,12,0,0,0;2,21,13,Times,65,12,0,0,0;6,14,10,Times,2,12,0,0,0;11,18,11,Times,66,10,0,0,0;1,12,8,Times,2,10,0,0,0;2,18,12,Times,32,10,0,0,0; :[font = text; inactive; preserveAspect] Example: Solve dx1/dt = x1 + 2x2 dx2/dt = 3x1 + 2x2 with x1(0) = x2(0) = 1. ;[s] 38:0,0;26,1;27,2;28,1;29,0;30,1;35,2;37,1;39,0;40,1;41,2;42,1;53,0;54,1;55,2;56,1;57,0;58,1;62,0;63,1;64,2;65,1;68,0;69,1;70,2;71,1;72,0;77,1;78,2;79,1;80,0;81,1;86,2;87,1;88,0;89,1;93,0;94,1;96,-1; 3:11,13,9,Times,0,12,0,0,0;19,14,10,Times,2,12,0,0,0;8,18,12,Times,64,10,0,0,0; :[font = input; preserveAspect] x1[0]={1,1}; A1={{1,2},{3,2}}; :[font = input; preserveAspect] {vals1,vecs1}=Eigensystem[A1] :[font = input; preserveAspect] {lvals1,lvecs1}=Eigensystem[Transpose[A1]] :[font = input; preserveAspect] Do[lvecs1[[i]]=lvecs1[[i]]/(vecs1[[i]].lvecs1[[i]]), {i,2}] :[font = input; preserveAspect] c1=lvecs1.x1[0] :[font = input; preserveAspect; endGroup] x1[t]=Sum[c1[[i]] Exp[vals1[[i]] t] vecs1[[i]],{i,2}] :[font = subsubsection; inactive; Cclosed; noPageBreak; preserveAspect; startGroup] Local stability analysis of nonlinear differential equations :[font = text; inactive; preserveAspect] 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: ;[s] 111:0,0;56,1;58,0;59,1;67,0;111,1;112,0;113,1;114,0;115,1;118,0;219,1;234,0;314,1;330,0;463,1;464,2;466,0;498,1;501,2;502,1;504,0;540,1;542,0;570,1;575,2;576,1;580,0;588,1;590,0;626,1;631,0;666,1;667,2;668,0;680,1;690,2;691,1;696,3;697,1;700,2;701,1;706,0;721,1;723,0;724,1;728,0;729,1;731,2;732,1;738,0;739,1;741,0;744,1;746,0;747,1;749,0;755,1;756,2;758,0;815,1;817,0;818,1;824,3;825,1;828,2;829,1;833,0;861,1;862,0;863,1;864,0;867,1;869,0;872,1;873,0;876,1;878,3;879,1;882,2;883,1;888,0;953,1;955,3;956,1;959,2;960,1;961,0;981,1;982,0;983,1;984,0;985,1;986,0;1030,1;1031,2;1032,1;1033,0;1115,1;1118,0;1175,1;1176,2;1177,0;1323,1;1324,0;1325,1;1327,0;1328,1;1334,0;1335,1;1344,0;1412,-1; 4:41,13,9,Times,0,12,0,0,0;53,14,10,Times,2,12,0,0,0;13,18,11,Times,34,10,0,0,0;4,7,5,Times,2,6,0,0,0; :[font = input; preserveAspect] f[n_]:=r (1 - n/k) n :[font = input; preserveAspect] equil=Solve[f[n]==0,n] :[font = input; preserveAspect] D[f[n],n] /.equil :[font = text; inactive; preserveAspect] 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: ;[s] 182:0,0;22,1;23,0;25,1;28,0;32,1;34,0;136,1;138,0;164,1;166,0;188,1;189,7;190,0;192,1;197,7;198,1;200,7;201,1;204,7;205,1;213,8;214,1;225,0;227,1;228,7;229,0;231,1;236,7;237,1;239,7;240,1;243,7;244,1;252,8;253,1;289,0;290,1;291,2;292,0;294,1;299,8;300,1;302,7;303,1;306,7;307,1;315,8;316,1;318,0;324,3;326,0;390,1;391,8;392,1;393,4;394,1;403,0;415,1;419,0;432,4;433,9;435,0;509,4;511,0;513,4;514,9;515,4;516,0;518,4;519,0;527,4;529,0;574,1;575,8;576,1;577,3;578,1;582,6;583,8;584,2;585,1;586,6;587,1;588,8;589,1;590,3;591,9;592,1;594,6;595,1;596,8;597,1;600,8;601,1;602,0;698,3;699,1;700,0;701,1;705,3;708,0;714,4;716,0;745,1;747,0;761,6;762,1;763,8;764,1;765,3;766,9;767,1;769,6;770,1;771,8;772,0;867,4;869,0;901,6;902,0;905,1;911,0;916,1;918,0;1105,6;1107,0;1116,1;1118,0;1199,6;1201,0;1342,1;1343,0;1428,1;1429,7;1430,0;1432,1;1438,7;1439,2;1440,1;1442,5;1443,7;1444,1;1445,7;1446,1;1449,7;1450,1;1461,0;1462,1;1463,7;1464,0;1466,1;1473,7;1474,2;1475,1;1477,2;1478,5;1479,7;1480,1;1481,7;1482,1;1485,7;1486,2;1487,0;1536,1;1537,7;1538,1;1541,7;1542,1;1548,7;1549,1;1550,5;1551,7;1552,0;1554,2;1555,1;1556,7;1557,1;1558,5;1559,7;1560,1;1561,0;1618,-1; 10:40,13,9,Times,0,12,0,0,0;66,14,10,Times,2,12,0,0,0;7,22,14,Times,66,12,0,0,0;6,14,10,Times,3,12,0,0,0;9,13,9,Times,1,12,0,0,0;4,14,10,Symbol,2,12,0,0,0;8,14,10,Symbol,0,12,0,0,0;26,18,12,Times,64,10,0,0,0;12,18,11,Times,66,10,0,0,0;4,18,12,Times,32,10,0,0,0; :[font = input; preserveAspect] f[1][n1_,n2_]:=(r1-a1 n2) n1 f[2][n1_,n2_]:=(-r2+a2 n1) n2 :[font = input; preserveAspect] Jac=Table[D[f[i][n[1],n[2]],n[j]],{i,2},{j,2}] /. {n[1]->r2/a2,n[2]->r1/a1} :[font = input; preserveAspect; endGroup; endGroup] Eigenvalues[Jac] :[font = subsection; inactive; Cclosed; pageBreak; preserveAspect; startGroup] J Linear difference equations and local stability analysis in discrete time :[font = subsubsection; inactive; Cclosed; noPageBreak; preserveAspect; startGroup] Linear difference equations :[font = text; inactive; preserveAspect] 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: ;[s] 42:0,0;51,1;55,0;56,1;168,0;171,1;172,0;210,1;212,0;246,1;248,0;296,1;297,0;300,1;301,0;333,1;334,0;338,1;343,0;346,1;358,0;361,1;367,0;370,1;374,3;375,0;376,1;377,0;380,1;421,0;458,1;460,0;471,1;489,3;490,1;492,0;496,1;600,0;603,1;604,0;615,2;624,0;645,-1; 4:20,13,9,Times,0,12,0,0,0;19,14,10,Times,2,12,0,0,0;1,13,9,Times,1,12,0,0,0;2,18,11,Times,34,10,0,0,0; :[font = input; preserveAspect] f[x_]:=a x NestList[f,x[0],10] :[font = special1; inactive; preserveAspect] Solve the linear difference equation with a = .5, 1.5, -.5, Ð1.5 with n(0) = 100 in each case for 10 generations. ;[s] 7:0,0;47,1;65,0;72,1;75,0;76,1;81,0;115,-1; 2:4,14,10,Times,2,12,0,0,0;3,13,9,Times,0,12,0,0,0; :[font = text; inactive; preserveAspect] 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 ;[s] 111:0,0;34,1;36,0;67,1;69,0;90,1;91,5;92,1;94,0;96,1;101,5;103,2;104,1;105,5;106,1;113,5;116,1;117,5;118,1;129,5;130,6;131,1;133,6;134,1;149,5;150,1;152,0;154,1;159,5;161,2;162,1;163,5;164,1;171,5;174,1;175,5;176,1;187,5;188,6;189,1;191,6;192,1;246,0;249,1;285,6;286,1;288,0;290,1;295,6;296,5;297,2;298,1;299,5;300,1;307,6;308,5;309,7;310,1;311,5;312,1;323,6;325,1;327,6;328,1;332,0;387,3;388,1;390,0;392,1;396,3;397,4;399,1;502,0;505,1;506,0;512,3;513,1;517,0;532,1;533,5;534,1;540,5;541,1;552,6;553,1;556,0;562,3;564,0;594,1;595,6;597,0;670,3;671,0;674,1;675,0;751,4;752,1;758,3;759,8;760,4;762,0;765,1;867,0;913,1;914,6;915,1;918,3;919,6;920,8;921,4;923,0;1043,-1; 9:20,13,9,Times,0,12,0,0,0;43,14,10,Times,2,12,0,0,0;3,22,14,Times,66,12,0,0,0;7,13,9,Times,1,12,0,0,0;4,14,10,Times,3,12,0,0,0;18,18,12,Times,64,10,0,0,0;13,18,11,Times,66,10,0,0,0;1,21,13,Times,64,12,0,0,0;2,18,11,Times,34,10,0,0,0; :[font = text; inactive; preserveAspect] x(t) = S ci lit ei [7] ;[s] 16:0,0;10,1;13,2;14,1;22,3;23,6;24,1;25,7;26,3;28,7;29,8;30,4;31,5;32,7;34,0;131,1;133,-1; 9:2,13,9,Times,0,12,0,0,0;4,14,10,Times,2,12,0,0,0;1,14,10,Times,3,12,0,0,0;2,14,10,Symbol,0,12,0,0,0;1,21,13,Times,65,12,0,0,0;1,13,9,Times,1,12,0,0,0;1,14,10,Symbol,2,12,0,0,0;3,18,11,Times,66,10,0,0,0;1,18,11,Times,34,10,0,0,0; :[font = text; inactive; preserveAspect] (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: ;[s] 87:0,0;17,2;20,7;21,8;22,3;23,5;24,7;25,5;27,0;30,3;31,0;34,2;37,7;38,4;39,0;40,5;41,7;42,0;45,5;46,7;47,9;48,0;50,5;52,0;58,2;61,7;62,8;63,6;64,5;65,7;66,0;70,5;71,7;72,9;73,0;74,5;76,0;82,1;83,2;85,1;86,7;87,2;89,7;90,8;91,6;92,5;93,7;94,3;95,0;100,1;105,0;172,1;173,10;174,0;175,1;176,0;179,1;183,10;184,1;189,0;192,1;193,10;194,1;209,10;210,0;211,1;212,0;215,1;218,0;219,1;221,10;222,1;227,0;229,1;231,10;232,1;318,0;321,1;322,0;346,1;347,10;348,0;351,1;355,10;356,0;363,1;365,0;493,-1; 11:24,13,9,Times,0,12,0,0,0;19,14,10,Times,2,12,0,0,0;5,14,10,Symbol,0,12,0,0,0;3,21,13,Times,64,12,0,0,0;1,22,14,Times,34,12,0,0,0;9,13,9,Times,1,12,0,0,0;2,21,13,Times,65,12,0,0,0;11,18,11,Times,66,10,0,0,0;3,18,11,Times,34,10,0,0,0;2,18,12,Times,32,10,0,0,0;8,18,12,Times,64,10,0,0,0; :[font = input; preserveAspect] x[t_]:=x[t]=A.x[t-1] x[0]={1,1}; A={{1,2},{3,2}}; :[font = input; preserveAspect] Table[x[t],{t,0,10}] //MatrixForm :[font = text; inactive; preserveAspect] They can also be solved by using Eqn [5]: :[font = input; preserveAspect] Table[MatrixPower[A,t].x[0],{t,0,10}] //MatrixForm :[font = text; inactive; preserveAspect] Finally, they can be solved using Eqn [7]: :[font = input; preserveAspect] {vals,vecs}=Eigensystem[A] :[font = input; preserveAspect] {lvals,lvecs}=Eigensystem[Transpose[A]] :[font = input; preserveAspect] Do[lvecs[[i]]=lvecs[[i]]/(vecs[[i]].lvecs[[i]]), {i,2}] :[font = input; preserveAspect] c=lvecs.x[0] :[font = input; preserveAspect] lvecs :[font = input; preserveAspect] xform[t_]:=N[Sum[c[[i]] vals[[i]]^t vecs[[i]],{i,2}]] :[font = input; preserveAspect] MatrixForm[xform[t]] :[font = special1; inactive; preserveAspect] Verify that the formula given by xform generates the same solutions as found above :[font = text; inactive; preserveAspect] 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. ;[s] 46:0,0;152,4;153,1;154,0;195,4;196,0;379,2;381,0;425,2;426,5;427,3;428,0;560,2;561,0;690,2;692,0;982,2;983,5;984,0;985,2;986,0;991,2;993,5;994,2;999,0;1000,2;1001,5;1002,2;1017,5;1018,0;1019,2;1020,0;1026,2;1027,0;1028,2;1029,5;1030,2;1038,5;1039,2;1125,0;1149,2;1150,5;1151,0;1156,2;1158,5;1159,0;1169,-1; 6:16,13,9,Times,0,12,0,0,0;1,21,13,Times,32,12,0,0,0;17,14,10,Times,2,12,0,0,0;1,22,14,Times,66,12,0,0,0;2,18,11,Times,34,10,0,0,0;9,18,12,Times,64,10,0,0,0; :[font = input; preserveAspect] A={{1,2},{-2,1}}; x[0]={1,1}; :[font = input; preserveAspect] {vals,vecs}=Eigensystem[A] :[font = input; preserveAspect] {lvals,lvecs}=Eigensystem[Transpose[A]] :[font = input; preserveAspect] Do[lvecs[[i]]=lvecs[[i]]/(vecs[[i]].lvecs[[i]]), {i,2}] :[font = input; preserveAspect] c=lvecs.x[0] :[font = text; inactive; preserveAspect] 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) ;[s] 21:0,0;86,2;88,0;93,3;94,0;115,1;116,4;117,0;119,2;121,4;122,5;123,0;127,2;128,3;130,2;131,0;134,2;135,0;140,3;142,2;143,0;145,-1; 6:8,13,9,Times,0,12,0,0,0;1,14,10,Symbol,0,12,0,0,0;6,14,10,Times,2,12,0,0,0;3,14,10,Symbol,2,12,0,0,0;2,18,11,Times,34,10,0,0,0;1,18,12,Times,32,10,0,0,0; :[font = input; preserveAspect] r=Abs[vals] //N :[font = input; preserveAspect] theta=Arg[vals] //N :[font = input; preserveAspect] func[i_,t_]:=r[[i]]^t (Cos[theta[[i]] t] + I Sin[theta[[i]] t]) :[font = input; preserveAspect] x[t_]:=Sum[c[[i]] func[i,t] vecs[[i]],{i,2}]; :[font = input; preserveAspect] x[t]=Chop[Simplify[x[t]]] :[font = special1; inactive; preserveAspect] 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. ;[s] 3:0,0;72,1;75,0;234,-1; 2:2,14,10,Times,2,12,0,0,0;1,13,9,Times,0,12,0,0,0; :[font = text; inactive; preserveAspect; endGroup] 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. ;[s] 21:0,0;165,1;167,0;177,1;178,0;281,2;295,0;348,2;353,0;364,3;366,0;482,3;483,0;492,3;493,0;545,3;547,0;761,3;762,0;771,3;772,0;779,-1; 4:11,13,9,Times,0,12,0,0,0;2,13,9,Times,1,12,0,0,0;2,14,10,Times,2,12,0,0,0;6,14,10,Symbol,0,12,0,0,0; :[font = subsubsection; inactive; pageBreak; dontNoPageBreak; preserveAspect; startGroup] Local stability analysis of nonlinear difference equations :[font = text; inactive; preserveAspect; endGroup; endGroup] 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). ;[s] 205:0,0;55,1;56,0;57,1;58,0;63,1;74,0;89,1;90,0;91,1;92,0;95,1;102,0;112,1;120,0;135,1;136,7;137,1;138,0;168,1;176,7;177,1;186,0;202,1;203,0;204,1;205,0;210,1;224,7;226,1;239,7;240,1;245,10;246,1;249,7;250,1;260,7;261,1;265,10;266,1;269,7;270,1;277,0;325,1;326,0;327,1;328,0;334,1;336,10;337,1;340,7;341,1;349,0;392,1;399,0;400,1;401,10;402,1;405,7;406,1;407,0;408,7;409,6;410,1;411,0;414,1;415,0;528,1;529,10;530,1;533,7;534,1;536,0;606,1;608,0;632,1;634,0;653,1;654,8;655,0;656,1;657,0;662,1;664,8;665,1;667,8;668,1;674,8;675,1;686,9;687,1;703,8;704,0;705,1;706,0;711,1;713,8;714,1;716,8;717,1;723,8;724,1;735,9;736,1;776,9;777,0;778,1;779,0;784,1;786,9;787,1;789,8;790,1;796,8;797,1;808,9;809,1;814,0;851,1;852,9;853,1;854,3;855,1;859,4;860,1;864,0;875,1;880,0;915,1;916,9;917,0;918,1;919,0;924,1;926,9;927,1;928,4;929,1;937,9;938,1;939,4;940,7;941,1;944,4;945,1;954,4;955,7;956,1;960,5;961,9;962,2;963,1;964,5;965,1;966,9;967,1;968,4;969,7;970,1;972,5;973,1;974,9;975,1;978,9;979,1;984,7;985,1;988,5;989,9;990,2;991,1;992,5;993,1;994,9;995,1;996,4;997,7;998,1;1000,5;1001,1;1002,9;1003,1;1006,9;1007,1;1008,0;1051,4;1052,0;1053,1;1054,0;1059,1;1060,3;1061,4;1063,1;1067,0;1073,3;1074,4;1075,0;1198,3;1200,0;1572,5;1574,0;1669,5;1670,0;1719,5;1722,0;1820,-1; 11:43,13,9,Times,0,12,0,0,0;88,14,10,Times,2,12,0,0,0;2,22,14,Times,66,12,0,0,0;4,13,9,Times,1,12,0,0,0;10,14,10,Times,3,12,0,0,0;9,14,10,Symbol,0,12,0,0,0;1,22,14,Times,34,12,0,0,0;16,18,11,Times,34,10,0,0,0;10,18,12,Times,64,10,0,0,0;17,18,11,Times,66,10,0,0,0;5,7,5,Times,2,6,0,0,0; :[font = subsection; inactive; Cclosed; pageBreak; preserveAspect; startGroup] K The Poisson distribution :[font = text; inactive; preserveAspect] 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. ;[s] 20:0,0;193,4;194,0;195,4;196,0;201,1;202,2;203,3;204,5;205,0;206,4;207,0;221,4;223,0;238,4;239,0;240,4;241,0;345,3;347,0;531,-1; 6:9,13,9,Times,0,12,0,0,0;1,21,13,Times,32,12,0,0,0;1,22,14,Symbol,32,12,0,0,0;2,14,10,Symbol,0,12,0,0,0;6,14,10,Times,2,12,0,0,0;1,22,14,Times,34,12,0,0,0; :[font = special1; inactive; preserveAspect] 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. ;[s] 5:0,0;97,1;98,0;118,1;119,0;297,-1; 2:3,14,10,Times,2,12,0,0,0;2,14,10,Symbol,2,12,0,0,0; :[font = text; inactive; preserveAspect] 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: ;[s] 139:0,0;122,2;123,0;140,2;141,0;144,2;145,0;160,2;161,0;242,2;243,0;591,2;592,0;593,2;594,0;599,2;600,0;602,2;603,0;605,2;606,0;609,2;610,0;614,2;615,5;616,0;620,2;622,0;623,5;624,6;625,5;626,0;637,2;638,0;656,2;657,1;658,0;704,2;705,0;720,2;721,0;724,2;725,0;735,2;756,0;760,2;761,5;762,0;767,2;768,0;769,5;770,6;771,5;772,1;773,0;786,2;788,0;790,2;791,0;793,2;795,0;797,2;798,0;931,2;932,0;976,2;978,0;1005,3;1007,0;1009,2;1011,0;1086,2;1087,0;1139,2;1142,0;1182,2;1183,0;1201,2;1202,0;1216,2;1217,0;1232,3;1234,0;1259,3;1261,0;1274,2;1275,0;1307,2;1308,0;1310,2;1311,0;1313,2;1315,0;1316,2;1318,0;1325,2;1326,5;1327,0;1328,2;1329,0;1351,2;1352,5;1353,0;1356,3;1357,5;1358,0;1359,2;1360,5;1361,0;1387,2;1388,0;1389,5;1390,6;1391,5;1392,0;1400,3;1401,0;1402,2;1403,0;1404,5;1405,6;1406,5;1407,0;1412,1;1413,4;1414,0;1441,2;1442,0;1443,2;1444,0;1452,1;1453,4;1454,3;1455,5;1456,0;1457,2;1458,0;1501,3;1502,0;2380,-1; 7:60,13,9,Times,0,12,0,0,0;4,21,13,Times,32,12,0,0,0;47,14,10,Times,2,12,0,0,0;7,14,10,Symbol,0,12,0,0,0;2,22,14,Symbol,32,12,0,0,0;15,18,11,Times,34,10,0,0,0;4,18,12,Times,32,10,0,0,0; :[font = input; preserveAspect] data={144,91,32,11,2}; mean = Sum[(i-1)*data[[i]],{i,5}]/280 //N :[font = input; preserveAspect] expected=Table[280 Exp[-mean] mean^x/x!,{x,0,10}] :[font = text; inactive; preserveAspect] 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: ;[s] 118:0,0;194,1;196,0;214,2;215,0;217,1;218,0;366,1;375,0;384,1;385,0;386,1;387,0;388,1;390,0;425,1;426,0;456,1;457,0;503,1;504,0;531,1;532,0;536,1;537,0;556,1;557,0;579,1;581,0;598,1;599,0;604,1;605,0;609,1;610,0;616,1;617,0;643,1;645,0;678,1;680,0;683,1;685,0;689,1;690,0;758,1;759,0;760,1;761,0;763,1;765,0;768,1;769,0;773,1;774,0;775,1;776,0;777,1;779,0;786,2;787,0;789,1;790,0;794,1;795,0;796,1;797,0;802,1;804,0;806,2;807,0;809,1;810,0;829,1;830,0;831,1;832,0;833,1;835,0;838,1;839,0;842,2;843,0;845,1;846,0;847,1;848,0;854,1;855,0;859,1;860,0;861,1;862,0;864,1;865,0;887,1;889,0;904,1;905,0;909,1;910,0;913,1;914,0;918,2;919,3;920,1;921,0;925,1;926,0;963,1;964,0;965,1;966,0;1045,1;1047,0;1183,2;1185,1;1186,0;1188,-1; 4:58,13,9,Times,0,12,0,0,0;53,14,10,Times,2,12,0,0,0;6,14,10,Symbol,0,12,0,0,0;1,14,10,Symbol,2,12,0,0,0; :[font = input; preserveAspect] p[n_,t_]:=Exp[-mu*t] (l t)^n/n! :[font = input; preserveAspect] p[n,0] :[font = input; preserveAspect] Limit[p[0,t],t->0] :[font = input; preserveAspect] D[p[0,t],t] + mu*p[0,t] :[font = input; preserveAspect] Simplify[D[p[n,t],t] - mu*(p[n-1,t]-p[n,t])] :[font = input; preserveAspect] %/.n/n!->1/(n-1)! :[font = text; inactive; preserveAspect] 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: ;[s] 138:0,0;140,4;141,0;142,4;143,0;184,4;185,0;227,4;228,0;231,4;232,0;236,4;237,0;238,4;239,0;265,4;266,0;267,4;268,0;274,4;275,5;276,0;280,1;281,6;282,4;284,5;285,4;286,0;287,7;288,0;329,4;331,0;332,4;333,0;338,1;339,6;340,0;342,8;343,9;344,0;345,4;346,1;347,0;348,8;349,0;350,10;351,0;356,8;357,9;358,0;359,8;360,3;361,1;363,3;364,0;365,8;366,5;367,8;370,9;371,1;381,0;401,1;402,0;426,4;427,0;441,4;442,0;542,4;543,11;544,4;545,0;546,4;547,0;548,2;552,6;553,2;555,0;558,1;559,6;560,0;561,4;563,5;564,8;566,4;567,0;568,4;569,0;571,2;574,6;575,2;577,0;580,1;581,6;582,0;583,4;585,0;586,4;587,0;593,4;594,0;629,4;630,0;649,4;652,11;653,4;655,0;656,4;657,0;658,2;661,6;663,2;666,0;670,4;671,0;672,4;673,0;684,4;685,8;686,0;692,4;693,0;761,4;762,8;763,0;768,8;769,0;770,4;771,0;777,4;778,0;779,4;780,0;791,4;792,0;797,8;798,0;799,4;800,0;851,-1; 12:53,13,9,Times,0,12,0,0,0;8,14,10,Symbol,0,12,0,0,0;6,21,13,Times,64,12,0,0,0;2,22,14,Symbol,32,12,0,0,0;40,14,10,Times,2,12,0,0,0;4,18,11,Times,34,10,0,0,0;7,22,14,Times,66,12,0,0,0;1,12,8,Times,2,10,0,0,0;11,18,12,Times,32,10,0,0,0;3,18,12,Symbol,32,10,0,0,0;1,12,9,Times,0,10,0,0,0;2,7,5,Times,2,6,0,0,0; :[font = input; preserveAspect] G=Exp[(s-1)*mu]; :[font = input; preserveAspect] mean=D[G,s] /.s->1 :[font = input; preserveAspect] variance=D[G,{s,2}] + mean - mean^2 /.s->1 :[font = text; inactive; preserveAspect] 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: ;[s] 13:0,0;82,1;84,0;100,2;101,0;140,2;141,0;158,2;159,0;161,2;162,0;163,2;164,0;235,-1; 3:7,13,9,Times,0,12,0,0,0;1,14,10,Symbol,0,12,0,0,0;5,14,10,Times,2,12,0,0,0; :[font = input; preserveAspect] G[s_]:=P[0]+s*P[1]+s^2*P[2]+s^3*P[3]+s^4*P[4]+s^5*P[5] :[font = input; preserveAspect] Table[D[G[s],{s,r}]/r!,{r,0,5}] /.s->0 :[font = text; inactive; preserveAspect] 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: ;[s] 195:0,0;342,1;343,0;348,1;349,0;394,1;395,2;396,0;397,1;398,0;404,1;405,2;406,0;407,1;408,0;433,1;434,0;437,1;438,0;454,1;455,2;458,0;459,1;460,0;466,1;467,3;470,0;476,1;477,3;478,1;479,3;480,0;505,1;506,2;507,0;508,1;509,0;511,1;512,2;513,0;514,1;515,0;520,1;521,0;526,1;528,0;577,4;578,5;579,0;585,4;586,5;588,0;637,6;642,7;643,8;645,0;646,6;651,7;652,8;654,0;657,6;663,7;664,8;665,6;666,5;667,7;668,8;669,6;671,0;729,4;730,10;731,0;733,4;734,10;735,0;878,1;879,0;939,1;940,2;941,0;951,1;952,2;953,0;954,1;955,0;962,1;963,0;966,1;967,5;968,0;971,1;972,5;973,0;980,1;981,2;982,0;992,1;993,2;994,0;995,1;996,0;1001,1;1002,2;1003,0;1004,1;1005,0;1007,3;1008,0;1045,1;1046,2;1047,0;1095,1;1096,0;1127,1;1128,2;1129,0;1130,1;1131,0;1141,1;1142,2;1143,0;1144,1;1145,0;1159,1;1160,0;1163,1;1164,5;1165,0;1168,1;1169,5;1170,0;1177,1;1178,2;1179,8;1180,0;1182,8;1193,1;1194,2;1195,0;1196,1;1197,0;1201,1;1202,2;1203,0;1209,1;1210,2;1211,0;1215,1;1216,2;1217,0;1218,1;1219,0;1223,1;1224,2;1225,0;1230,1;1231,2;1232,0;1233,1;1234,0;1236,6;1238,0;1240,1;1241,2;1242,0;1247,1;1248,2;1249,0;1250,1;1251,0;1253,6;1254,9;1255,0;1283,1;1284,2;1285,0;1286,1;1287,2;1288,0;1289,1;1290,0;1372,4;1373,0;1545,1;1546,0;1547,1;1548,0;1549,1;1550,0;1560,1;1561,0;1562,1;1563,0;1719,4;1720,0;1837,-1; 11:73,13,9,Times,0,12,0,0,0;61,14,10,Times,2,12,0,0,0;24,18,11,Times,66,10,0,0,0;4,18,11,Times,34,10,0,0,0;6,14,10,Symbol,0,12,0,0,0;7,18,12,Times,64,10,0,0,0;7,18,12,Times,32,10,0,0,0;4,18,12,Symbol,32,10,0,0,0;6,12,9,Times,0,10,0,0,0;1,21,13,Times,32,12,0,0,0;2,21,13,Times,64,12,0,0,0; :[font = input; preserveAspect] G[s_]:=Exp[(s-1)] :[font = input; preserveAspect] Table[Exp[-1]/r!,{r,0,10}] //N :[font = input; preserveAspect] Table[D[G[G[s]],{s,r}]/r!,{r,0,10}] /.s->0 //N :[font = text; inactive; preserveAspect] 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. ;[s] 71:0,0;55,2;56,0;57,2;58,0;163,2;164,3;165,0;166,2;167,0;171,2;172,0;173,2;174,0;175,2;176,0;180,2;181,0;182,2;183,0;206,2;210,0;304,2;305,3;306,0;307,2;308,0;312,2;313,0;314,2;316,3;317,4;319,0;320,2;321,0;342,2;343,0;344,2;345,0;400,2;402,0;470,2;471,0;472,2;473,0;477,2;478,3;479,0;525,2;528,0;532,2;533,0;534,2;535,0;536,2;537,0;548,2;549,0;550,2;551,0;626,2;627,0;628,2;629,0;634,5;635,6;636,5;637,6;642,1;675,7;676,1;685,-1; 8:30,13,9,Times,0,12,0,0,0;2,21,13,Times,32,12,0,0,0;29,14,10,Times,2,12,0,0,0;4,18,11,Times,66,10,0,0,0;1,18,12,Times,64,10,0,0,0;2,18,11,Times,34,10,0,0,0;2,18,12,Times,32,10,0,0,0;1,22,14,Times,34,12,0,0,0; :[font = input; preserveAspect] q[0]=0.0; q[t_]:=q[t]=Exp[q[t-1]-1] :[font = input; preserveAspect; endGroup; endGroup] Table[{t,q[t]},{t,10,100,10}] //TableForm ^*)