# Lab 3
## Introduction
In this lab we will analyse population dynamics under the logisitic model with managed harvesting.

First import the modules we need.


```python
import seaborn as sns
import matplotlib.pyplot as plt
import pandas as pd
from numpy import meshgrid, linspace, sqrt, arange
from scipy.integrate import odeint
```

## Harvesting of fish
A population of fish in a lake, left to its own devices, is modelled by the logistic differential equation
\begin{align}
\frac{\mathrm{d}y}{\mathrm{d}{t}} = 4y(1-y),
\end{align}
where the population $y$ is in units of thousands of fish and time $t$ is measured in years.

First define a function for $\mathrm{d}y/\mathrm{d}x$ in terms of $y$ and $x$.


```python
def diff_eq(y, x):
    return 4 * y * (1 - y)
```

Next define a function that creates a Plotly Figure object that contains a slope field and, optionally, a few solutions to initial value problems.

It automates a few things we did in the last lab.

- `diff_eq` is the differential equation to be plotted
- `x` and `y` should be outputs from `meshgrid`. 
- `args` is any additional arguments to `diff_eq` (we will use that below).
- `initial_values` is a list (or array) of starting $y$ values from which approximate solutions will start. The corresponding $x$ value is the minimum element of `x`.

Note that the numerical solutions will plotted for the whole range of $x$ values in `x`, so if they blow up you will probably get a warning and less-than-useful plot.


```python
def create_slope_field(diff_eq, x, y, args=(), initial_values=()):    
    S = diff_eq(y, x, *args)
    L = sqrt(1 + S**2)
    fig, ax = plt.subplots(figsize=(5, 5))
    q = ax.quiver(x, y, 1/L, S/L, scale=25, headwidth=0, headlength=0, color='grey')
    x = linspace(x.min(), x.max())
    data = {'x': x}
    for y0 in initial_values:
        y = odeint(diff_eq, y0, x, args)[:,0]
        data[f'y({x[0]}) = {y0}'] = y
    data = pd.DataFrame(data)
    data = data.melt(id_vars=['x'], value_name='y', var_name='initial value')
    sns.lineplot(data=data, x='x', y='y', hue='initial value')
    return fig, ax
```

The slope field below should hopefully give you some idea for the fish population dynamics.

Note that we use `arange` rather than `linspace` this week so that we can carefully control the increments between our grid points. `arange(0, 1.1, 0.25)` returns an array that starts with 0 and increments by 0.25 until it exceeds 1.1.

The plot also contains the solution curves for 
$y(0) = 1$ and $y(0) = 0.4$. Edit the cell to also include the solution curve for $y(0)=1.4$.


```python
x, y = meshgrid(arange(0, 1.1, 0.05), arange(-0.4, 1.41, 0.05))
fig, ax = create_slope_field(diff_eq, x, y, initial_values=(0.4, 1))
```

### Equilibrium solutions
Looking back to our differential equation, $\mathrm{d}y/\mathrm{d}t = 0$ when $y(t) = 0$ or $y(t) = 1$. Looking at the slope field, we see that the equilibrium solution $y(t) = 1$ is stable (this is the carrying capacity here, corresponding to 1000 fish), whereas the equilibrium solution $y(t) = 0$ is unstable. Any non-zero initial population will eventually stabilise at 1000 fish.

### What will happen if harvesting is now commenced at a steady rate?
For the simplest harvesting model, assume that $H$ units (thousands) of fish are taken
continuously (smoothly) over the year, rather than at one instant each year.
Note that the units of $H$ are the same as those of $\mathrm{d}y/\mathrm{d}t$, thousands of fish per year, so we simply subtract $H$ from the RHS of our existing equation to give the DE with harvesting as
\begin{align}
\frac{\mathrm{d}y}{\mathrm{d}{t}} = 4y(1-y) - H.
\end{align}
Again, the (constant) equilibrium solutions are found by setting $\mathrm{d}y/\mathrm{d}t = 0$, giving from the quadratic formula (check this),
\begin{align}
y(t) = \frac{4\pm\sqrt{16-16H}}{8} = \frac{1\pm\sqrt{1-H}}{2}.
\end{align}
What happens after harvesting starts will depend on the equilibrium solutions, their
stability and the initial number of fish $y(0)$.

Start by redefining `diff_eq` to include the `H` parameter. Note that defining `diff_eq` again overides our original definition.


