Figure 5.7:

Full model solution for k_1=1, k_{-1}=0.5, k_2=k_{-2}=10.

Code for Figure 5.7

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
# Converted from rls_2.m
# A <-> B <-> C,  k1=1, k_1=0.5, k2=10, k_2=10, nscheme=0 (full model)
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
nscheme = 0

def rhs(t, x, nscheme):
    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:
        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])

x0 = np.array([c_A0, c_B0, c_C0])
tfinal = 5.0
ntimes = 200
tout = np.linspace(0, tfinal, ntimes)
sqeps = np.sqrt(np.finfo(float).eps)

sol = solve_ivp(lambda t, x: rhs(t, x, nscheme), [0, tfinal], x0,
                method='Radau', t_eval=tout, rtol=sqeps, atol=sqeps)
x = sol.y.T

RA = -k1*x[:, 0] + k_1*x[:, 1]
RC =  k2*x[:, 1] - k_2*x[:, 2]
RB = -RA - RC

table = np.column_stack([tout, x, RA, RB, RC])
save_ascii('rls_2.dat', table)