source code

/* density dependent logistic growth simulation */

#include	<stdlib.h>  /* for drand48() */
#include	<stdio.h>

#define	MAX_TIME 	1000
#define	PRNT_INTERVAL	10
#define	DATA_POINTS		(MAX_TIME/PRNT_INTERVAL)
#define	BIRTH_RATE	0.2
#define	DEATH_RATE	0.1
#define	KK		20.0  /* This is NOT carrying capacity */
#define INIT_POP_SIZE   5

/*
  Model 1 ---   dN/dt = b N (1 - N / k) - d N
  Model 2 ---   dN/dt = b N - d N (1 + b N / (k d))
*/

double	birthR, deathR, kk;

int	PopChange(int, int);

int main(void)
{
  int	   run, ttt, n1, n2, n3, outputCnter, seed, numRuns;
  double   pops1[DATA_POINTS], pops2[DATA_POINTS];
  
  birthR = BIRTH_RATE;
  deathR = DEATH_RATE;
  kk = KK;

  seed = 1456739853;
  srand48(seed);
  
  fprintf(stderr, "Enter number of runs:\n");
  scanf("%d", &numRuns);

  /* initialize the array */
  for(ttt=0;ttt<DATA_POINTS;ttt++)
    {
      pops1[ttt] = pops2[ttt] = 0.0;
    }
  
  /* run simulations */
  for(run=0;run<numRuns;run++)
    {
      outputCnter = 0;
      n1 = n2 = INIT_POP_SIZE;
      
      for(ttt=0; ttt<MAX_TIME; ttt++)
	{
	  n1 += PopChange(n1,1);
	  n2 += PopChange(n2,2);
	  
	  if(ttt%PRNT_INTERVAL == 0)
	    {
	      pops1[outputCnter] += n1;
	      pops2[outputCnter] += n2;
	      outputCnter++;
	    }
	  
	  if(n1==0 && n2==0)  /* extinction */
	    break;
	}
    }
  
  /* print out the results */
  printf("time  Ave(n1)  Ave(n2)\n");
  for(ttt=0; ttt<DATA_POINTS; ttt++)
    printf("%d  %lf  %lf \n",ttt*PRNT_INTERVAL,
	   pops1[ttt]/(double)numRuns,
	   pops2[ttt]/(double)numRuns);
  
  return 0;
}

/*
   This function calculate the changes in the population size
   The 2nd argument is model number (1 for d.d. birth, 2 for d.d. death)
*/
int PopChange(int curPopSize, int model)
{
  int popIncrement, newborn=0, death=0, i;
  double popD;
  
  popD = (double) curPopSize;

  /* model 1 : density dep birthR */
  if(model==1)
    {
      for(i=0; i<curPopSize; i++)
	{
	  if(drand48() < birthR*(1.0 - popD/kk))
	    newborn ++;
	}
      for(i=0;i<curPopSize;i++)
	{
	  if(drand48() < deathR)
	    death ++;
	}
    }
  
  /* model 2 : density dep death */
  if(model==2)
    {
      for(i=0; i<curPopSize; i++)
	{
	  if(drand48() < birthR)
	    newborn++;
	}
      for(i=0; i<curPopSize; i++)
	{
	  if(drand48() < deathR + birthR*popD/kk)
	    death++;
	}
    }

  popIncrement = newborn - death;
  return popIncrement;
}