3 is required.) In some cases a sequence of period-doublings may occur and this may lead to chaotic behaviour (see section A.2.3) . Another way in which a periodic orbit may exchange stabil ity is through bifurcation into a torus (Hopf bifurcation of the corresponding Poincare map). This occurs when a complex pair of Floquet multipliers moves into or out of the unit circle. Aga in m > 3 is required. Details of this can be found in [111, 124]. From a practical point of view it is probably more informative to return to generating numerical solutions of a model over t ime when these complicated phenomena are encountered in a bifurcation analysis so that the behaviour of the system near these points can be seen explicitly. I w i l l not present the mathematical details of these bifurcations here. Chaos A brief description of chaos was given in section A .2 .3 . Further mathematical details can be found in [11, 16, 124, 128]. Currently there is wide debate as to the practical application of chaos. In ecological systems it is very difficult (probably impossible) to distinguish between stochastic noise and chaotic behaviour (see the preface in [72]). In Stone [114] it is shown that the addition of a single small term to a logistic type model Appendix A. Dynamical systems theory 251 removes the chaotic behaviour. However, in [82, 84] it is demonstrated that chaos is prevalent in many discrete ecological models and that higher order systems display chaotic behaviour more readily than one-dimensional systems. A good overview o f t h e current debate is given in the collection of papers in [72]. A . 3 . 5 M a p s (systems of difference equations) The discussion so far has been restricted to models consisting of systems of ordinary differential equations. A very similar theory can be developed for maps given by (A.2) or (A.3). I w i l l briefly mention some of the differences that occur. More detailed discussions can be found in [124, 128]. Referring to equations (A.2) and (A.3) , an equil ibr ium point occurs when f (x) = x or x i + i = x«. A linear stabil ity analysis can again be done for hyperbolic equi l ibr ium points. For maps a hyperbolic equi l ibr ium point is one for which none of the eigenvalues of the matr ix of part ial derivatives has unit modulus (that is, no eigenvalues have a magnitude of 1). Aga in it is the nonhyperbolic equi l ibr ium points which result in interesting bifurcations. If the linearised matr ix has a single eigenvalue equal to 1 then a l imi t point, trans-cr it ical or pitchfork bifurcation may occur. The bifurcation diagrams are the same as for continuous models. However, it must be remembered that the 'phase portraits' or dia-grams in state space consist of discrete points rather than continuous curves. Examples corresponding to a spiral sink are shown in figure A .25. Consecutive points are labelled in the discrete case. Appendix A. Dynamical systems theory 252 Figure A . 2 5 : State space diagrams of a spiral sink for (a) a continuous model and (b) a discrete model. For maps the special case of a single eigenvalue equal to -1 introduces another type of bifurcation called a period-doubling bifurcation. Figure A.26(a) shows a stable equi l ibr ium point undergoing a, period-doubling bifurcation to become a stable period-2 orbit as fx is increased. This period-2 orbit undergoes a further period-doubling to produce a stable period-4 orbit. Figure A.26(b) gives examples of diagrams in state space corresponding to these situations. Notice that this situation is analogous to the period-doublings of periodic orbits for continuous models. This is not surprising since that theory is based on Poincare maps which are discrete. It is important to note that periodic orbits for discrete systems are different from those for continuous systems. In particular they have integral periods. Consider the second iterate, f2, of the map (A.3): f 2 (x f ) = f ( f ( x 0 ) = f ( x , + 1 ) - x t + 2 . In general, the kth iterate of the map is given by f f c(x t) = xt+Jfe. Suppose there exists a value of x , x , such that f f c(x) = X Appendix A. Dynamical systems theory 253 (a) Xi - « O O 0 c oo • •••••• • • B O O o o o - f i—I—I—I Mi Mi* M2 M2 M3 M (b) X2 1 • 9 , «4 s . A 3« ii) ^2 1. 2* 4». * • • \ lll) 1 z 2 1 2 » 3' 5* 7' 11* 4 6 10 • . • ,14 12 V .9 +-.13 Figure A .26: (a)One-parameter bifurcation d iagram showing per iod-doubl ing bifurcations at // = p\ and p = p*2 for a discrete system. (b)State space diagrams showing the dynamics at (i) \x — \x\ (stable equ i l i b r ium) , (ii) p = \in (stable period-2 orbit) and (i i i) u = fi3 (stable period-4 orb i t ) . but f j ( x ) ^ x f o r j = l , 2 , . . . , k - l . Then x is called an equil ibrium or fixed point of period k. This means that the system has a cycle or periodic orbit whose period is equal to k t ime units. Suppose k = 2 and we have a one-dimensional system. Then / 2(x) = /(/(*)) = x but f(x) + X. Appendix A. Dynamical systems theory 254 Hence, if we start at x then after the first t ime unit we are at f(x) but after the second t ime unit we have returned to x. This is shown in figure A.27. A state space diagram i i i \ i \ — i — -0 1 2 3 4 5 6 7 T i m e Figure A.27: T i m e plot of a period -2 orbi t for a discrete system. for such a situation in two dimensions (that is, for two state variables) is shown in figure A.26(b)( i i ) . In this example both xi and x2 wi l l have t ime plots resembling figure A.27. A n analogue of the Hopf bifurcation also exists for maps. It is sometimes referred to as the Naimark-Sacker bifurcation [124] but is also simply called the Hopf bifurcation for maps. This bifurcation corresponds to a pair of eigenvalues of modulus 1. Instead of a periodic orbit, however, an invariant circle is init iated at this bifurcation point. Whi le geometrically similar to a periodic orbit, the dynamics are different. A state space diagram of an invariant circle is shown in figure A.28 together with a t ime plot. The stabil ity of this circle is intuitively similar to a periodic orbit but the methods of analysis are quite different [124]. There are two possibilities for an invariant circle. Either point 10 coincides with point 1 in figure A.28(a), point 11 coincides with point 2 and so on, or subsequent points are distinct from all the earlier ones but stil l lie on the circle. If the map is allowed to iterate for a long t ime in the latter case, then what is eventually observed in state space wil l Appendix A. Dynamical systems theory 255 (a) X2 (b) X l 1» »9 4t • 6 5« 6' .8 5« «9 2 • • 4 .3 T i m e Figure A .28: (a) State space diagram showing an invariant circle, (b) T i m e plot of the s i tua t ion in (a) in terms of x\. appear to be a continuous circle. In general, the behaviour associated with discrete models is more complicated than that for continuous models because of the bui l t - in t ime delays in the feedback relation-ships [128]. Even one-dimensional maps can exhibit chaotic behaviour. May [82] shows how the behaviour of the discrete analogue of the logistic equation changes from stable equi l ibr ium behaviour, to periodic behaviour, and finally to chaos as the growth rate is increased. The reader is referred to the literature that has been cited for further details on the dynamics of discrete models. A . 3 . 6 Stabi l i t y of bifurcations under perturbat ions The question of robustness or structural stability of a model is an important one. In order to determine how robust a model is, we need to see whether or not perturbing the model alters its qualitative structure. It turns out (see [124]) that l imi t point bifurca-tions, Hopf bifurcations and hysteresis phenomena are stable under small perturbations but transcrit ical and pitchfork bifurcations are not unless constraints or symmetries are preserved by perturbations. That is, small perturbations of the model do not affect whether or not the former three bifurcations occur (they only affect properties such as Appendix A. Dynamical systems theory 256 the parameter value at which the bifurcations occur and the positions of the equi l ibr ium points and periodic orbits) but can affect the occurrence of the latter two. Figure A.29 illustrates the destruction of transcritical and pitchfork bifurcations graphically. It turns out that all bifurcations of one parameter families of equations which have an equi l ibr ium point with a single zero eigenvalue (or single eigenvalue of modulus 1 for maps) can be perturbed to l imi t point bifurcations [49]. Figure A.29: Possible results of per turbing t ranscr i t ical and pitchfork bifurcations (adapted from [111], p.83). However, we cannot make such rapid conclusions when the model contains more than one free parameter. In these circumstances the idea of codimension of a bifurcation becomes important. However, for m > 3 the analysis is very complicated. A n intro-duction to this theory can be found in [124]. The only point I would like to make here is that, because the models in this dissertation have more than three parameters, none of the bifurcations that have been described in this appendix can be dismissed as being unimportant. Appendix A. Dynamical systems theory 257 A.3.7 Mult iple degeneracy A l l the cases discussed so far have assumed simple zero eigenvalues or a single com-plex conjugate pair of eigenvalues with zero real part (or the analogous situations for maps). Higher order singularities, or mult iple degeneracies do, of course, occur but^the behavioural dynamics associated with them can be difficult to interpret and there is st i l l much research being done in this area. Thus, at this stage, it is probably best to either solve the system of equations numerically or to generate phase portraits for parameter combinations in a region surrounding these complex points rather than struggling with the details of the bifurcation structure. For those interested some of these higher order degeneracies are investigated in [124]. A . 4 Conclusion ( The introduction to dynamical systems theory given in this appendix has been mainly intuit ive and not mathematical ly rigorous. Although the concepts are fairly simple, the mathematics involved in studying a particular system of equations can be quite formidable. Appendix B describes how computers can be useful in this regard. Appendix B Numerical details B . l Introduction This appendix describes some of the computer software that is available for analysing systems of equations. Packages such as A U T 0 8 6 [28], A U T 0 9 4 [31], Interactive A U T O [117], X P P A U T [35] and some others that wi l l be mentioned later, enable particular solu-tions to be 'continued' as a parameter is varied in order to produce a bifurcation diagram. In other words, these continuation programs trace out the location of equi l ibr ium points and periodic orbits as a parameter is varied. Bifurcation points that are encountered along the way are also detected and classified. In this way a whole range of different modes of model behaviour can be obtained with much greater ease than if the mathe-matics had to be done by hand. The process can be repeated for different parameters in order to obtain a more comprehensive picture of the qualitative behaviour associated with different regions in the parameter space. Section B.2 introduces some of the techniques used by continuation programs for continuing equi l ibr ium points and determining the stability of the solution branches (see section A.2.4) . Methods for detecting bifurcations along these branches are also described wi th particular reference to A U T O 1 . The section is intended as a brief overview. More H w i l l use the abbrevia t ion A U T O to refer to A U T 0 8 6 , A T J T 0 9 4 and Interactive A U T O as the lat ter is i n essence jus t a graphica l interface for A U T 0 8 6 and A U T 0 9 4 is an updated version of A U T 0 8 6 together w i t h a graphica l interface. Even when using X P P A U T I w i l l refer to A U T O when t a lk ing about bifurcat ion diagrams since X P P A U T generates these diagrams through a graphica l interface w i t h A U T 0 8 6 . In cases where I need to dis t inguish between the packages, I w i l l use their fu l l names. 258 Appendix B. Numerical details 259 detailed discussions as well as a comprehensive list of references can be found in Allgower and Georg [4] and Seydel [111]. Following this theoretical introduction section B.3 describes the capabilities and l i m -itations of the packages that I used and section B.4 describes how to obtain t ime plots, phase portraits and bifurcation diagrams. Final ly , section B.5 gives a few pointers and warnings regarding the use of some of the packages. These are based on my experiences gained through analysing the models in the main body of the thesis. B .2 Theory B.2.1 Continuation methods As has already been mentioned, continuation or path-following methods generate a chain of solutions (equil ibrium points, periodic orbits or bifurcation points) as a parameter is varied. A typical path-following method is the predictor-corrector method. This involves the repetition of two different steps. The predictor step approximates the next point on the curve, often by using the direction of the tangent to the curve (Euler predictor) [4]. A number of iterative steps (called corrector steps) then a im to improve this approximation and bring it back to the actual curve [4]. Typical ly Newton or gradient type methods are used in this step [4]. Some form of parameterisation of the curve is required for these steps. The obvious choice is the control parameter (that is, the parameter being varied) as it has physical significance. However, this leads to difficulties at l imit points (see section A.2.13) [111]. A n alternative is to choose another variable which involves adding another equation to the system. This extended system can then be solved using the predictor-corrector methods mentioned above [111]. Popular choices for this alternative parameter are arclength or a pseudo-arclength parameter proposed by Keller [67]. A U T O uses the latter choice and Appendix B. Numerical details 260 details are given in [30, 67]. The accuracy of the predictor-corrector method depends on the choice of steplength. In general, shorter steplengths lead to greater accuracy (provided that they are not so small that computer round-off errors become large) but they are more costly in terms of t ime. In some cases the objective is jusf to follow the curve as rapidly and safely as possible unt i l a crit ical point (such as a bifurcation point) is reached [4] and then a smaller steplength is needed for greater accuracy. Thus, for an efficient algorithm, the steplength needs to be adaptive and not fixed [4, 111]. Ideally a continuation method should also allow the user to have some control over the choice of steplength. A U T O fulfills both criteria. Stepsizes are changed automatically in the program depending on the speed of convergence of Newton's method (that is, depending on the number of iterations required to fulfi l l the stopping criteria). M a x i m u m and m i n i m u m stepsizes and convergence criteria are given by the user. There are many different ways in which choices of predictor, corrector, parameteri -sation and step control can be combined to produce a continuation method. Because of this no numerical comparison of different path-following methods has so far been done and, hence, no particular method can be recommended exclusively [111]. Simple E u -ler predictors together with Newton-type correctors have been found to be satisfactory in many circumstances [4]. Because of the stabil ity of Newton correctors, more stable higher order predictors based on polynomial interpolation (instead of on the direction of the tangent to the curve as with an Euler predictor) are often advantageous [4]. No continuation method can guarantee that all possible solutions wi l l be found in a given example [111]. Isolated branches (branches which are not attached to other branches v ia bifurcation points) are very difficult to detect. It is suggested in the A U T 0 8 6 manual [30] and by Seydel [111] that t ime integration of the governing system of equations using random in i t ia l data may be worthwhile for generating a starting point on an as yet Appendix B. Numerical details 261 undetected solution branch. D S T O O L [10] can also be useful for locating equi l ibr ium points which lie on these isolated branches. B .2 .2 Detection of bifurcations As a continuation method traces out a path, it needs to be able to detect bifurcation points. Techniques for doing this can be divided into direct and indirect methods [111]. Direct methods involve enlarging the original system of equations by including addit ional equations which characterise the bifurcation point. Indirect methods on the other hand utilise data obtained during a continuation together with a test function. The latter are generally recommended in practical computations (and are used in A U T O ) as direct methods involve solving much larger systems of equations and have higher storage re-quirements. Indirect methods do have more difficulty achieving high accuracy but when discretisation errors are present and when bifurcation points are unstable to perturba-tions, this greater accuracy is not needed [111]. If higher accuracy is required then a direct method can be applied once an indirect method has obtained an approximate result. A n indirect method detects bifurcation points using 'test functions' which are evalu-ated during branch tracing [111]. A bifurcation point is indicated by a zero of the test function, r. (In practice, an algorithm checks for a change of sign of r.) For simple bifurcation points a natural choice for r is the max imum of all real parts a ; of the eigen-values of the Jacobian matr ix A (see section A.3.2) . However, this choice may not be smooth and it does not signal bifurcations in which only unstable branches coalesce [111]. This problem can be overcome by choosing r = with |afc| = m i n . . . , |am|}. This choice also detects Hopf bifurcations. In general, the accuracy of T depends on the accuracy with which A is evaluated. There are many possibilities for the choice of test function. However, not much Appendix B. Numerical details 262 attention has been paid to which drawbacks are significant or which test function is best for which purpose [111]. Some authors have studied classes of test functions for various types of singular points (see [4] for references). Once a bifurcation or branch point has been detected, a method needs to be found for switching branches. A l l that is needed is a single point on the new branch as then the continuation method can be used to trace out the rest of the branch. Predictors and correctors can again be used for switching branches. Seydel [111] and Allgower and Georg [4] describe a few different approaches and give further references. A U T O uses a method suggested by Keller [67]. Doedel [30] notes that this method performs well in most applications although difficulties can occur if the angle of intersection of the branches is very small . A t a Hopf bifurcation point periodic orbits are introduced. Once a solution point on the branch of periodic solutions has been located, a continuation procedure can be used to trace out the branch. After imposing an integral condition in order to fix the phase of the orbit, the continuation procedure becomes a special case of the path-following techniques that have already been discussed. Details of the method as well as further references can be found in the A U T 0 8 6 manual [30]. B.2.3 Stability Stationary branches The stabil ity of stationary branches (that is, of continuations of equi l ibr ium points) is de-termined from the real parts of the eigenvalues of the Jacobian matr ix A. The calculation of eigenvalues (for example, v ia Q R factorisation or L U decomposition methods [4, 111]) can be very t ime-consuming especially when accurate approximations are required, such as near bifurcation points, and when the system of equations is large. However, the Appendix B. Numerical details 263 fact that the eigenvalues play an important role in the detection of bifurcations, as well as in determining stability, provides justification for using more accurate techniques. A U T O uses the I M S L subroutine E I G R F for computing the eigenvalues of a general real matr ix [30]. Interactive A U T O [117] assigns different colours to continuation branches depending on the number of unstable eigenvalues while X P P A U T [35] uses thick lines to indicate stable branches and thin lines to indicate branches having one or more unstable eigenvalues. Periodic branches It is the values of the Floquet multipliers (see page 248) that determine the stabil ity of periodic orbits. These multipliers are eigenvalues of a particular matr ix and, thus, an eigenvalue solver is again required. Since Floquet multipliers can be very large or very small in value, there may be a loss of accuracy when evaluating them numerically. This is especially true for unstable orbits [111] and near orbits of infinite period [30]. In the A U T 0 8 6 manual [30] it is noted that orbits normally retain their accuracy even when the computation of Floquet multipliers (and hence stabil ity determination) breaks down. In X P P A U T and A U T 0 9 4 the routine for calculating Floquet multipl iers has been improved and is more accurate than that in A U T 0 8 6 , and hence in Interactive A U T O . Floquet multipl iers are also used for detecting higher order periodic bifurcations such as period-doubling bifurcations or bifurcations to tori . B.3 Available computer packages In this section I list the capabilities of a number of computer packages and also list some advantages and disadvantages related to choosing one package over another. Appendix B. Numerical details 264 B.3.1 A U T O 8 6 Since Interactive A U T O , X P P A U T and A U T 0 9 4 are all based on A U T 0 8 6 [28], their basic capabilities are very similar. Although I have not used A U T 0 8 6 directly (it runs as a batch process), I wi l l describe its capabilities and l imitations since most of them apply to the abovementioned packages. I highlight the differences between the packages in subsequent sections. Capabilities (These are taken directly from the A U T 0 8 6 manual [30].) A U T O can do a l imited bifurcation analysis of algebraic systems of the form f ( x , A ) = 0 , x , f e 7 e m (B.6) where A denotes one or more free parameters, and of systems of ordinary differential equations of the form x = f ( x ,A ) , x , f e 7 t m (B.7) It can also do certain continuation and evolution computations for the diffusive system x t = D x u u + f ( x , A ) , x , f e 7 l m (B.8) where x = x ( u , £ ) and D denotes a diagonal matr ix of diffusion constants. For the algebraic system (B.6) A U T O can: • trace out branches of solutions. • locate bifurcation points and compute bifurcation branches. • locate l imit points (saddle-node bifurcations) and continue these in two parameters. Appendix B. Numerical details 265 • do al l of the above for fixed points of the discrete dynamical system x < + 1 = f (x , ,A) . (B.9) • optimisation: find extrema of an objective function along the solution branches and successively continue such extrema in more parameters. For the ordinary differential equations (B.7) A U T O can: • compute branches of stable or unstable periodic solutions and Floquet mult ipl iers. Starting data for the computation of periodic orbits are generated automatical ly at Hopf bifurcation points. • locate l imi t points, transcritical and pitchfork bifurcations, period-doubling bifurca-tions and bifurcations to tori along branches of periodic solutions. Branch switching is possible at transcrit ical, pitchfork and period-doubling bifurcations. • continue Hopf bifurcation points, l imit points and orbits of fixed period in two parameters. • compute curves of solutions to (B.7) on the interval [0,1] subject to general nonlinear boundary or integral conditions. • locate l imi t points and bifurcation points for such boundary value problems. Branch switching is possible at bifurcation points. Curves of l imit points can be computed in two parameters. For the parabolic system (B.8) A U T O can: • detect bifurcations to wave train solutions of given wave speed from spatially ho-mogeneous solutions. These are detected as Hopf bifurcations along fixed point branches of a related system of ordinary differential equations. Appendix B. Numerical details 266 • trace out the branches of wave solutions to (B.8) and detect bifurcations. The wave speed c is fixed but the wave length L wi l l normally vary. • trace out branches of waves of fixed wave length L in two parameters. If L is large, then one gets a branch of approximate solitary wave solutions. • do t ime evolution calculations for (B.8) with periodic boundary conditions on [0,L]. In this thesis systems of the form (B.7) and (B.9) are analysed. The discretisation used in A U T O to approximate ordinary differential equations and for calculating periodic solutions is orthogonal collocation with 2 , . . . , 7 Gauss collocation points per mesh interval. The mesh automatically adapts to the solution so that a measure of the local discretisation error is equi-distributed. Also, the adaptive mesh guards to some extent against computing spurious solutions. When spurious solutions do occur they are often easy to recognise by the jagged appearance of the solution branch. Some general l imitations of A U T 0 8 6 are listed below. Limitations The following difficulties are noted in the manual [30]: • degenerate (multiple) bifurcations cannot be detected in general. A lso , bifurcations that are close together may not be noticed when the pseudo-arclength step size is not sufficiently small . • Hopf bifurcation points may go unnoticed if no clear crossing of the imaginary axis takes place. This may happen when other real or complex eigenvalues are near the imaginary axis and when the pseudo-arclength step is large compared to the rate of change of the crit ical eigenvalue pair. A n often occurring case is a Hopf bifurcation close to a l imit point (saddle-node bifurcation). Appendix B. Numerical details 267 • similarly, Hopf bifurcations may go undetected if switching from real to complex conjugate, followed by crossing of the imaginary axis, occurs rapidly with respect to the pseudo-arclength step size. • secondary periodic bifurcations may not be. detected for very similar reasons. • for periodic orbits the numerical output should be checked to make sure that points labelled as bifurcation points, l imit points, period-doubling bifurcations or bifurca-tions to tori have been classified correctly. Some of the above problems may be solved by decreasing the m i n i m u m step size, dsmin, to allow A U T O to take smaller steps. This is particularly helpful when two continuation branches or two bifurcation points are very close together. In the former situation A U T O needs to take small steps to ensure that it does not switch branches during the calculation. Decreasing dsmax may also help as this prevents A U T O from taking large steps and missing important bifurcations. As noted above A U T O may have some difficulty with identifying bifurcations on periodic orbits. Problems arise when Floquet multipliers are close to the unit circle. Also, unstable orbits can be difficult to locate, especially when they are close to a stable orbit or equi l ibr ium. Accuracy may be increased by decreasing dsmin or increasing ntst, the number of mesh points used in the discretisation of the periodic orbits. Some other l imitations of A U T O which have resulted from doing bifurcation analyses on a few models are as follows: • transcrit ical and pitchfork bifurcations cannot be continued in two parameters be-cause the determinant of the Jacobian matr ix A is zero at these bifurcation points. Other continuation and bifurcation packages wi l l have the same l imitat ion. A possi-ble solution is to do a number of one-parameter continuations with different (fixed) Appendix B. Numerical details 268 values of the second parameter and then to plot the bifurcation points obtained from each continuation in two-parameter space. A n approximate curve can be drawn through these points. This is done in chapter 2. • error messages from A U T O can be misleading (as with most computer packages!). In many cases the problem is solved by checking the driving program for typing errors and variable or parameter names which begin with letters between h and o as such quantities are assumed to be integers by default. B.3.2 Interactive A U T O Capabilities In addition to the capabilities listed for A U T 0 8 6 , Interactive A U T O [117] allows the user to change the program constants interactively and to observe the development of a bifurcation diagram while a calculation is in progress. The corresponding eigenvalues are shown simultaneously in a separate graphics window. The advantages and disadvantages listed below relate mainly to the suitabil ity of this package for analysing ecological models and are included for comparison wi th the other packages. Since the packages were not set up specifically for the models that are analysed in this thesis, these comments are not necessarily criticisms. Advantages • This package has good on-screen graphics. Different colours are used to indicate the number of unstable eigenvalues corresponding to a particular branch. The locations of the eigenvalues relative to the real and imaginary axes and the unit circle are also shown. Bifurcation diagrams can be viewed in three dimensions if desired and the mouse can be used to shift diagrams and to zoom in and out. Appendix B. Numerical details 269 Disadvantages • Bifurcat ion diagrams cannot be printed out or saved as postscript files for later printing. • Each t ime the model equations are altered the driving routine needs to be recom-piled. • The package only runs on Personal Iris and Iris 4D workstations (instead of on all systems supporting X-windows as for the other packages). • Interactive control of the graphics output is l imited. Once a one-parameter bifur-cation diagram has been generated it is not possible to view it with a different variable on the y-axis. The scales of the axes are also not visible. B . 3 . 3 X P P A U T This package (as well as the tutorial) can be obtained v ia anonymous ftp from ftp.math.pitt .edu and is in the directory pub /hardware. Capabilities The A U T O interface allows most of the capabilities of A U T 0 8 6 to be enjoyed. How-ever, the a im of making A U T O easier to use has led to some restrictions which wi l l be mentioned below. This package also solves systems of equations numerically (there is a choice of a l -gorithms) and generates time plots and phase portraits (in two and three dimensions). Appendix B. Numerical details 270 Hardcopies of these'diagrams can be obtained. W i t h reference to phase portraits it is pos-sible to obtain nullclines, arrows indicating the direction of flow, as well as the positions and stabilities of equi l ibr ium points (singular points). Other capabilities include curve-fitting of data, spreadsheet type data manipulat ion and generation of histograms. Advantages • T i m e plots, phase portraits and bifurcation diagrams can all be generated using the same computer package. • The A U T O interface is easy to use—the number of constants that need to be altered has been reduced. • Hardcopies of bifurcation diagrams can be obtained through saving them as postscript files. The data can also be saved so that it can be read into other graphics packages. • The package has good on-screen graphics. The development of a bifurcation d i -agram can be seen while a calculation is in progress. The eigenvalues are shown simultaneously relative to the real and imaginary axes and the unit circle. It is possible to zoom in on specific regions of a diagram and the most recent one-and two-parameter bifurcation diagrams viewed are kept in memory. Once a one-parameter bifurcation diagram has been generated the variable on the y-axis can be altered. It is also possible to plot the period or the frequency of a periodic orbit as a function of the bifurcation parameter. M i n i m a as well as m a x i m a of periodic orbits can be plotted simultaneously. • A U T 0 8 6 ' s Floquet multipl ier routine has been improved. • It is possible to start continuations from a numerically calculated periodic orbit. Appendix B. Numerical details 271 • Three-dimensional phase portraits are possible. • The driving program is easy to write for most systems and is compiled automatical ly when the X P P A U T command is given. • W h e n calculating t ime plots and phase portraits, data information for auxil iary variables (variables which are subcomponents or composite functions of the state variables) can also be observed. • The package runs on any X-windows system as well as on L I N U X . • There is a comprehensive tutorial on the World Wide Web to help the user become acquainted with the package. The address is: ftp://mthbard. math, pitt.edu/pub/bardware/xpptut/start. h tml . D i s a d v a n t a g e s • W h e n generating phase portraits only one equi l ibr ium point is located at a t ime and the in i t ia l point often has to be quite close by. • The A U T O interface is not set up for discrete equations and error tolerances for A U T O cannot be altered. Also, the automatic detection of l imi t points cannot be turned off in this version of A U T O . This may cause difficulties when calculat-ing periodic orbits as at a l imit point of a periodic orbit two Floquet mult ipl iers equal 1 and this affects the convergence properties of the continuation algorithm. Decreasing dsmin and dsmax may help. • There is l imited control of the appearance of printouts. The stabil ity nature of a particular branch in a bifurcation diagram is often obscured when consecutive continuation points are very close together. Dashed lines for unstable branches then Appendix B. Numerical details 272 appear as solid lines and, instead of indiv idual dots, periodic branches become thick lines. In some cases this can be overcome by increasing the m a x i m u m step size, dsmax. A lso, the data for a diagram can be saved separately and then read into another graphics package. When two-parameter bifurcation diagrams are printed out, data from the previously generated one-parameter diagram is superimposed. B.3.4 A U T 0 9 4 Capabilities This package is an updated version of and graphical interface for A U T 0 8 6 . The graphical interface allows program constants to be altered interactively and bifurcation diagrams can be viewed. Advantages • There is a help menu which describes the functions of the program constants. • Demonstration examples show how the various capabilities of the package can be used. • Model equations can be changed interactively. • The package runs on any X-windows system. • It is possible to start a continuation from a numerically calculated periodic orbit. • The routine for calculating Floquet multipliers has been improved. Disadvantages • Printouts of bifurcation diagrams cannot be obtained. Appendix B. Numerical details 273 • The on-screen graphics are more cumbersome to use than those of Interactive A U T O and X P P A U T . Changes to a bifurcation diagram can only be viewed at the end of a calculation using a separate command instead of while the program is running. The way in which bifurcation points are labelled also makes the diagrams difficult to read. The disadvantages were the main reason why I did not make use of A U T 0 9 4 . B .3.5 D S T O O L This package can be obtained v ia anonymous ftp from macomb.tn.cornell.edu and is in the pub/dstool directory. Capabilities This package generates t ime plots and two-dimensional phase portraits for both discrete and continuous systems of equations. It calculates equil ibr ium points and their stabilities as well as stable and unstable manifolds for saddle points. There is provision for exten-sions to the package—three-dimensional graphics as well as continuation and bifurcation routines may be incorporated in the near future. Advantages • A mouse can be used for specifying starting points for t ime plots and phase por-traits. This is very convenient and speeds up the generation of these diagrams considerably. • The package is good at locating periodic points for discrete models. Appendix B. Numerical details 274 • Printouts of diagrams can be obtained. • For ordinary differential equation models trajectories can be calculated forwards or backwards. • The package runs on any X-windows system as well as on L I N U X . Disadvantages • Three subroutines need to be modified each t ime a new model is entered. • The package sometimes crashes when calculating equi l ibr ium points for difficult parameter values. B.3.6 Other packages Some other packages are available for analysing dynamical systems. Part of the following list can be found in [4]. 1. A L C O N [27]. This is a continuation method for algebraic equations f ( x , A) = 0. L i m i t points and simple bifurcations can be computed on demand. 2. B I F P A C K [112]. This is an interactive program for continuation of large systems of nonlinear equations. Bifurcation points are also detected. 3. D Y N A M I C S [96]. This package iterates maps, solves differential equations, and plots trajectories. It runs on both I B M P C ' s and U N I X workstations which support X-windows. 4. L O C B I F [68]. This is an interactive program designed for multiparameter bifur-cation analysis of equi l ibr ium points, l imit cycles and fixed points of maps. A t Appendix B. Numerical details 275 present it is set up for I B M P C ' s but a U N I X version is in process and is being incorporated into D S T O O L . 5. P A T H [66]. This software package for dynamical systems can apparently handle much larger systems of ordinary differential equations than A U T O . 6. P H A S E R [54]. This package generates phase portraits for continuous and discrete dynamical systems. It runs on I B M P C ' s . 7. P I T C O N [100, 101]. This is a Fortran subprogram for continuation of equi l ibr ium points and for detecting l imit points. As has already been mentioned, no comprehensive comparison of different techniques has been done because of the enormity of the task. None of the abovementioned packages analyses a model at the touch of a button and complementary analytical techniques, as well as other numerical techniques, are st i l l needed in order to obtain a complete understanding of any model. Parameter studies are as much an art as a science [111]. In the next section I explain how to generate t ime plots, phase portraits and bifurcation diagrams using some of the packages I have mentioned. B.4 U s i n g the packages Although this section gives some guidelines for using the various packages, it is important that anyone wishing to use these packages reads the relevant user manuals to find out the exact commands for performing various tasks. Many of the manuals have introductory examples for the reader to work through and these are invaluable for getting acquainted with the capabilities of the package. There is a comprehensive tutorial for using X P P A U T on the Wor ld Wide Web. The address is: f tp://mthbard. math, pitt.edu/pub/bardware/xpptut/start. h tml . Appendix B. Numerical details 276 Time plots These can be obtained using either D S T O O L or X P P A U T . The user s imply chooses time as the variable to be plotted on the x-axis and one of the state variables for the y-axis. After entering in i t ia l values for the state variables and the t ime period over which to simulate the model , the ' run' or 'go' command can be chosen. The process may be repeated after changing one or more in i t ia l values or parameter values. These steps can al l be done interactively. Both packages have a number of choices for the numerical algorithm that is used to calculate solutions. M i n i m u m and max imum or fixed stepsizes can be chosen as well as error tolerances. In most cases the default choices are quite adequate. If the mouse is used to choose an in i t ia l point, then only the in i t ia l t ime and the value of the state variable on the y-axis wi l l be altered automatically. Values for the other state variables wi l l remain unchanged unless new values are typed in . When analysing a discrete system using D S T O O L , changing the t ime increment to 1 in the 'orbit ' window and using the 'continue' icon allows the user to determine the period of a cycle. The amplitude of the cycle can be deduced by looking at the state variable values in the 'settings' window that correspond to the max imum and m i n i m u m of the cycle. Phase portraits These can again be generated using either D S T O O L or X P P A U T . A l though three-dimensional portraits are possible in X P P A U T , I w i l l only discuss two-dimensional por-traits as I find them much easier to interpret and, hence, find them more useful in most situations.- (Note that the dimension of the system can be greater than two but the results are projected into a two-dimensional plane.) Appendix B. Numerical details 277 In order to generate a phase portrait one of the state variables needs to be chosen for plott ing on the x-axis and another for the y-axis. For systems of dimension greater than two the solution trajectories wi l l be projected into this plane. Init ial values can be typed in manually or the mouse may be used. Only the values of the state variables shown on the axes wi l l be altered when using the mouse. The user can again choose between the various numerical methods for calculating solution trajectories. In addit ion to solution trajectories the positions of equi l ibr ium points can be shown in phase diagrams. When asked to find equil ibr ium points (also called fixed points or singular points) both D S T O O L and X P P A U T automatically calculate the eigenvalues and hence the stabil ity of these points. X P P A U T tends to locate one singular point at a t ime and often the in i t ia l point has to be fairly close to the fixed point if the search is to be successful. D S T O O L uses a Monte Carlo technique to generate a specified number (the default is 10) starting points and then searches for fixed points beginning at these starting values. This method is fairly efficient and choosing the 'f ind' option in the 'fixed point' window a few times generally locates all the relevant points. If desired, the user can choose the in i t ia l point from which to begin a calculation instead of using the Monte Carlo technique. For discrete models D S T O O L also finds periodic points. For continuous systems both packages also calculate one-dimensional stable and un-stable manifolds associated with saddle points in the plane. These help del imit domains of attraction. X P P A U T can calculate nullclines (curves showing where each differential equation is equal to zero) and display these in the phase portrait. Equ i l ib r ium points are located at intersections of nullclines corresponding to different state variables. One-parameter bifurcation diagrams For systems of ordinary differential equations one-parameter bifurcation diagrams can be generated using Interactive A U T O , A U T 0 9 4 and X P P A U T . For discrete systems either Appendix B. Numerical details 278 of the former two packages can be used. The first step is to choose the bifurcation parameter to be plotted on the x-axis. One of the state variables can be chosen for the y-axis. This choice must be specified in the driving program for Interactive A U T O but can be done interactively in the other two packages. In addit ion to choosing the scales for the axes there are a number of other program constants which need to be set. These govern, for example: the length and type of continuation; the output to the screen; the detection of l imit points; error tolerances, steplengths and mesh intervals for the various numerical routines; and a number of other aspects of the computation. The various possibilities are listed in the A U T 0 8 6 manual [30] as well as in the H E L P menu in A U T 0 9 4 . The example or demonstration programs also give an idea of appropriate values. In X P P A U T some of these constants have been preset to simplify the use of A U T O . This is very convenient in most situations but can be restrictive in others. In order to begin generating a bifurcation diagram a fixed point (equil ibrium point) , corresponding to a particular parameter set, is required. Such a point can be determined analytical ly (where possible) or numerically using X P P A U T or D S T O O L . A continuation can then be started from this point. A t transcrit ical, pitchfork and period-doubling bifurcations A U T O automatically calculates the various intersecting branches. A t Hopf bifurcation points the user must specifically initiate the calculation of periodic orbits by choosing the relevant restart label. In Interactive A U T O and A U T 0 9 4 one of the program constants also needs to be altered. In order to extend any continuation branch across a wider range of parameter values, the endpoint of the branch can be chosen as the restart value. For discrete systems period-doubling bifurcations are labelled as Hopf bifurcations by A U T O but the period-2 branches emanating from this bifurcation cannot be calculated directly. The second iterate of the model must be used for this purpose. Since a (period-1) Appendix B. Numerical details 279 equi l ibr ium point of the original model is also an equi l ibr ium point of the second iterate of the model , both the period-1 and the period-2 equi l ibr ium point branches wi l l be traced out using this latter model. However, for higher order period-doublings higher order iterates of the model are required and the process is more tedious. In such cases as well as for complicated discrete models whose second iterate is difficult to calculate, D S T O O L can be used to generate an approximate bifurcation diagram. Equ i l ib r ium points can be calculated at regular intervals across a range of parameter values and their coordinates recorded. These points can then be plotted using some other graphics package, such as G N U P L O T [125], to obtain an approximate bifurcation diagram. This approach is used in chapter 7. Two-parameter bifurcation diagrams Once a l imi t point or Hopf bifurcation point has been detected in a one-parameter con-t inuation, it can be continued in a second parameter. This is done by choosing the appropriate restart value corresponding to the bifurcation point, a second parameter for the y-axis, and either altering the relevant program constant in Interactive A U T O and A U T 0 9 4 or choosing the two-parameter option in X P P A U T . L i m i t cycles of fixed period can also be continued in two parameters. This is useful for approximating homoclinic or heteroclinic bifurcation curves since homoclinic and heteroclinic orbits have infinite period. Designating a U S Z R function to locate an orbit of high period, say period=100 or 1000, when constructing a one-parameter bifurcation diagram allows a two-parameter continuation of this orbit to be done 2 . The resulting two-parameter curve gives the required approximation to the curve of homoclinic or heteroclinic bifurcation points. This technique is used in section 2.4. When continuing 2 A po in t sat isfying the U S Z R function w i l l be given a label wh ich can be used as a restart value for the two-parameter cont inuat ion. Appendix B. Numerical details 280 these curves in two parameters, detection of l imit points should be turned off as it may lead to spurious bifurcation points. This is done automatically in X P P A U T . Any bifurcation points should also be checked against the numerical output to see whether an eigenvalue or Floquet mult ipl ier has, in fact, changed sign. B.5 Pointers and warnings As with al l computer packages, the ones that I have discussed do not work exactly as one might wish in a given situation. Each model is unique and requires a slightly different approach whereas a computer package is built for more general use. In this section I list some problems and suggestions arising from my experiences with the packages. Many of these wi l l only make sense once the examples in the main chapters of the thesis have been studied, but I have included them here so that most of the computer technicalities are confined to one place for ease of reference. • W h e n calculating fixed points using D S T O O L , the package sometimes hangs and wi l l not respond to input. I have not been able to locate the cause of this but quitt ing D S T O O L and restarting the package, although frustrating, does solve the problem. • W h e n starting a fixed point continuation in A U T O , the error tolerances cannot be set lower than the accuracy of the state variables that have been given as the starting point. If the starting points are only known to low accuracy, then a short continuation can be done on low accuracy (high error tolerance). The continuation can then be restarted at higher accuracy using the values calculated by A U T O . (Accuracy is increased by decreasing dsmax, the m a x i m u m step length, decreasing the error tolerances, or increasing ntst, the number of discretisation points.) Appendix B. Numerical details 281 • A U T O wi l l generate bifurcation diagrams faster if the accuracy is lower. In many cases the results are st i l l sufficiently accurate, however some bifurcation points may need to be checked with greater accuracy calculations. It is a good idea to check the eigenvalues and Floquet multipliers in the numerical output files created by A U T O to see whether a bifurcation point has been correctly labelled. A lso , if there is a sudden jump in a continuation, or a curve becomes very jagged, then the continuation should be repeated using greater accuracy. • The choice of step size depends on the extent of the parameter range over which changes in behaviour occur. In most cases ds = 0.02 and dsmax = 0.05 are good choices. However, I used ds = 0.0001 and dsmax between 0.001 and 0.01 for the population genetics models as the bifurcations take place within fairly small parameter ranges. Scaling the equations and parameters (as is done in chapters 3 and 4) circumvents this problem. • Bifurcation diagrams can become very complicated even when studying simple models. However, not al l continuation branches are always of interest. For example, some may correspond to negative or zero values of the state variables. A lso , when studying a particular aspect of a model, some branches may be superfluous. It is important to look at the numerical output which is printed to the screen during a continuation so that one can keep track of which branches are relevant and which are not. • It is a good idea to generate a number of starting points at a variety of parameter values to ensure that a complete bifurcation diagram is obtained and that isolated branches have not been overlooked. D S T O O L is convenient for this purpose. Appendix B. Numerical details 282 • Results can often be checked using more than one package. For example, a b i -furcation diagram generated by A U T O can be checked by choosing a number of parameter combinations corresponding to different qualitative regions and generat-ing phase portraits for these combinations (using X P P A U T or D S T O O L ) to check the dynamics. This is most convenient in X P P A U T where the bifurcation dia-grams and phase portraits can be generated within the same package. T ime plots and phase portraits can be checked by choosing different numerical methods to calculate solutions. Both D S T O O L and X P P A U T offer a variety of methods. B y looking at the eigenvalues corresponding to an equi l ibr ium point one can check that the stabil ity of the point has been correctly labelled. • It is easy to get side-tracked into studying the bifurcation structure of a model and into studying complicated behavioural changes instead of concentrating on phenomena which are of biological interest. In general, sharp boundaries and the exact values at which bifurcations occur are not important as biological parameters are hardly ever known with certainty. The general behaviour and the types of changes that can occur are of greater practical interest. For example, l imi t cycles of small amplitude wi l l be indistinguishable from equi l ibr ium points due to the natural variation of field data. Also, a system wi l l easily be perturbed from a stable node with a small domain of attraction and, hence, such a phenomenon may be of minor interest. It is also important to look at t ime plots corresponding to fixed parameter combinations to see how quickly a system approaches the l imi t ing behaviour indicated in a bifurcation diagram. If a system takes a long t ime to attain its equi l ibr ium configuration then the transient dynamics may be more important than the final equi l ibr ium values. Appendix B. Numerical details 283 • It is sometimes enlightening to generate bifurcation diagrams for a larger range of parameter values than is of direct interest as there may be hysteresis phenomena or bifurcation points outside the range of interest which affect the behaviour inside this range. Other suggestions which are best described with reference to a particular model are included in the relevant chapters. Clearly, the more one utilises the techniques and packages that have been introduced, the more familiar one becomes with them, and the more creative one can be in their use. The above comments wi l l be most useful to those researchers actively involved in analysing their own models. Appendix C Mathematical details for the sheep-hyrax-lynx model C . l Modelling delays in system dynamics models A system does not always respond immediately to a change in one of its components. Often there is a time lag between the ini t ia l change and the response of the system. For example, a change in hyrax population density wi l l only affect the size of the hyrax population the following season. It takes time for an increase in population density to affect the reproductive success of the hyrax and the survival of their offspring. In order to represent this in the model, averaged or smoothed versions of certain quantities are used in calculating growth rates. Three quantities are averaged in this mode l—hy rax density (Hp), prey abundance (Ap) , and the grazing mult ipl ier (GM)- For example, when sheep have been grazing more than their usual amount their condition is expected to improve resulting in a decline in juvenile mortal ity and a rise in fecundity (Swart and Hearne [116]). However, this wi l l only occur after a prolonged increase in the average amount of pasture consumed. Hence, sheep fecundity and juvenile mortal i ty are functions of the grazing multipl ier average, GM, instead of GM- A first order distributed delay (see MacDonald [77] or May [81]) is used to calculate GM- That is, <1GM _ GM — GM dt tdei where t^ei is the average delay t ime. Examples of this type of delay equation can be found in Forrester [39] and a detailed explanation of the mathematics underlying the above differential equation can be found in MacDonald [77]. The above ordinary differential 284 Appendix C. Mathematical details for the sheep-hyrax-lynx model 285 equation is added to the original system of model equations thus increasing the dimension of the system to be solved. Fortunately this increase in dimension does not pose a problem when solving the system numerically. Further discussions on delays in biological systems and their effects on system be-haviour can be found in [78, 81]. For the sheep-hyrax-lynx model the delays were not found to affect stabil ity even for large average delay times. C .2 Rescaling model equations Suppose the differential equation for the state variable u,- is given by dvi Let Vi = S j U j , then and Vi = — Si dvi 1 dvi dt Si dt = fi(siVi, • • • , SiVi, . . . , SmVm). Si Dropping the bars for convenience we get ^1 - Lf.( \ . — Ji\S^V\, . . . , SiVi, . . . , SmVm I dt Si as required. The state variables have been replaced by s,-u,- and the differential equation for Vi is divided through by s,- as stated in section 3.4. Appendix D Mathematical details for the budworm-forest model D . l Derivation of new foliage equation in spruce budworm model We expect the amount of new foliage consumed by an indiv idual larva to depend on the availabil ity of new foliage per larva, that is, to depend on FbL^F (Fi/Kp represents the amount of foliage that is new foliage). The larvae prefer new foliage to old and thus, in most circumstances, wi l l eat al l the new foliage before moving on to old foliage. Assuming that larval densities do not get so low that there is an overabundance of new foliage, we have new foliage consumption/larva = — ^ — — . Lb However, when larval densities are high, competition among budworm for new foliage becomes significant and larvae eat old foliage more readily than before. We can model this competit ion (see Starfield and Bleloch [113]) by including the factor L d0Lb \ V Fb/KF) where do is the m a x i m u m foliage consumption rate per larva. This factor is close to 1 when there is abundant new foliage per larva (that is, when F b ^ F is large) and close to 0 when there is an overabundance of larvae resulting in intraspecific competit ion. The new equation is n . .. n Fb ( doKpLh\ new foliage consumption/larva = — — — 1 — . KpLb V Fb 286 Appendix D. Mathematical details for the budworm-forest model 287 This factor is rather severe and could lead to negative values for very large larval densities. Instead we can use a negative exponential function [113] which approaches zero when there is very l i tt le new foliage available per larva. This gives Ff, ( d O K F L h new foliage consumption/larva = — — — I 1 — e F>> KFLI, V The total amount of new foliage consumed is then Fb / „ d O K F L b new foliage consumption and this gives KFLh 1 - e Fb remaining new foliage = —r ;T~~(1 ~ e A ) KF KF -A Fb where _ d^KpLb D.2 Summary of model equations The subscript b denotes the in i t ia l or base value for the state variable and subscript e denotes the new value after one year. Foliage (F) F, ^ where Fi = -^[e-A + (KF-l)e-B}Fb, KF Appendix D. Mathematical details for the budworm-forest model A B — doLt KF Cl(A-l + e-A), cj = 0.357. Branch surface area (S) where Budworm (L) where Se rsSi Si = [1 - ds(l - ^)2}Sb. •Tb Le — dsh-rrG2 L& Ks G = H = Le = u = H(2-H), Fi -n. A P m I 1 - 1 , P m . )LS> E = 1 + Er< AsrL\ Athr (EiW1/s - E2)AsrL4, Ei = 165.64 and E2 = 328.52, W — Api (1 + A F 2 ( / ^ F - 1) (1 A F I = 34.1, AF2 = 24.9 and BF = - 3 . 4 , Appendix D. Mathematical details for the budworm-forest mo> L4 = (Ap + Bp-^-)L3, Ap = 0.473 and Bp = 0.828, L3 = e~DL2, Sb(PsatFb2 + LIY kL = 0.425, Li = (1 - qmaxe~C)Lb, C = 0.003L b . Reference values for the parameters are given in table 7.1. Listings Examples of the calling programs for the various models are listed here. The entire calling program does not need to be rewritten each t ime. Once one file has been created, only those lines which define the model equations and set the parameter values and program constants need to be altered. Bazykin model with prey competition—XPPAUT # Mode l equations. dx/dt=a*x-b*x*y/(l+alp*x)-eps*x*x dy/dt=-c*y+d*x*y/(l+alp*x) # Init ial values. x(0)=0.2 y(0)=0.2 # Parameters and nominal values. param a=0.6,b=0.3,c=0.4,d=0.2,alp=0.3,eps=0.001 done 290 System dynamics model—XPPAUT # Sheep, hyrax, lynx and pasture model with pasture l imit ing function. # # Differential equations for the state variables. dpas/dt=(temp*ppn*area*plm(pai)-tssu*gn*gm(pai))/pmax dhj/dt=(hfmax*hf*hfn*hdfm(hda)-hjmax*hj*hjdn*hjdm(hda)-ahjdm(pai) -hjmax*hj*hjmn-hp*hjmax*hj/(hjmax*hj-r-hfmax*hf-|-hmmax*hm))/hjmax dhf/dt=(0.5*hjmax* hj*hjmn-hfmax*hf*hfdn-hfdm(pai)-hfmax*hf*hcn -hp*hfmax*hf/(hjmax*hj-(-hfmax*hf-|-hmmax*hm))/hfmax dhm/dt=(0.5*hjmax*hj*hjmn-hmmax*hm*hmdn-hmdm(pai) -hmmax*hm*hcn -hp*hmmax*hm/(hjmax*hj+hfmax*hf-(-hmmax*hm))/hmmax dlj/dt=(lfmax*lf*lfn*lfm(paa)-ljmax*lj*ljdn*ljdm(paa)-ljmax*lj*ljmn)/ljmax dlf/dt=(0.5*lj*ljmax*ljmn-lfmax*lf*lfdn-lfmax*lf*lcn)/lfmax dlm/dt=(0.5*ljmax*lj*ljmn-lmrnax*lm*lmdn-lmmax*lm*lcn)/lmmax dsj/dt=^(sfmax*sf*sfn*sdfm(gma)-sjmax*sj*sjdn*sjdm(gma)-sjmax*sj*sjmn-max(0,l-lpm(pa))*lut*sjpn-sjmax*sj*sjcn*sjcm)/sjmax dsf/dt=(0.5*sjmax*sj*sjmn-sfmax*sf*sfdn-sfmax*sf*sfcn)/sfmax dsm/dt=(0.5*sjmax*sj*sjmn-smmax*sm*smdn-smmax*sm*smcn*smcm)/smmax dhda/dt=(hd-hda)/del l dpaa/dt=(pa-paa) / del2 dgma/dt=(gm(pai)-gma)/del3 dtr/dt=(-tr - f (mutwool-culling-cssu*ssu+ssu*ssuval)/trmax) / tau # # State variables and ini t ia l values. pas(0)=0.66,hj(0)=0.7,hf(0)=0.525,hm(0)=0.525 lj(0)=0.4,lf(0)=0.6,lm(0)=0.6,sj(0)=0.80670,sf(0)=0.75567,sm(0)=0.50379 hda(0)=1.0,paa(0)=1.0,gma(0)=1.0,tr(0)=3.9069 # # F ixed variables—these are quantities which are used repeatedly # in the model equations. Mak ing them fixed variables simplifies # the appearance of the calculations considerably. hut=hjmax*hjr*hj+hfmax*hf+hmmax*hm lut=ljr*ljmax*lj-|-lmmax*lm+lfmax*lf ssu=sjr*sjmax*sj+sfmax*sf+smmax*sm hd=hut/hun pa=(hut/lut)/(hun/lun) 291 pai=pmax*pas/pav tssu=ssu+max(0,(hut-hun)/hs) hp=lut*lpn*lpm(pa) sjcm=(2.01-argf(1.01,0.01,1.2,gma))/2 smcm=2.01-argf(1.01,0.01,1.2,gma) # # F ixed variables needed to define revenue (mutton and wool sales and cost # of cull ing). mutwool=sjmv*sjmax*sj*sjcn*sjcm+sfmv*sfmax*sf*sfcn-|-smmv*smmax *sm*smcn*smcm-|-sjwv*sjmax*sj+sfwv*sfmax*sf-f-smwv*smmax*sm culling=ccl*(lfmax*lf*lcn+lmmax*lm*lcn)+cch*(hfmax*hf*hcn+hmmax*hm* # # Aux i l ia ry variables—these are quantities other than the state # variables whose values we would like to appear in X P P A U T ' s # data window. Here we have revenue and total revenue. # Total revenue is the quantity referred to in the analysis as revenue, aux rev=mutwool-culling-cssu*ssu aux totrev=mutwool-culling-cssu*ssu+ssu*ssuval # . . . # In order to view the values of the fixed variables in the data window # we need to have the following statements: fhut=hut f lut=lut fssu=ssu fhd=hd fpa=pa fpai=pai ftssu=tssu fhp=hp fsjcm=sjcm fsmcm=smcm # # Parameters and nominal values. param hjr=0.5,ljr=0.5,sjr=0.67,hun=700000,lun=:700,pav=6.6e7 param hs=18,ppn=332,area=200000,gn=365,sub=0.000005,temp=1.0 param hfn=1.5,hjdn=0.5,hfdn=0.1,hmdn=0.1,hjmn=1.0,lpn=84.11,hcn=0.0 param lfn=0.7,l jdn=0.5,lfdn=0.13,lmdn=0.13,ljmn=1.0,lcn=0.0 292 param sfn=0 .75,sjdn=0.1,sfdn=0.02,smdn=0.02,sjmn=0.5 param sjpn=90.0,sjcn=0.09,sfcn=0.28,smcn=0.28 param sjmv=55.0,sfmv=79.0,smmv=75.0,sjwv=5.0,sfwv=7.0,smwv=9.0 param ccl=30.0,cch=1.0,cssu=2.0,ssuval=200.0,pmax=1.0e8 param hjmax=500000,hfmax=500000,hmmax=500000,ljmax=500,lfmax=500 param lmmax=500,sjmax=100000,sfmax=100000,smmax=100000 param dell=0.9,del2=0.6,del3=0.9,tau=0.05,trmax=1.0e7 # # User functions—these are the multipl ier functions. a r g E ( A , B ) = ( A / B ) - l a rgC(A ,B )= ln (a rgE (A ,B )/ (A - l ) ) argr(A,B,slope)=slope*A/((A- l )*argC(A,B)) argf (A,B,s lope,x )=A/( l+argE(A,B)*exp( -argC(A,B) *exp(argr(A,B,slope)*ln(max(x,sub))))) plm(x)=3.25*(x+0.001)*exp(-1.2*(x+0.001)) gm(x)=argf(1.5,0.1,0.8,x) hdfm(x)=2.5-argf(2.4,0.1,1.3,x+0.4) hjdm(x)=argf(1.8,0.1,1.0,x) lpm(x)=argf(1.6,0.8,0.5,x) lfm(x)=argf(1.5,0.01,0.7,x) l jdm(x)=2.0-argf(1.9,0.05,l . l ,x+0.2) sdfm(x)=argf(2.8,0.4,1.3,x) sjdm(x)=12.0-argf(11.5,0.1,2.7,x+1.6) ahjdm(x)=max(0,exp(-3*x)*(hut-hun)*hj*hjmax/(hjmax*hj-|-hfmax*hf-|-hmmax*hm)) hfdm(x)=max(0,exp(-3*x)*(hut-hun)*hfmax*hf/(hjmax*hj+hfmax*hf-t-hmmax*hm)) hmdm(x)=max(0,exp(-3*x)*(hut-hun)*hmmax*hm/(hjmax*hj-|-hfmax*hf-|-hmmax*hm)) # done 293 Ratio-dependent model—XPPAUT # Mode l equations. d M l / d t = g a m l * ( l - e x p ( - o m e g a l / M l ) - M l A b l ) * M l -( l -exp(-omega2*Ml/M2))*phi2*M2/omega2 dM2/dt=gam2*( l -exp( -omega2*Ml/M2) -M2Ab2)*M2-(l -exp(-omega3*M2 / M3)) *phi3*M3/omega3 dM3/dt=gam3*(l -exp( -omega3*M2/M3)-M3Ab3)*M3 # Init ial values. M l (0 )=0 .8 M2(0)=0.2 M3(0)=0.2 # Nondimensionalised parameters and values, param omegal=14.0,omega2=20.0,omega3=16.67 param gaml=0.65,gam2=0.4,gam3=0.4 param phi2=0.07,phi3=0.06 param bl=0.02,b2=0.02,b3=0.0 done 294 Population genetics model I—Interactive A U T O . c Note: lines beginning with c are comments, c c Populat ion genetics model with exponential fitness functions, c P R O G R A M A U T O c I M P L I C I T D O U B L E P R E C I S I O N ( A - H , 0 - Z ) c c N O T E : parameters l iw and lw are V E R Y I M P O R T A N T . They set aside c space for A U T O in the work arrays IW and W . In Interactive A U T O c they must be "hardwired" into the code. If you begin to have problems c with large continuations (e.g. periodic solutions using a big N T S T ) , try c setting l iw and lw larger and recompiling the executable, c P A R A M E T E R ( l i w = 10000) P A R A M E T E R ( l w = 250000) dimension IW(l iw) ,W(lw) . dimension ipar(50),rpar(50),icp(20) character*10 params(20) character*50 name c cal l dfinit c c N D I M (number of state variables) ipar( 1)=2 c IPS (+1 for ode's, -1 for maps) ipar( 2 )= - l c IRS ipar( 3)=0 c I L P ipar( 4)=0 c N T S T ipar( 5)=15 c N C O L ipar( 6)=4 c I A D 295 ipar( 7 =3 ISP ipar( 8^ =1 I S W ipar( 9^ =1 I P L T ipar(10 )=0 N B C i p a r ( l T =0 N I N T ipar(12} =0 I A D S ipar(13 =1 N M X ipar(14^ =100 N U Z R ipar(15^ =0 N P R ipar(16} =50 M X B F ipar(17) =5 IID ipar(18) =2 I T M X ipar(19) =8 I T N W ipar(20) =5 N W T N ipar(21) =3 J A C ipar(22) =0 ICP( i ) icp( l ) = 1 icp(2) = 2 D O 1 1= =1,2 ipar(30+i)=icp( 1 C O N T I N U E DS rpar( l )= :0.0001d0 D S M I N rpar(2)= :0.000002d0 D S M A X rpar(3)= :0.01d0 RLO rpar(4)= =0.0 R L 1 rpar(5) = 4 .0 AO rpar(6) = -10.0 A l rpar(7) = 250.0 E P S S rpar(8) = = l .d -6 E P S L ( i ) , i =1,20 rpar(9) = = l .d -6 E P S U rpar(10) = l .d -6 nparams=number of parameters that you want to vary nparams = 6 declaration of parameter names params(l ) = ' a l l ' params(2) = ' a l 2 ' params(3) = 'a22' params(4) = ' b l l ' params(5) = ' b l 2 ' params(6) = 'b22' name = ' Genetics model ' call autool(ipar,rpar,iw,liw,w,lw,params,nparams,name) stop end 297 c S U B R O U T I N E VECFLD(ndim,u , icp ,par , i jac , f , t ) c c This subroutine evaluates the right hand side of the first order system c c input parameters : c nd im - dimension of u and f. c u - vector containing u. c par - array of parameters in the differential equations. c icp - par( icp(l ) ) is the in i t ia l 'free' parameter. c par(icp(2)) is a secondary 'free' parameter, c for subsequent 2-parameter continuations. c ijac - =1 if the jacobians dfdu and dfdp are to be returned, c =0 if only f(u,par) is to be returned in this cal l . c t - current t ime. c c value to be returned : c f - f(u,par) the right hand side of the ode. c impl ic i t double precision (a-h,o-z) c dimension u(ndim),par(30) dimension f(ndim) c Parameters a l l = par ( l ) a l 2 = par(2) a22 = par(3) b l l = par(4) b l 2 = par(5) b22 = par(6) c State variables P = u ( l ) eN = u(2) c Fitness functions w l l = e x p ( a l l - b l l * e N ) w l 2 = exp(al2-bl2*eN) w22 = exp(a22-b22*eN) 298 c Mean fitnesses w l m a r g = P * w l l + ( 1 - P ) * w l 2 w2marg = P*wl2+(1-P)*w22 fitmean = P*wlmarg - f ( l -P )*w2marg c c D I F F E R E N C E E Q U A T I O N S c f ( l ) = P*wlmarg/fitmean f(2)= fitmean*eN c return end c S U B R O U T I N E PARDER(ndim,u, icp ,par , i jac ,dfdu, t ) c c this subroutine evaluates the derivatives c of the first order system and with respect to (u(l),u(2)). c Not included for this model hence ijac=0 in the first subroutine. c return end c S U B R O U T I N E DFDPAR(ndim,u, icp ,par , i jac ,dfdp, t ) c c this subroutine evaluates the derivatives c of the first order system and with respect to free parameters. c Not included for this model. c return end c S U B R O U T I N E STPNT(nd im,u ,par ) c c in this subroutine the steady state starting point must be defined, c (used when not restarting from a previously computed solution), c the problem parameters (par) may be init ial ized here or else in init . c 299 c n d i m - dimension of the system of equations, c u - vector of dimension ndim. c upon return u should contain a steady state solution c corresponding to the values assigned to par. c par - array of parameters in the differential equations. c impl ic i t double precision (a-h,o-z) c dimension u(ndim),par(30) c c init ial ize the problem parameters. par( l )=2.1d0 par(2)=1.9d0 par (3 )= l . ld0 par(4)=1.0d0 par(5)=0.904d0 par(6)=0.524d0 par(14) = D B L E ( 1 ) c init ial ize the steady state. u( l )=0.5d0 u(2)=2.1008d0 c return end c S U B R O U T I N E S P R O J (ndim, u, isw, icp, par, vaxis, pt) c c this subroutine can be used to define a special projection c in the bifurcation window. This subroutine is called c when the ' S P ' is toggled O N (issue the command successively c to turn the toggle from O N to O F F , and vice versa). c c input values: c nd im - dimension of u. c u - vector containing coordinates of current solution. c isw - the number of parameters being used in the current c continuation. 300 icp - par( icp(l ) ) is the in i t ia l 'free' parameter. par(icp(2)) is a secondary 'free' parameter, for subsequent 2-parameter continuations. par - array of parameters in the differential equations. vaxis - controlled by program constant I P L T (see A U T 0 8 6 User Manual) , this is the second number per line written in unit 7 (file fort.7). return values: pt - array whose 1st, 2nd and 3rd elements are plotted the x, y and z axes, respectively. impl ic i t double precision (a-h,o-z) dimension icp(20) dimension u(ndim), par(30), pt(3) if ( isw.eq.l) then p t ( l ) = par( 1) pt(2) = u(2) pt(3) = vaxis else if (isw.eq.2) then pt ( l ) = par( 1 ) pt(2) = par( 2 ) pt(3) = vaxis endif return end S U B R O U T I N E S P J A X S (ndim, isw, icp, axes ) this subroutine defines the names of the axes used in the projection defined in the subroutine sproj. input value: nd im - dimension of u. 301 c IABS( isw) - the number of parameters being used in the current c continuation. c icp - par( icp(l ) ) is the in i t ia l 'free' parameter, c par(icp(2)) is a secondary 'free' parameter, c for subsequent 2-parameter continuations, c return value: c axes - character string array with the x, y, and c z axes names, respectively. c integer*4 nd im, isw, icp(20) character* 10 axes(3) c if (isw.eq.l) then axes(l) = ' icp( l ) ' axes(2) = ' N ' axes(3) = ' ' else if (isw.eq.2) then axes(l) = ' icp( l ) ' axes(2) = 'icp(2)' axes(3) = ' ' endif c return end c c*** graphics initializations for interactive A U T O c S U B R O U T I N E G P H D F T ( ldebug, l intog, labtog, ldsplt, @ leigen, lfltog, lsavpt, @ lgraph, lvideo, lsproj, nproj , @ ndmplt , delay, sclbif, scldis, @ sclev, filext ) c logical ldebug, l intog, labtog, ldsplt, leigen, lfltog logical lsavpt, lgraph, lvideo, lsproj integer*4 nproj , ndmplt real*8 delay 302 real*8 sclbif(6), scldis(6), sclev(6) character*10 filext(2) c c*** toggles c ldebug: T => debugging output c l intog: T => plotting of lines between points c labtog: T => plots two-character identifier at bifurcation points c ldsplt : T => open optional U N I X Graph window c leigen: T => open eigenvalue plott ing window c lfltog: T => il luminates points temporarily as they are plotted c lgraph: T => use graphics (F => can then run jobs in the background) c lvideo: T => reposition windows in botton 1/4 of screen for videotaping c lsproj: T => plot special projection denned in subroutine sproj c in the bifurcation window c lsavpt: T => eigenvalues saved in fort. 11 ('svaut *' moves this to m.*) c D O N O T A L T E R the following lines. ldebug = .false. l intog = .false. labtog = .true. ldsplt = .false. leigen = .true. lfltog = .false. , lcomfl = .false. loutfl = .false. lgraph = .true. lvideo = .false. lsproj = .false. lsavpt = .false. c c*** default window scales c ... Bifurcation window scales c eg. when isw = 1 => default plot axes x,y,z = par( icp( l ) ) , u ( l ) , u(2); c when isw = 2 => default plot axes x,y,z = par( icp( l ) ) , par(icp(2)), u ( l ) c The following lines can also be changed interactively. sclbif ( l ) = 0.5d0 sclbif(2) = 0.56d0 sclbif(3) = O.OdO 303 sclbif(4) = l.OdO sclbif(5) = O.OdO sclbif(6) = lO.OdO c ... Eigenvalue window scales sclev(l) = -2.0d0 sclev(2) = 2.0d0 sclev(3) = -2.0d0 sclev(4) = 2.0d0 c c*** set files names c comfil = 'input ' outfil = 'output ' c c*** file strings for saving, deleting, etc. c fTlext(l) = 'gen ' filext(2) = 'gen2 ' c c*** other stuff c nproj = number of projections to be plotted (up to nine) c ndmplt = dimension of the bifurcation window plot c delay = factor for duration of flash display c nproj = 2 ndmplt = 3 delay = O.OdOO c return end 304 Population genetics model I—DSTOOL. # include