

In [36]:

    
%matplotlib inline

Fish model

Setting up the model backbone

Before we show/code equations... we should start by setting up a the framework... Kinda' a "dummy" empty model...

Copy-paste into Spyder's EDITOR



In [37]:

    
days = 365 * 2 # Two years
dt   = 0.01 # units: days

Copy-paste into Spyder's CONSOLE



In [38]:

    
days









    Out[38]:





730



In [39]:

    
dt









    Out[39]:





0.01



In [40]:

    
type(days)









    Out[40]:





int



In [41]:

    
type(dt)









    Out[41]:





float

Make time vector

Copy-paste into Spyder's EDITOR



In [42]:

    
days = 365 * 2 # Two years
dt   = 0.01 # units: days 

# Setup the framework    # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<
NoSTEPS = int(days / dt) # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<

Copy-paste into Spyder's CONSOLE



In [43]:

    
NoSTEPS









    Out[43]:





73000



In [44]:

    
type(NoSTEPS)









    Out[44]:





int

Copy-paste into Spyder's EDITOR



In [45]:

    
import numpy as np # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<

days = 365 * 2 # Two years
dt   = 0.01 # units: days 

# Setup the framework
NoSTEPS = int(days / dt)
time = np.linspace(0,days,NoSTEPS) # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<

Copy-paste into Spyder's CONSOLE



In [46]:

    
time









    Out[46]:





array([  0.00000000e+00,   1.00001370e-02,   2.00002740e-02, ...,
         7.29980000e+02,   7.29990000e+02,   7.30000000e+02])



In [47]:

    
time[0]









    Out[47]:





0.0



In [48]:

    
time[-1]









    Out[48]:





730.0



In [49]:

    
time[:10]









    Out[49]:





array([ 0.        ,  0.01000014,  0.02000027,  0.03000041,  0.04000055,
        0.05000068,  0.06000082,  0.07000096,  0.0800011 ,  0.09000123])



In [50]:

    
time[-10:]









    Out[50]:





array([ 729.90999877,  729.9199989 ,  729.92999904,  729.93999918,
        729.94999932,  729.95999945,  729.96999959,  729.97999973,
        729.98999986,  730.        ])



In [51]:

    
type(time)









    Out[51]:





numpy.ndarray



In [52]:

    
len(time)









    Out[52]:





73000



In [53]:

    
time.shape









    Out[53]:





(73000L,)

Lets make one "dummy" variable. For now, lets call it Biomass (B). Lets start with and empty vector.

Copy-paste into Spyder's EDITOR



In [54]:

    
import numpy as np

days = 365 * 2 # Two years
dt   = 0.01 # units: days 

# Setup the framework
NoSTEPS = int(days / dt)
time = np.linspace(0,days,NoSTEPS)

B = np.zeros((NoSTEPS,),float)  # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<

Copy-paste into Spyder's CONSOLE



In [55]:

    
B









    Out[55]:





array([ 0.,  0.,  0., ...,  0.,  0.,  0.])



In [56]:

    
type(B)









    Out[56]:





numpy.ndarray



In [57]:

    
B[0]









    Out[57]:





0.0



In [58]:

    
B[-1]









    Out[58]:





0.0



In [59]:

    
B[:10]









    Out[59]:





array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.])

Make our first plot

Copy-paste into Spyder's EDITOR



In [60]:

    
import numpy as np
import matplotlib.pyplot as plt # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<

days = 365 * 2 # Two years
dt   = 0.01 # units: days 

# Setup the framework
NoSTEPS = int(days / dt)
time = np.linspace(0,days,NoSTEPS)

B = np.zeros((NoSTEPS,),float)

# Make plot
fig, (ax) = plt.subplots(1,1) # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<

Make initial condition for B



In [61]:

    
import numpy as np
import matplotlib.pyplot as plt

days = 365 * 2 # Two years
dt   = 0.01 # units: days 

# Setup the framework
NoSTEPS = int(days / dt)
time = np.linspace(0,days,NoSTEPS)

B = np.zeros((NoSTEPS,),float)

# Make plot
fig, (ax) = plt.subplots(1,1)
ax.plot(time/365,B,'b-')      # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<









    Out[61]:





[<matplotlib.lines.Line2D at 0x8025748>]



