Logistic Population Growth -------- ---------- ------ Logistic population growth is given by the equation dN/dt = rN(1 - N/K) where t is the time N is the current population (a function of t) r is the population growth rate, a constant K is the carrying capacity of the environment, a constant (we allow this to be non-integral) The population starts from an initial value N(0) at t = 0 and grows or shrinks exponentially until it approximately equals K. The above equation has the solution N(t) = K / ( 1 + F * exp ( -rt) ) ) where F = ( K - N(0) ) / N(0) These equations define what is called a `deterministic model'. Real population growth curves tend to wander from deterministic model predictions unless N is large, as we will show in this problem. Population growth can be more carefully modeled by a `stochastic model'. To build this model we define B(N) = max (0, (b0 - b1 * N) * N) Birth rate D(N) = (d0 + d1 * N) * N Death rate p(N,t) probability that the population is N at time t b0, b1, d0, d1 >= 0 are constants of the model From this we get the stochastic differential equation dp(N,t)/dt = - p(N,t)*(B(N) + D(N)) + p(N+1,t)*D(N+1) + p(N-1,t)*B(N-1) Here the corresponding deterministic model is dN/dt = B(N) - D(N) = ((b0 - d0) - (b1 + d1) * N) * N so comparing to the above we get r = b0 - d0 r/K = b1 + d1 hence K = (b0 - d0) / ( b1 + d1 ) For the stochastic model all four constants b0, b1, d0, d1 are needed to define the model, but for the determin- istic model only r and K are needed. The deterministic initial condition N = N(0) is equiva- lent to the stochastic initial condition p(N,0) = 1 and p(n,0) = 0 for n != N. Integrating the stochastic differential equation to find p(N,t) turns out to be difficult, even numerically using a computer. But there is a steady state where all the dp(N,t)/dt = 0, and we can solve for the p(N,t) in this steady state. Call these steady state probabilities P(N), and then replacing dp(N,t)/dt by 0 and p(N,t) by P(N) in the above equation we get 0 = - P(N)*(B(N) + D(N)) + P(N+1)*D(N+1) + P(N-1)*B(N-1) This equation in turn can be rewritten as P(N+1)*D(N+1) - P(N)*B(N) = P(N)*D(N) - P(N-1)*B(N-1) = G a constant independent of N For N = 0 this equation is P(0)*D(0) - P(-1)*B(-1) = 0 and if D(0) = P(-1) = 0 we have G = 0. Thus P(N+1) = P(N)*B(N)/D(N+1) if D(N+1) != 0 We have B(0) = 0 and D(N) > 0 for N >= 1 and this gives the solution P(0) = 1, P(N) = 0 for N > 0, which is the `extinction solution' and is not very interes- ting. However, if we assume that extinction never actually occurs, this is equivalent to assuming that P(0) = 0 and throwing out the equation P(0)*D(0) - P(-1)*B(-1) = G Then we have P(1)*D(1) - P(0)*B(0) = G and as P(0)= 0, G = P(1)D(1). We thus get P(N+1) = ( P(N)*B(N) + P(1)*D(1) )/D(N+1) and if we set Q(1) = 1 Q(N+1) = ( Q(N)*B(N) + Q(1)*D(1) )/D(N+1) then P(N) = Q(N)/(sum Q(N) for all N >= 1) which normalizes P(N) so the sum of the probabilities is 1. Note that to actually sum the Q(N) in a computer you need to stop summing at some finite value of N. In our case B(N) = 0 when N >= b0/b1, and sum ( Q(1)/D(N) for N > L ) <= Q(1)/(d1*L) (because (sum 1/N**2 for N > L) is <= 1/L), so if we set L = max (b0/b1, 1/d1)) the sum of Q(N) for N > L will be at most Q(1), which we expect to be very small part of the total sum, so we restrict the summing to 1 <= N <= L and our equation becomes P(N) = Q(N)/(sum Q(N) for 1 <= N <= L ) Given this we define the steady state statistics of the model as MEAN = steady state mean of N = sum( N * P(N) for 1 <= N <= L ) VAR = steady state variance of N = sum( (N-MEAN)**2 * P(N) for 1 <= N <= L ) STD = steady state standard deviation of N = sqrt ( VAR ) Lastly a simulation of the stochastic model can be implemented by the following pseudo-code N = N(0) t = 0 loop: choose time s to next event choose whether next event is birth or death if next event is birth: N = N + 1 else if next event is death: N = N - 1 t = t + s Events occur at the rate R(N) = B(N)+D(N) so s is an exponentially distributed random variable such that probability {s' => 0 : s' <= s} = exp ( - R(N) * s ). Note that this equals probability {s' => 0 : exp(-R(N)*s') >= exp(-R(N)*s} so if we set Y = exp ( - R(N) * s ) we get probability {s' => 0 : exp(-R(N)*s') >= Y} = Y and Y = exp(-R(N)*s) is therefore uniformly distributed. So we can choose s by to choose s: pick a pseudo-random uniformly distributed number Y, 0 <= Y <= 1. set s = - (ln Y)/R(N) The relative probabilities of births and deaths are B(N) and D(N) so to choose whether the next event is a birth or a death: pick a pseudo-random uniformly distributed number Y, 0 <= Y <= 1. if Y <= B(N)/(B(N)+D(N)) the event is a birth otherwise it is a death If N == 0, then B(0) = D(0) = R(0) = 0, and this is a special case in which N is stuck at 0 forever. By now you must have guessed that you are going to be asked to compute all the above. Furthermore, you must get very precisely the same answers as the judge. To do this you need to use double precision floating point numbers and the following pseudo-random number genera- tor: C or C++: long long seed; double random_Y ( void ) { seed = 16807 * seed; seed = seed % 2147483647 return double(seed) / 2147483646; } JAVA: long seed; double random_Y ( void ) { seed = 16807 * seed; seed = seed % 2147483647 return double(seed) / 2147483646; } You will be given an initial value of seed, and for each pseudo-random number Y you need (including the first), you call random_Y(). Input ----- For each of several test case, first a line containing just the name of the test case. Then a line containing b0 b1 d0 d1 N(0) Tsize Nsize Tinc Ninc where b0 .. N(0) define the model, Tsize is the number of lines in the plot to be produced below, Nsize is the number of columns in each of these lines, Tinc is the amount time is incremented between these lines, and Ninc is the amount N is incremented between columns of these lines. b0, b1, d0, d1, Tinc, Ninc are floating point N(0), Tsize, Nsize are integers 0 <= b0, b1, d0, d1 0 < N(0) 0 < Tsize <= 100 0 < Nsize <= 80 0 < Tinc 0 < Ninc The lines of a test case between the second line of the test case and the last line of the test case each have the form C seed where C is the display character for plotting (see `Output' below) and seed is a pseudo-random number generator seed: 0 < seed < 2147483647 ( == 2**31-1 ) The last line of a test case contains just `.'. Input ends with an end of file. Output ------ For each test case, first output an exact copy of the first test case input line which names the test case. Then output a plot containing Tsize lines each with Nsize columns. The T+1'st line corresponds to the time t = T*Tinc (so the first line corresponds to t = 0). To plot a number x with a display character C on a line, put the character C in column round ( x / Ninc ) + 1 where `round' rounds to the nearest integer. For each seed you are to plot the simulation using that seed with the display character given on the same input line as the seed. Note that display characters from a simulation may overwrite display characters from a previous simulation. Then you are to plot the deterministic model N(t) using the display character `*'. Note that this display character may overwrite simulation display characters. There should be NO TABs in any plot line. You may end a plot line with single space characters. After the plot output the line: r = #, K = #, MEAN = #, STD = # where the #'s are as follows: r = b0 - d0 is the rate of population growth when N is small K = (b0 - d0) / ( b1 + d1 ) is the carrying capacity MEAN = mean of the stochastic steady state as computed above STD = standard deviation of the stochastic steady state as computed above Print all the #'s to at least 2 decimal places. Note the spacing required in this line: `='s are sur- rounded by whitespace and `,'s are followed by white- space but NOT preceded by whitespace. Also, you may NOT use TABs in this line. Failure to observe these rules may result in a `format error' score. Note that as decisions requiring comparison of floating point values are made, output in general will be sensi- tive to floating point accuracy. However, the judge has tuned the judging input so if you use double prec- ision and follow the above formulae exactly, you will get exactly the plot the judge gets, and you are in fact required to do so. One thing to be careful of is the order in which you use the values returned by random_Y(); specifically you must choose s BEFORE you choose whether the event is a birth or death. Note that a `format error' score might mean that you have plotted display characters in the wrong columns but have somehow managed to get the right display character overlays, so you may consistently be a column off. Notes ----- The simulations reveal that after a population reaches capacity it wanders enough that it does not appear stable. This is because the standard deviation of the stochastically stable solution is not that small. Therefore deterministic stability != apparent stochastic stability Sample Input ------ ----- -- SAMPLE 1 -- 2.2 0.2 0.1 0.1 1 15 50 0.5 0.25 x 838765873 # 098763498 @ 162738493 . -- SAMPLE 2 -- 10.1 0.1 0.0 0.1 10 15 50 0.10 1.2 x 898765873 # 098763498 @ 62738493 . Sample Output ------ ------ -- SAMPLE 1 -- * x* @ * @ * x @ * @ # # *@ x x * # * x # @ * # # * @ x # * @ # * @ x # * @ x # x * @ # @ * r = 2.10, K = 7.00, MEAN = 6.54, STD = 1.75 -- SAMPLE 2 -- * x * # x @ * # x @ * x @ * x # * @ x#*@ # * @ x *@ # @ * # x* @ x * # x *@# x @ * # *x @ # r = 10.10, K = 50.50, MEAN = 49.99, STD = 5.05 File: logistic.txt Author: Bob Walton Date: Wed Oct 14 20:50:13 EDT 2009 The authors have placed this file in the public domain; they make no warranty and accept no liability for this file. RCS Info (may not be true date or author): $Author: walton $ $Date: 2009/10/15 00:51:52 $ $RCSfile: logistic.txt,v $ $Revision: 1.16 $