Week 5 and Discrete Modeling¶
The first thing we should look at is what is the difference between continuous and discrete.
We know that continuous models usually exists in the "reals" $ \mathbb{R} $ and that there is an uncountable number of reals, that means they are not enumeralable, I cannot get $ x_1, x_2, x_3 ... $ etc etc to represent the reals. This is rather disappointing for computers, since we can only get an approximation of irrational numbers, for example $ \pi $. There is a value that numpy has to represent $ \pi $, but it is only a representation.
So when we look at discrete models we are looking not at an infinite approximation of the instantaneous change, like in Calculus, but we look at the step change that can happen over a finite time. The easiest example of this is the daily updates of a simple interest account.
$ X_{n+1} = X_{n} + \text{interest} $
In [2]:
import matplotlib.pyplot as plt
import numpy as np
%config InlineBackend.figure_format='retina'
In [20]:
X_0 = 150
rate = 0.14
Compound = 1
Years = 10
compound_rate= rate/Compound
Total_compounds = Compound*Years
Amount=[X_0]
Steps = [0]
X_old = Amount[-1]
for t in range(1,Total_compounds):
X_new = X_old + X_0*compound_rate
Amount.append(X_new)
Steps.append(t)
X_old = X_new
plt.plot(Steps,Amount, '.',alpha = 0.3, c='r')
plt.ylabel('Amount in $')
plt.xlabel('Time in compounds')
plt.title(' Yawn, simple interest')
plt.show()
Lets do compound interest¶
$$ X_{n+1} = X_n +X_n*\text{rate} $$$$ X_{n+1} = X_n(1 +\text{rate}) $$In [25]:
X_0 = 150
rate = 0.14
Compound = 20
Years = 1
compound_rate= rate/Compound
Total_compounds = Compound*Years
Amount=[X_0]
Steps = [0]
X_old = Amount[-1]
for t in range(1,Total_compounds):
X_new = X_old*(1+compound_rate)
Amount.append(X_new)
Steps.append(t)
X_old = X_new
plt.plot(Steps,Amount, '.',alpha = 0.3, c='b')
plt.ylabel('Amount in $')
plt.xlabel('Time in compounds')
plt.title(' Wow, compound interest')
plt.show()
In [29]:
X_0 = 150
rate = 0.14
Compound = 20
Years = 7
compound_rate= rate/Compound
Total_compounds = Compound*Years
Amount=[X_0]
Steps = [0]
X_old = Amount[-1]
for t in range(1,Total_compounds):
X_new = X_old*(1+compound_rate)
Amount.append(X_new)
Steps.append(t)
X_old = X_new
plt.plot(Steps,Amount, '.',alpha = 0.3, c='b', label = 'compound')
Amount=[X_0]
Steps = [0]
X_old = Amount[-1]
for t in range(1,Total_compounds):
X_new = X_old + X_0*compound_rate
Amount.append(X_new)
Steps.append(t)
X_old = X_new
plt.plot(Steps,Amount, '.',alpha = 0.3, c='r', label = 'simple')
plt.ylabel('Amount in $')
plt.xlabel('Time in compounds')
plt.title(' Comparison between simple and compound')
plt.legend()
plt.show()
Lets compare¶
In [44]:
X_0 = 150
rate = 0.74
Compound = 2
Years = 10
compound_rate= rate/Compound
Total_compounds = Compound*Years
Amount=[X_0]
Steps = [0]
X_old = Amount[-1]
for t in range(1,Total_compounds):
X_new = X_old*(1+compound_rate)
Amount.append(X_new)
Steps.append(t)
X_old = X_new
plt.plot(Steps,Amount, '.',alpha = 0.3, c='b')
def Compound_Exact(c):
return X_0*(1+compound_rate)**(c)
Steps = np.array(Steps)
plt.plot(Steps,Compound_Exact(Steps),'.', c='r',alpha = 0.3)
plt.ylabel('Amount in $')
plt.xlabel('Time in compounds')
plt.title(' Wow, compound interest')
plt.show()
Discrete Model Equation¶
Ricker's Equation¶
$$ X_{n+1} = X_{n} e^{r_0(1-\frac{X_n}{K})} $$
Models Fish populations
$r_0 = $ max per fish growth rate
$K = $ carrying capacity
In [57]:
r0 =0.5
K = 10000
X_0 = 100
Steps = 30
Pops = []
Times = []
X_old = X_0
for t in range(0,Steps):
X_new = X_old*np.exp(r0*(1-X_old/K))
Pops.append(X_new)
Times.append(t)
X_old = X_new
plt.plot(Times,Pops, '.-', alpha = 0.3, c='r')
Out[57]:
In [70]:
plt.figure(figsize=(13,5))
r = 0.1
while r <0.9:
r0 =r
K = 10000
X_0 = 100
Steps = 30
Pops = []
Times = []
X_old = X_0
for t in range(0,Steps):
X_new = X_old*np.exp(r0*(1-X_old/K))
Pops.append(X_new)
Times.append(t)
X_old = X_new
plt.plot(Times,Pops, '.-', alpha = 0.3, label = round(r,1))
r+=0.1
plt.legend(loc='center left', bbox_to_anchor= (1,0.5))
plt.title('Fish populations for various fish growth rates')
plt.xlabel('Time in step sizes')
plt.ylabel('Number of fish')
plt.show()
Is is possible to create a $(X_\text{new},X_\text{Old})$ graph
In [83]:
r0 =0.6
K = 1000
X_0 = 20
Steps = 30
X_old = X_0
for t in range(0,Steps):
X_new = X_old*np.exp(r0*(1-X_old/K))
plt.plot(X_old,X_new, '.', alpha = 0.3, c='r')
plt.ylim(0,K)
X_old = X_new
X = np.arange(0,K)
Y = X*np.exp(r0*(1-X/K))
plt.plot(X,Y,c='b',alpha=0.2)
plt.plot(X,X, c='r', alpha = 0.4)
Out[83]:
Lets create a stairway to represent the movement of the points, from X_old to X_new
In [90]:
plt.figure(figsize=(13,5))
r0 =0.6
K = 1000
X_0 = 200
Steps = 30
X_old = X_0
for t in range(0,Steps):
X_new = X_old*np.exp(r0*(1-X_old/K))
plt.plot(X_old,X_new, '.', alpha = 0.3, c='r')
plt.vlines(X_old, X_old,X_new,colors='g',alpha = 0.3)
plt.hlines(X_new,X_old,X_new, colors='g', alpha=0.3)
X_old = X_new
X = np.arange(0,K)
Y = X*np.exp(r0*(1-X/K))
plt.plot(X,Y,c='b',alpha=0.2)
plt.plot(X,X, c='r', alpha = 0.4)
plt.ylim(0,K)
plt.show()
In [95]:
plt.figure(figsize=(13,5))
r0 =0.6
K = 1000
X_0 = 200
Steps = 30
Pops = []
Times = []
X_old = X_0
plt.subplot(2,1,1)
for t in range(0,Steps):
Pops.append(X_old)
Times.append(t)
X_new = X_old*np.exp(r0*(1-X_old/K))
plt.plot(X_old,X_new, '.', alpha = 0.3, c='r')
plt.vlines(X_old, X_old,X_new,colors='g',alpha = 0.3)
plt.hlines(X_new,X_old,X_new, colors='g', alpha=0.3)
X_old = X_new
X = np.arange(0,K)
Y = X*np.exp(r0*(1-X/K))
plt.plot(X,Y,c='b',alpha=0.2)
plt.plot(X,X, c='r', alpha = 0.4)
plt.ylim(0,K)
plt.xlim(0,K)
plt.subplot(2,1,2)
plt.plot(Pops,Times,'.--')
plt.xlim(0,K)
plt.show()
Can we formailze the Cobwebbing Diagram¶
In [1]:
def Cobweb(old,new,color,alpha):
plt.vlines(old, old,new,colors=color,alpha =alpha)
plt.hlines(new,old,new, colors=color, alpha=alpha)
In [10]:
plt.figure(figsize=(13,5))
r0 =.5
K = 1000
X_0 = 200
Steps = 30
Max_x=K
X_old = X_0
for t in range(0,Steps):
X_new = X_old*np.exp(r0*(1-X_old/K))
plt.plot(X_old,X_new, '.', alpha = 0.3, c='r')
Cobweb(X_old,X_new,'g',0.2)
X_old = X_new
Max_x=max(Max_x,X_new)
X = np.arange(0,Max_x)
Y = X*np.exp(r0*(1-X/K))
plt.plot(X,Y,c='b',alpha=0.2)
plt.plot(X,X, c='r', alpha = 0.4)
plt.ylim(0,Max_x*1.1)
plt.show()
In [21]:
plt.figure(figsize=(13,5))
r0 =2
K = 1000
X_0 = 50
Steps = 30
Pops = []
Times = []
X_old = X_0
plt.subplot(2,1,1)
for t in range(0,Steps):
Pops.append(X_old)
Times.append(t)
X_new = X_old*np.exp(r0*(1-X_old/K))
plt.plot(X_old,X_new, '.', alpha = 0.3, c='r')
Cobweb(X_old,X_new,'g',0.2)
X_old = X_new
Max_x=max(Max_x,X_new)
X = np.arange(0,Max_x)
Y = X*np.exp(r0*(1-X/K))
plt.plot(X,Y,c='b',alpha=0.2)
plt.plot(X,X, c='r', alpha = 0.4)
plt.ylim(0,Max_x*1.1)
plt.xlim(0,Max_x*1.1)
plt.subplot(2,1,2)
plt.plot(Times,Pops,'.--')
plt.xlim(0,Times[-1])
plt.show()
In [28]:
for r in np.arange(1.5,2.3,0.1):
r0=round(r,1)
plt.figure(figsize=(13,5))
K = 1000
X_0 = 50
Steps = 30
Pops = []
Times = []
X_old = X_0
plt.subplot(2,1,1)
for t in range(0,Steps):
Pops.append(X_old)
Times.append(t)
X_new = X_old*np.exp(r0*(1-X_old/K))
plt.plot(X_old,X_new, '.', alpha = 0.3, c='r')
Cobweb(X_old,X_new,'g',0.2)
X_old = X_new
Max_x=max(Max_x,X_new)
X = np.arange(0,Max_x)
Y = X*np.exp(r0*(1-X/K))
plt.plot(X,Y,c='b',alpha=0.2)
plt.plot(X,X, c='r', alpha = 0.4)
plt.ylim(0,Max_x*1.1)
plt.xlim(0,Max_x*1.1)
plt.title("r0 = "+str(r0))
plt.subplot(2,1,2)
plt.plot(Times,Pops,'.--')
plt.xlim(0,Times[-1])
plt.show()
In [35]:
plt.figure(figsize=(13,5))
r0 =4.3
K = 1000
X_0 = 50
Steps = 300
Max_x=K
X_old = X_0
for t in range(0,Steps):
X_new = X_old*np.exp(r0*(1-X_old/K))
plt.plot(X_old,X_new, '.', alpha = 0.3, c='r')
Cobweb(X_old,X_new,'g',0.2)
X_old = X_new
Max_x=max(Max_x,X_new)
X = np.arange(0,Max_x)
Y = X*np.exp(r0*(1-X/K))
plt.plot(X,Y,c='b',alpha=0.2)
plt.plot(X,X, c='r', alpha = 0.4)
plt.ylim(0,Max_x*1.1)
plt.show()
In [34]:
r0 =3.3
K = 1000
X_0 = 50
Steps = 300
Pops = []
Times = []
X_old = X_0
for t in range(0,Steps):
X_new = X_old*np.exp(r0*(1-X_old/K))
Pops.append(X_new)
Times.append(t)
X_old = X_new
plt.plot(Times,Pops, '.-', alpha = 0.3, c='r')
Out[34]:
Logistic Equation¶
Bunny Growth on Magical Island $$ X_{n+1} = rX_n$$
On a non-magical island $$X_{n+1} = rX_n(1-X_n)$$
In [46]:
plt.figure(figsize=(13,5))
r0 =3
X_0 = 0.25
Steps = 10
Max_x=1.1
X_old = X_0
for t in range(0,Steps):
X_new = r0*X_old*(1-X_old)
plt.plot(X_old,X_new, '.', alpha = 0.3, c='r')
Cobweb(X_old,X_new,'g',0.2)
X_old = X_new
Max_x=max(Max_x,X_new)
X = np.linspace(0,Max_x,100)
Y = r0*X*(1-X)
plt.plot(X,Y,c='b',alpha=0.2)
plt.plot(X,X, c='r', alpha = 0.4)
plt.ylim(0,Max_x*1.1)
plt.show()
