IRREVERSIBLE unimolecular reaction `A -> B`,¶

with 1st-order kinetics.

Adaptive variable time steps compared with exact analytical solution

TAGS : "uniform compartment", "chemistry", "numerical", "quick-start", "basic", "under-the-hood"¶

In [1]:

LAST_REVISED = "Oct. 6, 2024"
LIFE123_VERSION = "1.0.0.beta.39"   # Library version this experiment is based on

In [2]:

#import set_path                    # Using MyBinder?  Uncomment this before running the next cell!

In [3]:

#import sys
#sys.path.append("C:/some_path/my_env_or_install")   # CHANGE to the folder containing your venv or libraries installation!
# NOTE: If any of the imports below can't find a module, uncomment the lines above, or try:  import set_path   

import ipynbname
import numpy as np

from life123 import check_version, UniformCompartment, ReactionKinetics, GraphicLog, PlotlyHelper

In [4]:

check_version(LIFE123_VERSION)

OK

In [5]:

# Initialize the HTML logging (for the graphics)
log_file = ipynbname.name() + ".log.htm"    # Use the notebook base filename for the log file
                                            # IN CASE OF PROBLEMS, set manually to any desired name

# Set up the use of some specified graphic (Vue) components
GraphicLog.config(filename=log_file,
                  components=["vue_cytoscape_2"],
                  extra_js="https://cdnjs.cloudflare.com/ajax/libs/cytoscape/3.21.2/cytoscape.umd.js")

-> Output will be LOGGED into the file 'react_2_d.log.htm'

In [ ]:

PART 1 - VARIABLE TIME STEPS (Numerical Approximation to Solution)¶

In [6]:

# Instantiate the simulator and specify the chemicals
# Here we use the "fast" preset for the variable steps, trying to push the envelope on speed
uc = UniformCompartment(preset="fast")

# Reaction A <-> B , with 1st-order kinetics in both directions
uc.add_reaction(reactants="A", products="B", 
                forward_rate=3., reverse_rate=0)      # Notice the zero reverse rate (irreversible)

uc.describe_reactions()

Number of reactions: 1 (at temp. 25 C)
0: A <-> B  (kF = 3 / kR = 0) | 1st order in all reactants & products
Set of chemicals involved in the above reactions: {'B', 'A'}

In [7]:

# Send a plot of the network of reactions to the HTML log file
uc.plot_reaction_network("vue_cytoscape_2")

[GRAPHIC ELEMENT SENT TO LOG FILE `react_2_d.log.htm`]

In [8]:

# Set the initial concentrations of all the chemicals
uc.set_conc({"A": 50., "B": 10.})
uc.describe_state()

SYSTEM STATE at Time t = 0:
2 species:
  Species 0 (A). Conc: 50.0
  Species 1 (B). Conc: 10.0
Set of chemicals involved in reactions: {'B', 'A'}

In [9]:

uc.get_history()

Out[9]:

	SYSTEM TIME	A	B	caption
0	0.0	50.0	10.0	Initialized state

In [10]:

uc.enable_diagnostics()   # To save diagnostic information about the call to single_compartment_react()
                          # Useful for insight into the inner workings of the simulation

Run the reaction¶

Passing True (default) to variable_steps automatically adjusts up or down the time steps¶

In [11]:

uc.single_compartment_react(initial_step=0.1, target_end_time=1.5,
                            variable_steps=True)

Some steps were backtracked and re-done, to prevent negative concentrations or excessively large concentration changes
46 total step(s) taken
Number of step re-do's because of elective soft aborts: 5
Norm usage: {'norm_A': 24, 'norm_B': 10, 'norm_C': 10, 'norm_D': 10}

In [12]:

history = uc.get_history()   # The system's history, saved during the run of single_compartment_react()
history

Out[12]:

	SYSTEM TIME	A	B	caption
0	0.000000	50.000000	10.000000	Initialized state
1	0.007776	48.833600	11.166400	1st reaction step
2	0.019440	47.124815	12.875185
3	0.028771	45.805621	14.194379
4	0.038102	44.523357	15.476643
5	0.047434	43.276988	16.723012
6	0.061430	41.459770	18.540230
7	0.072628	40.067040	19.932960
8	0.083825	38.721095	21.278905
9	0.095023	37.420364	22.579636
10	0.106220	36.163327	23.836673
11	0.123016	34.341112	25.658888
12	0.136453	32.956795	27.043205
13	0.149890	31.628281	28.371719
14	0.163327	30.353320	29.646680
15	0.176764	29.129754	30.870246
16	0.196919	27.368389	32.631611
17	0.213044	26.044499	33.955501
18	0.229168	24.784650	35.215350
19	0.253355	22.986291	37.013709
20	0.272704	21.651993	38.348007
21	0.292053	20.395149	39.604851
22	0.321077	18.619317	41.380683
23	0.344296	17.322350	42.677650
24	0.367515	16.115727	43.884273
25	0.402343	14.431866	45.568134
26	0.430206	13.225529	46.774471
27	0.472000	11.567277	48.432723
28	0.505436	10.407008	49.592992
29	0.555589	8.841178	51.158822
30	0.595711	7.776989	52.223011
31	0.655895	6.372845	53.627155
32	0.716078	5.222222	54.777778
33	0.776262	4.279344	55.720656
34	0.836446	3.506704	56.493296
35	0.896629	2.873565	57.126435
36	0.956813	2.354740	57.645260
37	1.016997	1.929589	58.070411
38	1.077181	1.581200	58.418800
39	1.137364	1.295713	58.704287
40	1.197548	1.061770	58.938230
41	1.257732	0.870067	59.129933
42	1.317915	0.712975	59.287025
43	1.378099	0.584247	59.415753
44	1.438283	0.478760	59.521240
45	1.498466	0.392320	59.607680
46	1.558650	0.321486	59.678514	last reaction step

