Figure 5.10:

Normalized concentration of C versus dimensionless time for the series reaction A -> B -> C.

Code for Figure 5.10

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

k1 = 1.0
t = np.arange(0, 5.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, 5], 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, 5], 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, 5], Initial,
                 method='Radau', t_eval=t, rtol=sqeps, atol=sqeps)
answer = np.column_stack([answer, sol3.y.T])

c_Css = 1 - np.exp(-k1 * t)
answer = np.column_stack([answer, c_Css])

save_ascii('schm1_conc.dat', answer)