In [62]:

    
import numpy as np
import matplotlib.pyplot as plt

days = 365 * 2 # Two years
dt   = 0.01 # units: days 

# Setup the framework
NoSTEPS = int(days / dt)
time = np.linspace(0,days,NoSTEPS)

B = np.zeros((NoSTEPS,),float)

# Make plot
fig, (ax) = plt.subplots(1,1)
ax.plot(time/365,B,'b-')      
ax.set_xlabel('Time (years)'); ax.set_ylabel('B') # <<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<









    Out[62]:





<matplotlib.text.Text at 0x8da6b38>



In [63]:

    
import numpy as np
import matplotlib.pyplot as plt

days = 365 * 2 # Two years
dt   = 0.01 # units: days

# Initial conditions
InitCond = {}       # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<
InitCond['B'] = 0.1 # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<

# Setup the framework
NoSTEPS = int(days / dt)
time = np.linspace(0,days,NoSTEPS)

B = np.zeros((NoSTEPS,),float)

# Make plot
fig, (ax) = plt.subplots(1,1)
ax.plot(time/365,B,'b-')
ax.set_xlabel('Time (years)'); ax.set_ylabel('B')









    Out[63]:





<matplotlib.text.Text at 0x7fbb048>

Copy-paste into Spyder's CONSOLE



In [64]:

    
InitCond









    Out[64]:





{'B': 0.1}



In [65]:

    
type(InitCond)









    Out[65]:





dict



In [66]:

    
InitCond.keys()









    Out[66]:





['B']



In [67]:

    
InitCond['B']









    Out[67]:





0.1

Initialize B

Copy-paste into Spyder's EDITOR



In [68]:

    
import numpy as np
import matplotlib.pyplot as plt

days = 365 * 2 # Two years
dt   = 0.01 # units: days

# Initial conditions
InitCond = {}
InitCond['B'] = 0.1 

# Setup the framework
NoSTEPS = int(days / dt)
time = np.linspace(0,days,NoSTEPS)

B = np.zeros((NoSTEPS,),float)

# Initilizing
B[0] = InitCond['B'] # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<

# Make plot
fig, (ax) = plt.subplots(1,1)
ax.plot(time/365,B,'b-')
ax.set_xlabel('Time (years)'); ax.set_ylabel('B')









    Out[68]:





<matplotlib.text.Text at 0x98ea5c0>

Copy-paste into Spyder's CONSOLE



In [69]:

    
B[:10]









    Out[69]:





array([ 0.1,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ])

MAIN MODEL LOOP

Time now for the "big enchilada!"

Copy-paste into Spyder's EDITOR

First lets ad a "for" loop. I'll just ad a "print" statement for demostration purposes.



In [70]:

    
import numpy as np
import matplotlib.pyplot as plt

days = 365 * 2 # Two years
dt   = 0.01 # units: days

# Initial conditions
InitCond = {}
InitCond['B'] = 0.1 

# Setup the framework
NoSTEPS = int(days / dt)
time = np.linspace(0,days,NoSTEPS)

B = np.zeros((NoSTEPS,),float)

# Initilizing
B[0] = InitCond['B']


# MAIN MODEL LOOP ************************************************************************
#for t in range(0,NoSTEPS-1):
    
    #print t
# end of main model LOOP******************************************************************



# Make plot
fig, (ax) = plt.subplots(1,1)
ax.plot(time/365,B,'b-')
ax.set_xlabel('Time (years)'); ax.set_ylabel('B')









    Out[70]:





<matplotlib.text.Text at 0x9b649e8>

Lets change the "print" for the an actual "dummy model". At its minimum... it needs a dBdt and a B[t+1] equations



In [71]:

    
import numpy as np
import matplotlib.pyplot as plt

days = 365 * 2 # Two years
dt   = 0.01 # units: days

# Initial conditions
InitCond = {}
InitCond['B'] = 0.1 

# Setup the framework
NoSTEPS = int(days / dt)
time = np.linspace(0,days,NoSTEPS)

B = np.zeros((NoSTEPS,),float)

# Initilizing
B[0] = InitCond['B']


# MAIN MODEL LOOP ************************************************************************
for t in range(0,NoSTEPS-1):
    
    dBdt = 0
    
    # Time stepping
    B[t+1] = B[t] + (dBdt * dt)
