# 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
from plotly.figure_factory import create_quiver
from plotly import graph_objs as go
from numpy import meshgrid, arange, sqrt, linspace
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)
    scale = 0.9*min(x[0][1]-x[0][0], y[1][0]-y[0][0])  # assume a regular grid
    fig = create_quiver(x, y, 1/L, S/L, scale=scale, arrow_scale=1e-16)
    fig.layout.update(yaxis=dict(scaleanchor='x',
                                 scaleratio=1,
                                 range=[y.min()-scale, y.max()+scale]),
                      xaxis=dict(range=[x.min()-scale, x.max()+scale]),
                      showlegend=False, width=500,
                      height=0.8*(y.max()-y.min())/(x.max()-x.min())*500)
    x = linspace(x.min(), x.max())
    for y0 in initial_values:
        y = odeint(diff_eq, y0, x, args).flatten()
        fig.add_trace(go.Scatter(x=x, y=y))
    return fig
```

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.1), arange(-0.4, 1.41, 0.1))
fig = create_slope_field(diff_eq, x, y, initial_values=(0.4, 1))
fig.show('png')
```

### 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.1), arange(-0.4, 1.41, 0.1))
fig = create_slope_field(diff_eq, x, y, args=(0.6,))
fig.show('png')
```

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
x, y = meshgrid(arange(0, 2.1, 0.1), arange(-0.4, 1.41, 0.1))
fig = create_slope_field(diff_eq, x, y, args=(0.6,), initial_values=(0.183, 0.25))
fig.show('png')
```

## Exercises

In this lab you will experiment with the population dynamics given by the logistic equation with harvesting that we started analysing in the lab.

This week the questions will be a combination of plots and written answers.

1. Assume that the harvest is 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
x, y = meshgrid(arange(0, 1.2, 0.1), arange(-0.4, 1.3, 0.1))
fig = create_slope_field(diff_eq, x, y, args=(0.6,), initial_values=(1, 0.3, 0.15))
fig.show('png')
```

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

If the initial value is y(0) = 1, the popuplation falls towards the stable equilirium at approx 0.8.
If the initial value is y(0) = 0.3, then the population tends out of the unstable equilibrium at approximately 0.2 and towards a stable equilibrium f If our initial value is at y(0) = 0.15, at after 0.6 years the population will be about 0. It zooms off to negative infinity as time goes on.

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
x, y = meshgrid(arange(0, 1.1, 0.05), arange(-0.4, 1.3, 0.05))
fig = create_slope_field(diff_eq, x, y, args=(0.8,), initial_values=(1, 0.45, 0.2, 0.7236, 0.2764))
fig.show('png')
```

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

We can see limiting behaviours for lines with initial values (1, 0.45 and 0.7236) is towards the stable equilibrium found at approx 0.7236 in infinity time. 

The limiting behaviour of the line at point 0.2764 stays at this point.

The bottom line (initial value of 0.2) tends to negative infinity.

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
x, y = meshgrid(arange(0, 1.1, 0.05), arange(-0.4, 1.41, 0.05))
fig = create_slope_field(diff_eq, x, y, args=(1,), initial_values=(1, 0.46, 0.2764))
fig.show('png')
```

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

Initial value of 1 declines and tends towards the equilibrium at approx. 0.46.

For any value below the equilibrium, it hits 0 and continues declining so it tends towards negative infinity.

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


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

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

This situation demonstrates no critical points, the bottom line (initial value 0.4) will have a limiting behaviour that tends off to negative infinity.

Y will never become 0.

The fish population will always be declining. The species will inevitably die out.

3. 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.

The higher the H value the less equilibrium points exist. Beyond 1 there exists critical points so the populaiton will inevitably decline and remain declining to negative infinity. 

4. 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.

A sustainable amount of fish to remove from the environment would be around 500. 

We have a margin of error of around 500 fish, and by taking 500 from the population we are approximately half way through this margin of error (so as to incorporate catastrophes/disease/temporary overfishing)

d1b07d3cf8c6086ed7027ee4081d38491d05ce1e
ipynb
Jupyter Notebook
notebooks/lab-03.ipynb
JoshE117/mm-labs
07558ee8dff76952a8afd6410419aaae2fffb893
BSD-3-Clause
notebooks/lab-03.ipynb
JoshE117/mm-labs
07558ee8dff76952a8afd6410419aaae2fffb893
BSD-3-Clause
notebooks/lab-03.ipynb
JoshE117/mm-labs
07558ee8dff76952a8afd6410419aaae2fffb893
BSD-3-Clause
Qwen/Qwen-72B
1. YES
2. YES