In [ ]:
import numpy as np
import sympy
import math
import matplotlib.pyplot as plt
from pyodesys.symbolic import SymbolicSys
sympy.init_printing()
%matplotlib inline
In [ ]:
# 2 HNO2 -> H2O + NO + NO2; MassAction(EyringHS.fk('dH1', 'dS1'))
# 2 NO2 -> N2O4; MassAction(EyringHS.fk('dH2', 'dS2'))
#
# HNO2 H2O NO NO2 N2O4
def get_odesys(scaling=1, T_as_param=False, par_by_name=False, dep_by_name=False, **kwargs):
names = 'HNO2 H2O NO NO2 N2O4'.split()
param_names = 'dH1 dS1 dH2 dS2 T'.split()
def rhs(t, y, p, backend=math):
HNO2, H2O, NO, NO2, N2O4 = [y[k] for k in (names if dep_by_name else range(5))]
dH1, dS1, dH2, dS2 = [p[k] for k in (param_names[:-1] if par_by_name else range(4))]
R = 8.314
if T_as_param:
T = p['T' if par_by_name else 4]
else:
T = 300 + 10*backend.sin(0.2*math.pi*t - math.pi/2)
kB_h = 2.08366e10
k1 = kB_h*T*backend.exp(dS1/R - dH1/(R*T))/scaling # bimolecular => scaling**-1
k2 = kB_h*T*backend.exp(dS2/R - dH2/(R*T))/scaling # bimolecular => scaling**-1
r1 = k1*HNO2**2
r2 = k2*NO2**2
exprs = [-2*r1, r1, r1, r1 - 2*r2, r2]
return dict(zip(names, exprs)) if dep_by_name else exprs
return SymbolicSys.from_callback(
rhs, 5, 4 + int(T_as_param), names=names,
param_names=param_names[slice(None) if T_as_param else slice(0,-1)],
dep_by_name=dep_by_name, par_by_name=par_by_name, **kwargs
)
In [ ]:
def integrate_and_plot(system, init_y, p, axes=None):
result = system.integrate(np.linspace(0, 60, 200), init_y, p, integrator='cvode', nsteps=5000)
if axes is None:
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
result.plot(ax=axes[0], indices=[0,2,4], title_info=True)
result.plot(ax=axes[1], indices=[3])
print({k: v for k, v in sorted(result.info.items()) if not k.startswith('internal')})
return result
In [ ]:
def compare_autonomous(scaling):
fig, all_axes = plt.subplots(3, 2, figsize=(10, 14))
odesys = get_odesys(scaling)
autsys = odesys.as_autonomous()
assert all(o - a == 0 for o, a in zip(odesys.exprs, autsys.exprs[:-1]))
assert autsys.exprs[-1] == 1
init_y = [1*scaling, 55*scaling, 0, 0, 0]
p = [85e3, 10, 70e3, 20]
argso = odesys.pre_process(*odesys.to_arrays(0, init_y, p))
argsa = autsys.pre_process(*autsys.to_arrays(0, init_y, p))
for idx, (argo, arga) in enumerate(zip(argso, argsa)):
assert np.allclose(argo, arga[:-1] if idx == 1 else arga)
fo = odesys.f_cb(*argso)
fa = autsys.f_cb(*argsa)
assert np.allclose(fo, fa[:, :-1])
assert np.allclose(fa[:, -1], 1)
res1 = integrate_and_plot(odesys, init_y, p, axes=all_axes[0, :])
res2 = integrate_and_plot(autsys, init_y, p, axes=all_axes[1, :])
assert np.allclose(res1.yout, res2.yout, atol=1e-6)
diff = res1.yout-res2.yout
res1.plot(ax=all_axes[2, 0], indices=[0,2,4], y=diff/scaling)
res1.plot(ax=all_axes[2, 1], indices=[3], y=diff/scaling)
In [ ]:
compare_autonomous(1)
In [ ]:
compare_autonomous(1000)
In [ ]:
from pyodesys import chained_parameter_variation
In [ ]:
odesys_T = get_odesys(T_as_param=True)
odesys_T.autonomous_interface and odesys_T.autonomous_exprs
In [ ]:
def integrate_and_plot_piecewise(N, axes=None):
tout = np.linspace(0, 60, N)
tcent = tout[:-1] + np.diff(tout)/2
T = 300 + 10*np.sin(0.2*np.pi*tcent - np.pi/2)
init_y = [1, 55, 0, 0, 0]
p = [85e3, 10, 70e3, 20, 300]
res = chained_parameter_variation(odesys_T, np.diff(tout), init_y, {4: T}, p,
integrate_kwargs={'integrator': 'cvode'})
if axes is None:
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
res.plot(ax=axes[0], indices=[0,2,4], title_info=True)
res.plot(ax=axes[1], indices=[3])
axes[2].plot(tcent, T)
In [ ]:
fig, axes = plt.subplots(3, 3, figsize=(14, 14))
integrate_and_plot_piecewise(10, axes=axes[0])
integrate_and_plot_piecewise(30, axes=axes[1])
integrate_and_plot_piecewise(90, axes=axes[2])
In [ ]:
def plot_results(results):
for axes, res in zip(plt.subplots(len(results), 2, figsize=(11, 11))[1], results):
res.plot(ax=axes[0], indices=[0,2,4], title_info=True)
res.plot(ax=axes[1], indices=[3], title_info=True)
In [ ]:
odesys_by_name_T = get_odesys(T_as_param=True, dep_by_name=True, par_by_name=True)
In [ ]:
tout = np.linspace(0, 60, 100)
tcent = tout[:-1] + np.diff(tout)/2
T = 300 + 10*np.sin(0.2*np.pi*tcent - np.pi/2)
H1vals = [80e3, 85e3, 90e3]
ps = dict(dH1=H1vals, dS1=10, dH2=70e3, dS2=20)
ic = dict(HNO2=1, H2O=55, NO=0, NO2=0, N2O4=0)
ikw = dict(integrator='cvode')
results = chained_parameter_variation(odesys_by_name_T, np.diff(tout), ic, dict(T=T), ps,
integrate_kwargs=ikw)
In [ ]:
plot_results(results)
In [ ]:
odesys_by_name = get_odesys(dep_by_name=True, par_by_name=True)
plot_results(odesys_by_name.integrate(tout, ic, ps, **ikw))
In [ ]: