In [ ]:
fig, ax = plt.subplots()
ax.plot(sol.t, sol.y[0])
ax.set_xlabel('t')
ax.set_ylabel('S(t)')
Out[ ]:
Text(0, 0.5, 'S(t)')
Find the solution to the difference equation: $x_{t+1} = x_t - 4$, where $x_0 = 0$
(EMEA exercise 11.8.2a)
Consider the difference equation $x_{t+1} = \frac{1}{2}x_t + 2$, where $x_0 = 6$. Find a solution of the differential equation. Does the dynamic system have a steady-state? If so, what is it?
(Based on EMEA 11.8.2b)
Suppose we have a differential equation $\dot{x} = .1x$, where $x_0 = 5$. Find a solution.
The following differential equation defines the biological growth of a fish stock in a lake
$$\dot{x} = 0.05(200-x(t))$$Find a solution to the differential equation.
(emea example 11.9.2)
Our focus has been looking at discrete difference equations and in the next lab we will continue to focus on computationally simulating difference equations and looking at some applications. But we also have some tools available in Python to numerically approximate solutions to continuous time differential equations. PPNM works through the math of some of these approximation methods, but we will jump direct to the tools implimented in python.
This will be based on PPNM Ch 22.6
We begin with a simple differential equation
$$\frac{dS(t)}{dt} = cos(t)$$Where we have an initial value of 0:
$$S(0) = 0$$To numerically solve this problem we use the solve_ivp function from scipy. ivp stands for initial value problem.
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
First we define the range of t to be between 0 and $\pi$ with a grid that has intervals h=.1
t_eval = np.arange(0, np.pi, 0.1)
t_eval
array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. , 1.1, 1.2,
1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2. , 2.1, 2.2, 2.3, 2.4, 2.5,
2.6, 2.7, 2.8, 2.9, 3. , 3.1])
Then we define our our differential equation as a function:
def F(t,s):
return(np.cos(t))
#Or equivalently using *lambda* notation
#F = lambda t, s: np.cos(t)
And then we use the solver to find the numerical solution
sol = solve_ivp(F, [0, np.pi], [0], t_eval=t_eval)
sol
message: The solver successfully reached the end of the integration interval.
success: True
status: 0
t: [ 0.000e+00 1.000e-01 ... 3.000e+00 3.100e+00]
y: [[ 0.000e+00 9.983e-02 ... 1.416e-01 4.177e-02]]
sol: None
t_events: None
y_events: None
nfev: 38
njev: 0
nlu: 0
And that is it. So the hard part here is really just setting up a dynamic model, then letting the computer solve for it. Now we can plot the solution
fig, ax = plt.subplots()
ax.plot(sol.t, sol.y[0])
ax.set_xlabel('t')
ax.set_ylabel('S(t)')
Text(0, 0.5, 'S(t)')
Here we get the solution to the problem. Recall from our lecture on numerical integration that the analytic solution to this problem is simply the sin function. We can then chart the "error" of the numeric approach
fig, ax = plt.subplots()
ax.plot(sol.t, sol.y[0] - np.sin(sol.t))
ax.set_xlabel('t')
ax.set_ylabel('S(t) - sin(t)')
Text(0, 0.5, 'S(t) - sin(t)')
If you set the step-size to h=.01, how does this change the magnitude of the errors?
Consider the differential equation $\frac{dS(t)}{t} = -S(t)$. Solve and plot the solution. What is the exact analytic solution? Show the error.
Now let's consider a differential equation with multiple variables
$$\frac{d x(t)}{dt} = -ty(t)$$$$\frac{d y(t)}{dt} = t^2x(t)$$Where the initial values are $x(0) = 1$, and $y(0)=1$
We can write this in matrix form defining:
$$S(t) = \begin{bmatrix} x(t) \\ y(t) \end{bmatrix}$$And then the system of differential equations as:
$$\frac{dS(t)}{dt} = \begin{bmatrix} 0 & -t \\ t^2 & 0 \end{bmatrix} S(t)$$We can imagine that these equations describe a path on a (x,y) coordinate chart starting at (1,1). When time changes a little bit, $\Delta t$, then x (horizontal direction) changes by $-\Delta t y(t)$ and y (vertical direction) changes by $\Delta t^2x(t)$.
Notice that the direction and magnitude of change in both the x and y direction is dependent on the state or position of the other variable.
Here we define the problem, with a range of t=0 to 10
F = lambda t, s: np.dot(np.array([[0, -t], [t**2, 0]]), s)
t_eval = np.arange(0, 10.01, 0.01)
sol = solve_ivp(F, [0, 10], [1, 1], t_eval=t_eval)
sol
message: The solver successfully reached the end of the integration interval.
success: True
status: 0
t: [ 0.000e+00 1.000e-02 ... 9.990e+00 1.000e+01]
y: [[ 1.000e+00 9.999e-01 ... 1.214e-01 -9.236e-02]
[ 1.000e+00 1.000e+00 ... 2.114e+00 2.128e+00]]
sol: None
t_events: None
y_events: None
nfev: 854
njev: 0
nlu: 0
fig, ax = plt.subplots()
ax.plot(sol.y.T[:, 0], sol.y.T[:, 1])
ax.set_xlabel('x')
ax.set_ylabel('y')
Text(0, 0.5, 'y')
Now we will explore a somewhat more involved model. Here we wish to model "production" or "harvesting" of a natural resource where we need to consider the biological dynamics in our model. To be concrete, let's think of this as a wild fishery, perhaps for cod outside the coast of Norway, and a company (or the government) wants to decide how much to harvest each year.
We start with a difference equation describing the growth of the biological "stock" of fish:
$$\frac{dB(t)}{dt} = rB(t)(1-\frac{B(t)}{k})$$Here we model the growth rate as two components. First is the $rB(t)$ where:
This is the easiest to interpret. If we assume that all else equal the fish stock will grow by 10\% per year (or whatever the frequency is), then this would be (.10*B(t)).
But a model which allows for unending growth in the stock of fish would be highly unrealistic. Natural environments and ecosystems have limits to how much fish they can support. That is what the second term, $(1-\frac{B(t)}{k})$, attempts to put into the model:
So let's say that the k=1000 (units are arbitrary, perhaps tons of fish), and the total population is only B=100. Then $1-\frac{B(t)}{k} = 1-\frac{100}{1000} = .9$. This means that the intrinsic rate of growth from the first term is not strongly hampered by the natural carrying capacity of the environment. But if the stock of fish were B=900, then $1-\frac{B(t)}{k} = 1-\frac{900}{1000} = .1$, which would strongly limit the growth rate of the fish stock.
Now let's consider the economics side of this bioeconomic model. We will devise a type of production function, here called a "harvesting" function:
$$Y(t)=qf(t)B(t)$$The definitions of the terms are:
Y(t): Catch rate at time t, Y(t)
f(t): fishing effort at time t
q: "catchability quotient"- fraction of the population fished by an effort unit
The idea here is that the amount we are able to harvest is the combination of how abundant fish are (B) and how much effort is put into fishing f(t), where q is a parameter we choose that determines the link between effort and how much of the fish population is harvested.
Then we can update our differencial equation to reflect the effect of a commercial harvest:
$$\frac{dB(t)}{dt} = rB(t)(1-\frac{B(t)}{k})-Y(t)$$$$\frac{dB(t)}{dt} = rB(t)(1-\frac{B(t)}{k})-qf(t)B(t)$$We now have a differential equation that combines a biological model of species growth with an economic "production function" or harvesting. We can start to simulate
We setup the bounds of our model, some of the values are chosen arbitrarily and some are chosen through trial-and-error.
t_eval = np.arange(0, 200, 0.1)
r=.10
f=1
q=.04
k=1000
B0 = [100]
Above we let the time go from 0 to 200 (we can think of these as weeks or months. We discretize our time into periods of h=.1
Our intrinsic rate of biological growth (assuming no environmental constraints) is 10 \%
For simplicity, we fix our effort level at a constant and normalized f=1
Our q is set at 4\%, so with a constant effort level of 1, this means we are harvesting 4\% of the biological stock every period.
The carrying capacity is set to 1000
Our initial value is $B(0)=100$.
We define our differential equation as a function:
def dB(t, B):
return(r*B*(1-B/k)-q*f*B)
Then solve using our python initial value problem routine:
fish_sol = solve_ivp(dB, [0, 200], B0, t_eval=t_eval)
And we can plot our results
fig, ax = plt.subplots()
ax.plot(fish_sol.t, fish_sol.y[0])
ax.set_xlabel('t')
ax.set_ylabel('B(t)')
Text(0, 0.5, 'B(t)')
What we see is a solution that goes towards a certain asymptotic value which represents a type of steady-state equilibrium. With the parameters we choose this ends up being at just under 600.
The intuition here is that this is the population level where the growth rate of the population exactly equals the harvesting so that we get a constant population stock
We can transform the population dynamics to see the harvesting dynamics by using the harvesting function:
Y = q*f*fish_sol.y[0]
fig, ax = plt.subplots()
ax.plot(fish_sol.t, Y)
ax.set_xlabel('t')
ax.set_ylabel('Y(t)')
Text(0, 0.5, 'Y(t)')
We could ask another question: what is the optimal amount of harvesting. The solution to this problem requires a somewhat more involved methodology called dynamic programming which is beyond the scope of this course. But we could experiment with various effort levels to see what leads to the highest equilibrium catch.
Now assume that effort f(t) is not fixed but is related to state of the fish stock (relative to carrying capacity, B/k). That is, when B/k is higher then we require less effort and when B/k is lower effort increases. Thus we can write our differential equation for B(t) the same:
$$\frac{dB(t)}{dt} = rB(t)(1-\frac{B(t)}{k})-qf(t)B(t)$$But now we also write a differential equation for f(t):
$$\frac{df(t)}{dt} = \gamma (\omega - \frac{ B(t)}{k})$$Which simply states that the change in effort is related to the relative size of the fish stock, where $\gamma$ is a parameter that tells us how sensitive this effect is and $\omega$ is the point where the sign changes direction (is effort increasing or decreasing.
We can then write this system of differential equations in matrix form:
$$S(t) = \begin{bmatrix} B(t) \\ f(t) \end{bmatrix}$$And then the system of differential equations as:
$$\frac{dS(t)}{dt} = \begin{bmatrix} dB(t)/ dt \\ df(t)/dt \end{bmatrix} = \begin{bmatrix} r(1-\frac{B(t)}{k}) & -qB(t) \\ \gamma (\frac{\omega}{B(t)}- \frac{1}{k}) & 0 \end{bmatrix} S(t)$$Confirm by matrix multiplication that the above matrix system of equations corresponds to the two equations we wrote above
We now impliment the system of differential equations in python
def dS(t,S):
return(np.dot(np.array([[r*(1-S[0]/k), -q*S[0]],[gamma*(omega/S[0]-1/k),0]]), S))
t_eval = np.arange(0, 200, .1)
r=.10
q=.04
k=1000
S0 = [100,1] #initial values of B and f
gamma = .005
omega = .4
Above we have set $\gamma =.5$ and $\omega=.3$, $r=.10$, $q=.04$ and $k=1000$ so that the differential equation for the stock, B, is:
$$ \frac{dB(t)}{dt} = .10 B(t)(1-\frac{B(t)}{1000})-.04f(t)B(t)$$and that for effort is
$$\frac{df(t)}{t} = .02*(.4 - \frac{ B(t)}{1000})$$fishSol2 = solve_ivp(dS, [0, 200], S0, t_eval=t_eval)
fishSol2.y[0]
array([100. , 100.50097037, 101.00388298, ..., 638.67468962,
638.72207712, 638.76950323])
fig, ax = plt.subplots(2)
ax[0].plot(fishSol2.t, fishSol2.y[0])
ax[0].set_xlabel('t')
ax[0].set_ylabel('B(t)')
ax[1].plot(fishSol2.t, fishSol2.y[1])
ax[1].set_xlabel('t')
ax[1].set_ylabel('f(t)')
Text(0, 0.5, 'f(t)')
For the range and parameters I specified above, this model has some intuition. At the starting point of B(0)=100 the effort is increasing, but the stock of fish is also increasing. At a certain point (B=400), the stock of fish is getting big enough where effort starts decreasing. Thus in this model the equilibrium stock doesn't go towards an equilibrium value but continues to increase as the fishing companies continue to use less effort. Eventually the stock approaches the ecological limit of 1000 and the fishing companies are using gradually less effort. This seems like a very benign path.
But changing some of the parameters (increasing $\gamma$ for example) can lead to some bad equilibriiums where a declining fish stock leads fishing companies to harvest ever more intensely, eventually leading to depletion of the entire stock (try it!).
Still, this is probably not a great model (I made it up - a true fisheries economist would probably be horrified). Try playing with some more of the parameters - you can get effort to go negative, which probably doesn't make sense.
But this is how modelling moves forward. You try something, test it out and then figure out weaknesses, improve the model and so on. And no model is perfect - the question is whether you can extract some insight about how the real world will work.
Solve the following differential equations
a.) $\dot{x} = -3x$
b.) $\dot{x} + 4x = 12$
c.) $\dot{x} - 3x = 12x^2$
d.) $5\dot{x} = -x$
e.) $3\dot{x} + 6x = 10$
(EMEA Review question 11.16)
Suppose that a firm's capital stock K(t) satisfies the differential equation $\dot{K}(t) = I-\sigma K(t)$, where investment I is constant and $\sigma K(t)$ denotes depreciation, with $\sigma$ a positive constant.
a.) Find the solution of the equation if the capital stock at time t=0 is $K_0$
b.) Let $\sigma = 0.05$ and $I=10$. Explain what happens to K(t) as $t \rightarrow \infty$, when i.) $K_0 = 150$, ii.) $K_0 = 250$
(EMEA Exercise 11.9.2)
Now let's introduce something in the model which as 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.