Fish Harvesting Revisited

Kurt Bryan and SIMIODE

Notebook for fish logistic-growth harvesting project in Section 3.5.2.

The Data: The population data from Table 3.10, starting with 1978, which we will call year 0.

udata = {72148, 73793, 74082, 92912, 82323, 59073, 59920, 48789, 70638, 67462, 68702, 61191, 49599, 46266, 34877, 28827, 21980, 17463, 18057, 22681, 20196, 25776, 23796, 19240, 16495, 12167, 21104, 18871, 21241, 22962}

However, Mathematica indexes from 1 to 30, so udata[[1]] = 72148.

The harvesting data is

hd = {0.18847, 0.149741, 0.21921, 0.17678, 0.28203, 0.34528, 0.20655, 0.33819, 0.14724, 0.19757, 0.23154, 0.20860, 0.33565, 0.29534, 0.33185, 0.35039, 0.28270, 0.19928, 0.18781, 0.19357, 0.18953, 0.17011, 0.15660, 0.28179, 0.25287, 0.25542, 0.08103, 0.087397, 0.081952, 0.10518}

A plot of each data set:

ListPlot[udata]

ListPlot[hd]

Fitting the Data to the ODE: Form finite difference approximations (u(k+1) - u(k))/1 from the population data, for k = 1 to k = 29 (= n-1). This approximates u'(t) when t = k:

udiffdata = Table[udata[[k+1]]-udata[[k]],{k,1,29}]

Next substitute u = udata[[k]] into the harvested logistic ODE u'(t) = r*u(t)*(1-u(t)/K) - h(t)*u(t) for k = 0 to k = n-1, to approximate the right side of the ODE at time t = k for k = 0 to k = n-1. This expression depends on r and K.

udiffdata2 = Table[r*udata[[k]]*(1-udata[[k]]/K)-udata[[k]]*hd[[k]],{k,1,29}];

Form a sum of squares that depends on r and K

SS = Sum[(udiffdata[[k]]-udiffdata2[[k]])^2,{k,1,29}];

Finding the optimal r and K: Start with a plot of SS(r, K), or log(SS).

Plot3D[Log[SS],{K,40000,300000},{r,0.1,0.5}]

Or perhaps a contour plot:

ContourPlot[Log[SS],{K,40000,300000},{r,0.1,0.5},Contours->50]

Something near r = 0.3, K = 200000 might be a good initial guess at a minimizer for SS.

If we use these values for r and K (but you should find the true minimizer) we could compare the predicted cod population to the data by following the suggestion of Modeling Exercise 5.2.3. Let U be an array to hold our numerical estimates of the cod population:

U = Table[0,{k,1,30}];

Now use the harvested logistic ODE with our estimates for r and K to march the predicted cod biomass out in time (this is essentially Euler's method):

U[[1]] = udata[[1]]
For[k=1,k≤29,k++,
U[[k+1]] =(U[[k]] + r*U[[k]]*(1 - U[[k]]/K) - hd[[k]]*U[[k]])/.{r->0.3,K->200000}
]

p1 = ListPlot[udata];
p2 = ListPlot[U,PlotStyle->{Red}];
Show[p1,p2]