Figure 5.9:

Comparison of equilibrium assumption to full model for k_1=1, k_{-1}=0.5, k_2=k_{-2}=10.

Code for Figure 5.9

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
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
# Converted from rls_5.m
# Compares nscheme=1 (rx1 rls) vs nscheme=0 (full) for k2=k_2=10
import numpy as np
from scipy.integrate import solve_ivp
from misc import save_ascii

k1   = 1.0
k_1  = 0.5
k2   = 10.0
k_2  = 10.0
c_A0 = 1.0
c_B0 = 0.4
c_C0 = 0.0

def rhs(t, x, nscheme, k2, k_2):
    c_A, c_B, c_C = x
    r1 = k1*c_A - k_1*c_B
    r2 = k2*c_B - k_2*c_C
    if nscheme == 0:
        return np.array([-r1, r1 - r2, r2])
    elif nscheme == 1 or nscheme == 5:
        K2 = k2/k_2
        return np.array([-r1, 1/(1+K2)*r1, K2/(1+K2)*r1])
    elif nscheme == 2:
        K1 = k1/k_1
        return np.array([-1/(1+K1)*r2, -K1/(1+K1)*r2, r2])
    elif nscheme == 3:
        c_B = k1/(k_1+k2)*c_A + k_2/(k_1+k2)*c_C
        r1 = k1*c_A - k_1*c_B
        r2 = k2*c_B - k_2*c_C
        dB = (k1*k2*(k_2 - k1)/(k_1+k2)**2 * c_A +
              k_1*k_2*(k1 - k_2)/(k_1+k2)**2 * c_C)
        return np.array([-r1, dB, r2])

sqeps = np.sqrt(np.finfo(float).eps)
tfinal = 5.0
ntimes = 200
tout = np.linspace(0, tfinal, ntimes)

# nscheme=1: IC for rx1 rls
nscheme = 1
K2 = k2/k_2
x0 = np.array([c_A0, 1/(1+K2)*(c_B0+c_C0), K2/(1+K2)*(c_B0+c_C0)])
sol1 = solve_ivp(lambda t, x: rhs(t, x, nscheme, k2, k_2), [0, tfinal], x0,
                 method='Radau', t_eval=tout, rtol=sqeps, atol=sqeps)
table = np.column_stack([tout, sol1.y.T])

# nscheme=0: IC for full model
nscheme = 0
x0 = np.array([c_A0, c_B0, c_C0])
sol0 = solve_ivp(lambda t, x: rhs(t, x, nscheme, k2, k_2), [0, tfinal], x0,
                 method='Radau', t_eval=tout, rtol=sqeps, atol=sqeps)
table = np.column_stack([table, sol0.y.T])

save_ascii('rls_5.dat', table)