Figure 10.12:

The two steady-state biomass and substrate concentrations versus dilution rate; stable (solid), unstable (dashed).

Code for Figure 10.12

Text of the GNU GPL.

main.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
# Converted from chemoss.m - chemostat steady states (Sf=5)
import numpy as np
from misc import octave_save

mum = 1.0; Ks = 1.0; Sf = 5.0; y = 1.0

def jacob(X, S, mum_, D, Ks_, Sf_, y_):
    mug = mum_ * S / (Ks_ + S)
    tmp = Ks_ * mum_ / (Ks_ + S)**2
    return np.array([[mug - D,    X * tmp],
                     [-mug / y_, -D - X/y_ * tmp]])

Dopt = mum * (1 - np.sqrt(Ks / (Ks + Sf)))
Dc   = mum * Sf / (Ks + Sf)
nDs  = 25

def compute_tables(Dvec):
    n   = len(Dvec)
    xs1 = np.zeros((n, 2)); xs2 = np.zeros((n, 2))
    lam1 = np.zeros((n, 2)); lam2 = np.zeros((n, 2))
    for i, D in enumerate(Dvec):
        xs1[i] = [0.0, Sf]
        Ss = D * Ks / (-D + mum)
        xs2[i] = [y*(Sf - Ss), Ss]
        lam1[i] = np.sort(np.linalg.eigvals(jacob(*xs1[i], mum, D, Ks, Sf, y)).real)
        lam2[i] = np.sort(np.linalg.eigvals(jacob(*xs2[i], mum, D, Ks, Sf, y)).real)
    return xs1, xs2, lam1, lam2

Dvec1 = np.linspace(0.001, Dc, nDs)
Dvec2 = np.linspace(Dc, 0.99, nDs)

xs11, xs21, lam11, lam21 = compute_tables(Dvec1)
xs12, xs22, lam12, lam22 = compute_tables(Dvec2)
DX1 = Dvec1 * xs21[:, 0]
DX2 = Dvec2 * xs22[:, 0]

table1 = np.column_stack([Dvec1, xs11, xs21, DX1])
table2 = np.column_stack([Dvec2, xs12, xs22, DX2])
octave_save('chemoss.dat', ('table1', table1), ('table2', table2))