```python
def diff_eq(y, x, H=0):
    return 4 * y * (1 - y) - H
```

Now set $H = 0.6$ and plot the slope field. This is done by setting `args=(0.6,)` when we call `create_slope_field`. This is exactly how you would pass additional arguments like this one to `odeint` if you were calling it directly.


```python
x, y = meshgrid(arange(0, 1.1, 0.05), arange(-0.4, 1.41, 0.05))
fig, ax = create_slope_field(diff_eq, x, y, args=(0.6,))
```

From the solutions to the quadratic equation above, the equilibrium solutions of the DE are found to be $y(t) \approx 0.184$ and $y(t) \approx 0.816$. The previous equilibrium solution with no harvesting at $y(t) = 0$ has moved up to $y(t) \approx 0.184$, while the previous equilibrium solution with no harvesting at $y(t) = 1$ has moved down to $y(t) \approx 0.816$.

From the slope field, we see that the equilibrium solution $y(t) \approx 0.184$ is unstable, whereas the equilibrium solution $y(t) \approx 0.816$ is stable. If the population ever falls below about 0.184, or 184 fish, it will then drop to 0. This is a new feature, introduced by harvesting.

In the cell below, use `create_slope_field` to experiment by plotting the solutions to the initial value problems $y(0)=0.183$ and $y(0)=0.25$. Extend the $x$ range of your slope field until the top line is close to equlibrium. Note that if you extend it too far you will break `odeint` (why?). You may also like to increase the increments in `arange` to make the plot clearer.


```python

```

## Exercises

In the following questions we will experiment with how different levels of harvesting affect outcomes for the fish population.

This week the questions will be a combination of plots and written answers. Written answers go in blank cells, but if a cell says "Type _Markdown_ and LaTeX: $\alpha^2$", you should double-click on it to edit it.

1. Assume that the harvest is now 600 fish per year. **On the same figure,**  
  a. plot the slope field,  
  b. plot the equilibrium solutions that we found in above, and  
  c. plot the solution curves for $y(0)=1$, $y(0)=0.3$, and $y(0)=0.15$.


```python

```

1.  d. In the cell below, describe the behaviour of the fish population for each of these five initial numbers of fish.



2. a. i. Assume that $H=0.8$. Plot the slope field and five solutions, one for each equilibrium solution and one for each region between, above, or below them. You can use the equation from the lab to calculate the equilibrium solutions.


```python

```

2. a. ii. In the cell below, describe the limiting behaviour of each line.



2. b. i. Assume that $H=1$. Plot the slope field and three solutions for the equilibrium solution and the regions above and below it.


```python

```

2. b. ii. Describe the limiting behaviour of each line.



2. c. i. Assume that  𝐻=1.2. Plot the slope field and two or three solutions.


```python

```

2. c. iii. Describe the limiting behaviour of the lines.



3. Use `sns.lineplot` to plot a bifurcation graph showing the reulationship between $H$ and the locations of the critical points $c$. $H$ should go on the horizontal axis but you may find it easier to calculate $H$ in terms of $c$ by rewriting the above equations giving the equilibrium solutions. 

  Hint: your plots won't work unless you add `sort=False` as a keyword argument to `sns.lineplot`, and you may observe some artifacts in your plot unless you also set `estimator=None`.


```python

```

4. Summarize what happens to the equilibrium solutions and their stability as $H$
is increased from 0 to beyond 1. Refer to your plots to support your answers.



5. What is a reasonable strategy for sustainable fishing in this case?
Don’t forget to allow qualitatively for minor catastrophes, such as disease or temporary overfishing.



56a75973fe37659f4f7bfeb52fbc342c8f76c929
ipynb
Jupyter Notebook
notebooks/lab-03.ipynb
BenKaehler/mm-labs
5409ba7f6a4d4edb802c96e4bfc47e477fc84329
BSD-3-Clause
notebooks/lab-03.ipynb
BenKaehler/mm-labs
5409ba7f6a4d4edb802c96e4bfc47e477fc84329
BSD-3-Clause
notebooks/lab-03.ipynb
BenKaehler/mm-labs
5409ba7f6a4d4edb802c96e4bfc47e477fc84329
BSD-3-Clause
2020-07-27T06:33:02.000Z
2020-07-27T07:30:41.000Z
Qwen/Qwen-72B
1. YES
2. YES