In [13]:

uc.plot_history(colors=['darkturquoise', 'orange'], show_intervals=True)    # Plots of concentration with time

Notice how the reaction proceeds in smaller steps in the early times, when [A] and [B] are changing much more rapidly¶

That resulted from passing the flag variable_steps=True to single_compartment_react()

In [14]:

uc.plot_step_sizes(show_intervals=True)        #  To see the sizes of the steps taken

Why the zigzag? It's because of the "fast" preset picked for the variable steps, in the instantiation of the class UniformCompartment: it's like a "high-strung driver" that tries to get away with excessive speeed - and periodically overdoes on acceleration, and then slams on the brakes!
Other presets (such as "mid") are more "mild-mannered" and more conservative about going too fast too soon.

In [ ]:

PART 2 - Comparison with exact analytical solution¶

In [15]:

t_arr = np.linspace(0., 1.5, 50)   # A relatively dense uniform grid across our time range
t_arr

Out[15]:

array([0.        , 0.03061224, 0.06122449, 0.09183673, 0.12244898,
       0.15306122, 0.18367347, 0.21428571, 0.24489796, 0.2755102 ,
       0.30612245, 0.33673469, 0.36734694, 0.39795918, 0.42857143,
       0.45918367, 0.48979592, 0.52040816, 0.55102041, 0.58163265,
       0.6122449 , 0.64285714, 0.67346939, 0.70408163, 0.73469388,
       0.76530612, 0.79591837, 0.82653061, 0.85714286, 0.8877551 ,
       0.91836735, 0.94897959, 0.97959184, 1.01020408, 1.04081633,
       1.07142857, 1.10204082, 1.13265306, 1.16326531, 1.19387755,
       1.2244898 , 1.25510204, 1.28571429, 1.31632653, 1.34693878,
       1.37755102, 1.40816327, 1.43877551, 1.46938776, 1.5       ])

In [16]:

# The exact solution is available for a simple scenario like the one we're simulating here

A_exact, B_exact = ReactionKinetics.exact_solution_unimolecular_irreversible(kF=3., A0=50., B0=10., t_arr=t_arr)

In [17]:

fig_exact = PlotlyHelper.plot_curves(x=t_arr, y=[A_exact, B_exact], title="EXACT solution", x_label="SYSTEM TIME", y_label="concentration",
                                     legend_title="Chemical", curve_labels=["A (EXACT)", "B (EXACT)"],
                                     colors=["darkturquoise", "orange"], show=True)

To avoid clutter, we'll just plot [A], as obtained from the variable-step approx solution and the exact analytical solution¶

In [18]:

fig_exact = PlotlyHelper.plot_curves(x=t_arr, y=A_exact, title="EXACT solution", x_label="SYSTEM TIME", y_label="concentration",
                                     curve_labels="A (EXACT)", legend_title="Chemical",
                                     colors="red", show=True)     # Repeat a portion of the diagram seen just before

In [19]:

fig_variable = uc.plot_history(chemicals=['A'], colors='darkturquoise', title="VARIABLE time steps", show=True)     # Repeat a portion of the diagram seen in Part 1

In [20]:

PlotlyHelper.combine_plots(fig_list=[fig_variable, fig_exact],
                           xrange=[0, 1.5], ylabel="concentration [A]",
                           title="Variable time steps vs. Exact soln, for [A] in irreversible reaction `A->B`",
                           legend_title="Simulation run")    # Both plots put together

A pretty good overlap!¶

In [ ]:

IRREVERSIBLE unimolecular reaction A -> B,¶

TAGS : "uniform compartment", "chemistry", "numerical", "quick-start", "basic", "under-the-hood"¶

PART 1 - VARIABLE TIME STEPS (Numerical Approximation to Solution)¶

Run the reaction¶

Passing True (default) to variable_steps automatically adjusts up or down the time steps¶

Notice how the reaction proceeds in smaller steps in the early times, when [A] and [B] are changing much more rapidly¶

PART 2 - Comparison with exact analytical solution¶

To avoid clutter, we'll just plot [A], as obtained from the variable-step approx solution and the exact analytical solution¶

A pretty good overlap!¶

IRREVERSIBLE unimolecular reaction `A -> B`,¶