# end of main model LOOP******************************************************************



# Make plot
fig, (ax) = plt.subplots(1,1)
ax.plot(time/365,B,'b-')
ax.set_xlabel('Time (years)'); ax.set_ylabel('B')









    Out[71]:





<matplotlib.text.Text at 0x9e94ba8>

Copy-paste into Spyder's CONSOLE



In [72]:

    
B[:10]









    Out[72]:





array([ 0.1,  0.1,  0.1,  0.1,  0.1,  0.1,  0.1,  0.1,  0.1,  0.1])

FISH MODEL: Partial differential equations

(1) $$ \frac{\partial B}{\partial t} = A - R $$

(2) $$ A = \epsilon \cdot IE \cdot AE \cdot Food $$

(3) $$ Food = \gamma \cdot B $$

(4) $$ R = (B \cdot m) + ( \beta \cdot A) $$

Symbol	----------Units-------------	Value	Description
$B$	$mmol N \cdot ind^{-3}$	Initial B = 1104	Fish biomass
$A$	$mmol N \cdot ind^{-3} \cdot d^{-1}$	-	Assimilation rate
$R$	$mmol N \cdot ind^{-3} \cdot d^{-1}$	-	Respiration rate
$Food$	$mmol N \cdot ind^{-3} \cdot d^{-1}$	-	Amount of food delivered to individual fish per day
$\gamma$	$d^{-1}$	0.0085	Food delivery rate
$\epsilon$	$dimensionless$	0.92	Fraction of food swallowed by fish (part of this is spit out)
$IE$	$dimensionless$	0.95	Ingestion Efficiency: Fraction of swallowed food that makes it to the stomach
$AE$	$dimensionless$	0.885	Absorption Efficiency: Fraction of ingested food that is digestible
$m$	$d^{-1}$	0.001225	Weight-specific maintenance respiration rate
$\beta$	$dimensionless$	0.44	Cost of growth coefficient

Copy-paste into Spyder's EDITOR



In [73]:

    
import numpy as np
import matplotlib.pyplot as plt

days = 365 * 2 # Two years
dt   = 0.01 # units: days

# Initial conditions
InitCond = {}
InitCond['B'] = 0.1 

# Setup the framework
NoSTEPS = int(days / dt)
time = np.linspace(0,days,NoSTEPS)

B = np.zeros((NoSTEPS,),float)

# Initilizing
B[0] = InitCond['B']


# MAIN MODEL LOOP ************
for t in range(0,NoSTEPS-1):
    
    A = 0    # Eq. 2 TO DO <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<<<<<<<<<<<
    R = 0    # Eq. 3 TO DO <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<<<<<<<<<<<
    
    dBdt = A - R # Eq.1<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<<<<<<<<<<<
    
    # Time stepping
    B[t+1] = B[t] + (dBdt * dt)
# end of main model LOOP******



# Make plot
fig, (ax) = plt.subplots(1,1)
ax.plot(time/365,B,'b-')
ax.set_xlabel('Time (years)'); ax.set_ylabel('B')









    Out[73]:





<matplotlib.text.Text at 0xa5b7748>

Ad parameters



In [74]:

    
import numpy as np
import matplotlib.pyplot as plt

days = 365 * 2 # Two years
dt   = 0.01 # units: days

# Parameters
par = {}         # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<<<<<<<<<<<
par['gamma'] = 0.0085 # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<<<<<<<<<<<
par['epsilon'] = 0.92 # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<<<<<<<<<<<
par['IE'] = 0.95 # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<<<<<<<<<<<
par['AE'] = 0.885 # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<<<<<<<<<<<
par['m'] = 0.001225 # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<<<<<<<<<<<
par['beta'] = 0.44 # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<<<<<<<<<<<

# Initial conditions
InitCond = {}
InitCond['B'] = 0.1

# Setup the framework
NoSTEPS = int(days / dt)
time = np.linspace(0,days,NoSTEPS)

B = np.zeros((NoSTEPS,),float)

# Initilizing
B[0] = InitCond['B']


# MAIN MODEL LOOP ************
for t in range(0,NoSTEPS-1):
    
    A = 0    # Eq. 2 TO DO 
    R = 0    # Eq. 3 TO DO 
    
    dBdt = A - R # Eq.1
    
    # Time stepping
    B[t+1] = B[t] + (dBdt * dt)
