Figure 5.12:

Fractional error in the QSSA concentration of C versus dimensionless time for the series-parallel reaction, A<-> B-> C.

Code for Figure 5.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
# Converted from schm2_error.m
import numpy as np
from scipy.integrate import solve_ivp
from misc import save_ascii

k1 = 1.0
t = np.arange(0, 3.01, 0.01)
t_tmp = t[1:]
Initial = np.array([1., 0., 0.])
sqeps = np.sqrt(np.finfo(float).eps)

def rxrate(t, x, k_1, k2):
    r1 = k1*x[0] - k_1*x[1]
    r2 = k2*x[1]
    return np.array([-r1, r1 - r2, r2])

def run_case(k_1, k2):
    sol = solve_ivp(lambda t, x: rxrate(t, x, k_1, k2), [0, 3], Initial,
                    method='Radau', t_eval=t, rtol=sqeps, atol=sqeps)
    alpha = k1 * k2 / (k_1 + k2)
    c_Css = 1 - np.exp(-alpha * t_tmp)
    c_Cexact = sol.y[2, 1:]
    EE = np.abs((c_Cexact - c_Css) / c_Cexact)
    return EE

EE1 = run_case(100, 100)
EE2 = run_case(1, 1000)
EE3 = run_case(1000, 1)
EE4 = run_case(1000, 1000)

answer3 = np.column_stack([t_tmp, EE1, EE2, EE3, EE4])
save_ascii('schm2_error.dat', answer3)