Figure 5.8:

Concentrations of B and C versus time for full model with increasing k_2 with K_2=k_2/k_{-2}=1.

Code for Figure 5.8

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
# Converted from rls_fast.m
# Full model (nscheme=0) for k2=k_2 = 10,20,50,100 at short time
import numpy as np
from scipy.integrate import solve_ivp
from misc import save_ascii

k1   = 1.0
k_1  = 0.5
c_A0 = 1.0
c_B0 = 0.4
c_C0 = 0.0
x0   = np.array([c_A0, c_B0, c_C0])
nscheme = 0

def rhs(t, x, 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
    return np.array([-r1, r1 - r2, r2])

tfinal = 0.25
ntimes = 100
sqeps = np.sqrt(np.finfo(float).eps)

table = None
for k2val in [10., 20., 50., 100.]:
    tout = np.linspace(0, tfinal, ntimes)
    sol = solve_ivp(lambda t, x: rhs(t, x, k2val, k2val), [0, tfinal], x0,
                    method='Radau', t_eval=tout, rtol=sqeps, atol=sqeps)
    if table is None:
        table = np.column_stack([tout, sol.y.T])
    else:
        table = np.column_stack([table, sol.y.T])

save_ascii('rls_fast.dat', table)