# end of main model LOOP******



# Make plot
fig, (ax) = plt.subplots(1,1)
ax.plot(time/365,B,'b-')
ax.set_xlabel('Time (years)'); ax.set_ylabel('B')









    Out[74]:





<matplotlib.text.Text at 0xac4fa20>



In [75]:

    
par









    Out[75]:





{'AE': 0.885,
 'IE': 0.95,
 'beta': 0.44,
 'epsilon': 0.92,
 'gamma': 0.0085,
 'm': 0.001225}



In [76]:

    
par.keys()









    Out[76]:





['AE', 'epsilon', 'm', 'beta', 'IE', 'gamma']



In [77]:

    
type(par)









    Out[77]:





dict

Add Assimilation Eq.



In [78]:

    
import numpy as np
import matplotlib.pyplot as plt

days = 365 * 2 # Two years
dt   = 0.01 # units: days

# Initial conditions
InitCond = {}
InitCond['B'] = 0.1 

# Setup the framework
NoSTEPS = int(days / dt)
time = np.linspace(0,days,NoSTEPS)

B = np.zeros((NoSTEPS,),float)

# Initilizing
B[0] = InitCond['B']


# MAIN MODEL LOOP ************
for t in range(0,NoSTEPS-1):
    
    Food = 1  # TO DO # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<<<<<<<<<
    A = par['epsilon'] * par['IE'] * par['AE'] * Food  # <<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<
    R = 0    # TO DO 
    
    dBdt = A - R 
    
    # Time stepping
    B[t+1] = B[t] + (dBdt * dt)
# end of main model LOOP******



# Make plot
fig, (ax) = plt.subplots(1,1)
ax.plot(time/365,B,'b-')
ax.set_xlabel('Time (years)'); ax.set_ylabel('B')









    Out[78]:





<matplotlib.text.Text at 0xae96e48>

Ad Food Eq.



In [79]:

    
import numpy as np
import matplotlib.pyplot as plt

days = 365 * 2 # Two years
dt   = 0.01 # units: days

# Initial conditions
InitCond = {}
InitCond['B'] = 0.1 

# Setup the framework
NoSTEPS = int(days / dt)
time = np.linspace(0,days,NoSTEPS)

B = np.zeros((NoSTEPS,),float)

# Initilizing
B[0] = InitCond['B']


# MAIN MODEL LOOP ************
for t in range(0,NoSTEPS-1):
    
    Food = par['gamma'] * B[t]  #  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<<<
    A = par['epsilon'] * par['IE'] * par['AE'] * Food  
    R = 0    # TO DO 
    
    dBdt = A - R 
    
    # Time stepping
    B[t+1] = B[t] + (dBdt * dt)
# end of main model LOOP******



# Make plot
fig, (ax) = plt.subplots(1,1)
ax.plot(time/365,B,'b-')
ax.set_xlabel('Time (years)'); ax.set_ylabel('B')









    Out[79]:





<matplotlib.text.Text at 0xb22f780>



In [80]:

    
import numpy as np
import matplotlib.pyplot as plt

days = 365 * 2 # Two years
dt   = 0.01 # units: days

# Initial conditions
InitCond = {}
InitCond['B'] = 0.1 

# Setup the framework
NoSTEPS = int(days / dt)
time = np.linspace(0,days,NoSTEPS)

B = np.zeros((NoSTEPS,),float)

# Initilizing
B[0] = InitCond['B']


# MAIN MODEL LOOP ************
for t in range(0,NoSTEPS-1):
    
    Food = par['gamma'] * B[t]
    A = par['epsilon'] * par['IE'] * par['AE'] * Food 
    
    R = (par['m']* B[t]) + (par['beta']*A)   #  <<<<<<<<<<<<<<<<<<< NEW LINE <<<<<<<<<<<<<<<<<< 
    
    dBdt = A - R 
    
    # Time stepping
    B[t+1] = B[t] + (dBdt * dt)
# end of main model LOOP******



# Make plot
fig, (ax) = plt.subplots(1,1)
ax.plot(time/365,B,'b-')
ax.set_xlabel('Time (years)'); ax.set_ylabel('B')









    Out[80]:





<matplotlib.text.Text at 0xb388710>



In [ ]: