Figure 5.11:

Fractional error in the QSSA concentration of C for the series reaction A -> B -> C.

Code for Figure 5.11

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
56
# Converted from schm1_error.m
import numpy as np
from scipy.integrate import solve_ivp
from misc import save_ascii

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

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

sol1 = solve_ivp(lambda t, x: rxrate(t, x, 1.0), [0, 1], Initial,
                 method='Radau', t_eval=t, rtol=sqeps, atol=sqeps)
answer = np.column_stack([t, sol1.y.T])

sol2 = solve_ivp(lambda t, x: rxrate(t, x, 10.0), [0, 1], Initial,
                 method='Radau', t_eval=t, rtol=sqeps, atol=sqeps)
answer = np.column_stack([answer, sol2.y.T])

sol3 = solve_ivp(lambda t, x: rxrate(t, x, 50.0), [0, 1], Initial,
                 method='Radau', t_eval=t, rtol=sqeps, atol=sqeps)
answer = np.column_stack([answer, sol3.y.T])

sol4 = solve_ivp(lambda t, x: rxrate(t, x, 1000.0), [0, 1], Initial,
                 method='Radau', t_eval=t, rtol=sqeps, atol=sqeps)
answer = np.column_stack([answer, sol4.y.T])

sol5 = solve_ivp(lambda t, x: rxrate(t, x, 10000.0), [0, 1], Initial,
                 method='Radau', t_eval=t, rtol=sqeps, atol=sqeps)
answer = np.column_stack([answer, sol5.y.T])

sol6 = solve_ivp(lambda t, x: rxrate(t, x, 100000.0), [0, 1], Initial,
                 method='Radau', t_eval=t, rtol=sqeps, atol=sqeps)
answer = np.column_stack([answer, sol6.y.T])

# Strip first row
answer = answer[1:, :]
t_tmp = t[1:]
c_Css = 1 - np.exp(-k1 * t_tmp)

# Columns 13, 16, 19 in 1-indexed = columns 12, 15, 18 in 0-indexed (C col of 4th, 5th, 6th solutions)
c_Cexact = answer[:, 12]
EE1 = np.abs((c_Cexact - c_Css) / c_Cexact)

c_Cexact = answer[:, 15]
EE2 = np.abs((c_Cexact - c_Css) / c_Cexact)

c_Cexact = answer[:, 18]
EE3 = np.abs((c_Cexact - c_Css) / c_Cexact)

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