Solution sketch, difference equation

3.) Computation: altering the fishery model

Now let's introduce something in the model which has been left out: Prices. We use prices as an informational connection between the fish stock, harvest and effort.

We can work with a simple inverse demand function of the form:

$$P(t) = a-cY(t)$$

And we can also assume that effort is positively associated with price in the form:

$$\frac{df(t)}{dt} = g(P(t)- \bar{P})$$

If you want to go a step further, you can also make the sensible assumption that effort is costly and define a profit function:

$$\pi(t) = (P(t) - wf(t))Y(t)$$

where w is the cost of effort (you can think of it as a wage to the fisher-people).

Create a system of differential equations with some or all of these elements. You may need to experiment with different parameters to get simulations that appear reasonable.

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp

As a reminder, here was our initial system where effort was dynamically related to effort.

$$\frac{dB(t)}{dt} = rB(t)(1-\frac{B(t)}{k})-qf(t)B(t)$$$$\frac{df(t)}{dt} = \gamma (\omega - \frac{ B(t)}{k})$$

Where we have a harvesting function related to effort:

$$Y(t)=qf(t)B(t)$$

We can now rewrite our system with fish stock, B, and effort as our dynamic variables, but now with harvest and price as intermediate variables:

$$\frac{dB(t)}{dt} = rB(t)(1-\frac{B(t)}{k})-qf(t)B(t)$$$$\frac{df(t)}{dt} = g(P(t)- P_{bar})$$$$Y(t)=qf(t)B(t)$$$$P(t) = a-cY(t)$$$$P(t) = a-cqf(t)B(t)$$$$\frac{df(t)}{dt} = g([a-cqf(t)B(t)]- P_{bar})$$$$\frac{df(t)}{dt} = gB(t)(\frac{a}{B(t)}-cqf(t)- \frac{P_{bar}}{B(t)})$$$$\frac{df(t)}{dt} = gB(t)(\frac{a-\bar{P}}{B(t)}-cqf(t))$$

def dS2(t,S):      
    return(np.dot(np.array([[r*(1-S[0]/k), -q*S[0]],[g*((a-P_bar)/S[0]-c*q*S[1]),0]]), S))

t_eval = np.arange(0, 200, .1)
r=.10
q=.04
k=1000
S0 = [100,1] #initial values of B and P
gamma = .005
omega = .4
P_bar = 1

a = 1
c = .01
g=.05 #relationship between effort and price

fishSol3 = solve_ivp(dS2, [0, 200], S0, t_eval=t_eval)

Y = q*fishSol3.y[0]*fishSol3.y[1]
P = a-c*Y

fig, ax = plt.subplots(2)
ax[0].plot(fishSol3.t, fishSol3.y[0])
ax[0].set_xlabel('t')
ax[0].set_ylabel('B(t)')
ax[1].plot(fishSol3.t, fishSol3.y[1])
ax[1].set_xlabel('t')
ax[1].set_ylabel('f(t)')

Text(0, 0.5, 'f(t)')

fig, ax = plt.subplots(2)
ax[0].plot(fishSol3.t, Y)
ax[0].set_xlabel('t')
ax[0].set_ylabel('Y(t)')
ax[1].plot(fishSol3.t, P)
ax[1].set_xlabel('t')
ax[1].set_ylabel('P(t)')

Text(0, 0.5, 'P(t)')