I introduced a recent paper, in the previous post, on the Tragedy of the Commons. The paper is published in *Sustainability*, an open access paper and is therefore available to everyone. One bit that I did not include in the paper, however, is the code that I used for the simulations therein. On the suggestion of a friend, I’ll publish them here in a series of posts. (Thanks Mauro)

The simulations are very simple, being based on difference equations, and were all coded in Octave. The code is not pretty, but I cleaned it up and added some comments. If you have questions, just drop me a comment on the blog. This first simulation is of a basic tragedy involving two users, and corresponds to Figures 1B-C in the paper. A resource is simulated to a steady state based on a Ricker model, for 100 steps, and then users begin to utilize it. They increase their utilization per time step at a steady rate (acceleration of utilization is constant), and the simulation is run until the resource is exhausted. The basic output of resource, one user’s benefit, and total utilization, are plotted in the figure here. The equations simulated are outlined in the excerpt from pg. 754 of the paper:

“*We now introduce a term to represent consumption or utilization by our human TOC agents.
where U(t) is the total fraction of R utilized by human users at time t, and ranges from 0 to 1. The term in square brackets on the right hand side of the equation is the modified growth rate of R. The resource is stable when this term is positive or zero, that is, the rate of resource renewal exceeds or is equal to utilization. Since U is the total standardized utilization rate of all users, it may be expanded to
where there are N users and is the standardized utilization rate of the user at time t. The benefit gained from utilization is modelled as
where is the benefit to user i, f is a factor that converts resources utilized into another commodity, for example converting harvested food to energy or currency, and B is the total benefit to all users.*”

The code follows, and can be executed by running it in an interactive Octave session. Please be aware that WordPress has somewhat limited options for formatting code. I’ve used Matlab formatting here, but wrapped lines might not be obvious. Single lines are all terminated with a semicolon, “;”, so just look out for those.

#RESOURCE DATA #initial resource level R0 = 1; R = zeros(1,2); R(1, 1) = 1; R(1, 2) = 1; #carrying capacity K = 50; #intrinsic growth rate r r = 1.5; #USER DATA #There are 2 users #initial acceleration of utilization du0 = 1.05; #individual user parameters, users 1 and 2 #initialize arrays u1 = zeros(1,2); u2 = zeros(1,2); b1 = zeros(1,2); b2 = zeros(1,2); u1(1,1) = 1; #initial utilization rate u1(1,2) = .001; #conversion factor of resource to benefit f1 = 1; b1(1,1) = 1; b1(1,2) = 0; u2(1,1) = 1; #initial utilization rate u2(1,2) = .001; #conversion factor of resource to benefit f2 = 1; b2(1,1) = 1; b2(1,2) = 0; U = zeros(1,2); U(1,1) = 1; #total utilization U(1,2) = b1(1,2) + b2(1,2); du = zeros(1,2); du(1,1) = 0; du(1,2) = du0; #BEGIN SIMULATION #RUN RESOURCE TO STEADY STATE, NO UTILIZATION for count2 = 2:100 R(count2, 1) = count2; R(count2, 2) = R(count2-1, 2) * (exp(r*(1-(R(count2-1, 2)/K)))); u1(count2, 1) = count2; u1(count2, 2) = u1(1,2); u2(count2, 1) = count2; u2(count2, 2) = u2(1,2); du(count2,1) = count2; du(count2,2) = du0; endfor #INITIATE UTILIZATION for count1 = 101:265 R(count1, 1) = count1; #modified Ricker model with utilization R(count1, 2) = R(count1-1, 2) * (exp(r*(1-(R(count1-1, 2)/K))) - (u1(count1-1,2)+u2(count1-1,2))); u1(count1, 1) = count1; #increase utilization u1(count1, 2) = u1(count1-1, 2) * du(count1-1,2); #benefits (equal for both users in this case) b1(count1, 1) = count1; b1(count1, 2) = f1 * u1(count1-1, 2) * R(count1-1, 2); b2(count1, 1) = count1; b2(count1, 2) = f2 * u2(count1-1, 2) * R(count1-1, 2); u2(count1, 1) = count1; u2(count1, 2) = u2(count1-1, 2) * du(count1-1,2); #total benefits U(count1, 1) = count1; U(count1, 2) = b1(count1, 2) + b2(count1, 2); du(count1,1) = count1; du(count1,2) = du(count1-1,2); endfor