```python
import numpy as np
import sympy
import math
import matplotlib.pyplot as plt
from pyodesys.symbolic import SymbolicSys
sympy.init_printing()
%matplotlib inline
```


```python
# 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
    )
```


```python
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
```


```python
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)
```


```python
compare_autonomous(1)
```


```python
compare_autonomous(1000)
```


```python
from pyodesys import chained_parameter_variation
```


```python
odesys_T = get_odesys(T_as_param=True)
odesys_T.autonomous_interface and odesys_T.autonomous_exprs
```


```python
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)
```


```python
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])
```


```python
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)
```


```python
odesys_by_name_T = get_odesys(T_as_param=True, dep_by_name=True, par_by_name=True)
```


```python
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)
```


```python
plot_results(results)
```


```python
odesys_by_name = get_odesys(dep_by_name=True, par_by_name=True)
plot_results(odesys_by_name.integrate(tout, ic, ps, **ikw))
```


```python

```

6cb90b99e77b6610de543ff4e870b8de677bda78
ipynb
Jupyter Notebook
examples/_convert_to_autonomous_eyring_sinusoidal.ipynb
slayoo/pyodesys
8e1afb195dadf6c6f8e765873bc9dd0fae067c39
BSD-2-Clause
2015-09-29T16:51:03.000Z
2022-02-02T13:26:50.000Z
examples/_convert_to_autonomous_eyring_sinusoidal.ipynb
slayoo/pyodesys
8e1afb195dadf6c6f8e765873bc9dd0fae067c39
BSD-2-Clause
2015-09-29T14:40:45.000Z
2021-09-18T19:29:50.000Z
examples/_convert_to_autonomous_eyring_sinusoidal.ipynb
slayoo/pyodesys
8e1afb195dadf6c6f8e765873bc9dd0fae067c39
BSD-2-Clause
2016-03-18T14:00:39.000Z
2021-09-17T13:54:29.000Z
Qwen/Qwen-72B
1. YES
2. YES