Modul 3 Kodeeksempler forskerforedrag

import matplotlib
if not hasattr(matplotlib.RcParams, "_get"):
    matplotlib.RcParams._get = dict.get

Modul 3 Kodeeksempler forskerforedrag#

Opsætning#

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import matplotlib
if not hasattr(matplotlib.RcParams, "_get"):
    matplotlib.RcParams._get = dict.get

Solve system 1#

# --------------------------------------------------------------------
# Tragedy of the commons generic ODE system - math version
# --------------------------------------------------------------------
# state vector y = (B,A)
# B: biomass
# A: actors
# --- biological parameters ---
# r: growth rate              
# K: carrying capacity      
# --- economical parameters ---
# s: sector adaptation rate 
# p: catch price             
# c: operating cost
#
# In this version parameter values are close to unity and not so
# closely reflecting real world systems - it has 5 parameters.
# By non-dimensionalization the system can further be reduced to
# 3 essential parameters. Parameter/variable units are not spelled out.
#
# It has 2 parameter settings - one with an interior equilibrium, and one with
# a boundary equilibrium (sector crash)
# --------------------------------------------------------------------
def ToC(t,y, r,K,p,c,s):
    B,A = y
    dBdt = r*B*(1-B/K) - A*B
    dAdt = s*(p*B/c - 1)
    if A<0:
        dAdt = max(dAdt, 0)     # dAdt > 0 is OK if A<0, allow fleet to rebuild
    return np.array([dBdt,dAdt])

# ---- this parameter combination corresponds to an interior equilibrium ----
r  = 1    # growth rate
K  = 1    # carrying capacity
p  = 1    # catch price  
c  = 0.6  # operating cost
s  = 1.5  # sector adaptation rate 
y0 = (K, 0) # initial condition

# ---- this parameter combination is in sector crash regime; uncomment line below
# c   = 1.4; y0  = (1.3*K, 1) 

tspan  = (0, 20)
t_eval = np.linspace(tspan[0], tspan[1], 1000)
args   = (r,K,p,c,s)
sol1   = solve_ivp(ToC, tspan, y0, t_eval=t_eval, args=args)

# -------- plotting --------

fig, (axB,axA,axphase) = plt.subplots(1,3)

axB.plot(sol1.t, sol1.y[0])
axB.set_title('B ')
axB.set(xlabel="t")

axA.plot(sol1.t, sol1.y[1])
axA.set_title('A ')
axA.set(xlabel="t")

axphase.plot(sol1.y[1], sol1.y[0])
axphase.set_title('phaseplot')
axphase.set(xlabel="A", ylabel="B")
axphase.yaxis.set_label_position("right")
axphase.yaxis.tick_right()

plt.show()

Solve system 2#

# --------------------------------------------------------------------
# Tragedy of the commons ODE system - course 25328 version
# --------------------------------------------------------------------
# state vector y = (B,A)
# B: fish / km2
# A: boats / km2
# --- biological parameters ---
# r: growth rate            [1/year]     
# K: carrying capacity      [fish/km2]
# b: habitat clearance rate [km2/boat/year] 
# --- economical parameters ---
# s: sector adaptation rate [boats/km2/year]
# p: catch price            [$/fish]       
# c: operating cost         [$/year/boat]
#
# In this version parameter values reflect a semireal situation,
# fishermen exploiting a fish population, and parameters/variables
# are with units and strictly interpretable/observable, making them
# easier to guesstimate. This is the version we use in course 25328 exercises.
# It has 2 parameter settings - one with an
# interior equilibrium, and one with a boundary equilibrium (sector crash)
# By non-dimensionalization the system can be reduced to 3 essential
# parameters. 
# --------------------------------------------------------------------
def ToC(t,y, r,K,p,b,c,s):
    B,A = y
    dBdt = r*B*(1-B/K) - b*A*B
    dAdt = s*(p*b*B/c - 1)
    if A<0:
        dAdt = max(dAdt, 0)     # dAdt > 0 is OK if A<0, allow fleet to rebuild
    return np.array([dBdt,dAdt])

# ---- this parameter combination corresponds to an interior equilibrium ----
r  = 1            # [1/year]          growth rate
K  = 1e3          # [fish/km2]        fish carrying capacity
p  = 1            # [$/fish]          catch price
b  = 50*5*1e-2*86 # [km2/boat/year]   habitat clearance rate (sail_days * fish_hours_pr_day * gear_width * sail_speed)   
c  = 1e5          # [$/year/boat]     operating cost
s  = 0.2/3/(3-1)  # [boats/km2/year]  sector adaptation rate (boat_dens / adapt_time / incentive_factor)
y0 = (K, 0) # initial condition

# ---- this parameter combination is in sector crash regime; uncomment line below
# c   = 3e5; y0  = (1.3*K, 0.001) 

tspan  = (0, 30)
t_eval = np.linspace(tspan[0], tspan[1], 1000)
args   = (r,K,p,b,c,s)
sol1   = solve_ivp(ToC, tspan, y0, t_eval=t_eval, args=args)

# -------- plotting --------

fig, (axB,axA,axphase) = plt.subplots(1,3)

axB.plot(sol1.t, sol1.y[0])
axB.set_title('B [fish/km2]')
axB.set(xlabel="t [years]")

axA.plot(sol1.t, sol1.y[1])
axA.set_title('A [boats/km2]')
axA.set(xlabel="t [years]")

axphase.plot(sol1.y[1], sol1.y[0])
axphase.set_title('phaseplot')
axphase.set(xlabel="A [boats/km2]", ylabel="B [fish/km2]")
axphase.yaxis.set_label_position("right")
axphase.yaxis.tick_right()

plt.show()
Contents