In [ ]:

    

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
from pyodesys.tests._robertson import get_ode_exprs
from pyodesys.symbolic import ScaledSys, PartiallySolvedSystem
sp.init_printing()
%matplotlib inline


    

In [ ]:

    

linf, linj = get_ode_exprs()
logf, logj = get_ode_exprs(True, True)

linsys = ScaledSys.from_callback(linf, 3, 3, dep_scaling=1e8, linear_invariants=[[1, 1, 1]], names='ABC')


    

In [ ]:

    

psysA = PartiallySolvedSystem.from_linear_invariants(linsys, preferred=[0], description='A ')
psysC = PartiallySolvedSystem.from_linear_invariants(linsys, preferred=[2], description='C ')
psysA.exprs, psysC.exprs


    

In [ ]:

    

tend, iv, pars = 1e18, [1, 0, 0], [0.04, 1e4, 3e7]
integrate_kw = dict(integrator='cvode', record_rhs_xvals=True, record_jac_xvals=True, nsteps=3000,
                    atol=1e-8, rtol=1e-8, return_on_error=True)

results = []
fig, axes = plt.subplots(1, 3, figsize=(16, 4))
def plot_res(res, ax):
    res.plot(xscale='log', yscale='log', ax=ax)  #, info_vlines_kw=True)
    ax.set_title((res.odesys.description or '') + 
                 ('%d steps, ' % res.info['n_steps']) +
                 ('success' if res.info['success'] else 'failed'))
    ax.set_xlim([1e-12, tend])
    ax.set_ylim([1e-30, 1e9])
    ax.legend(loc='best')

for odesys, ax in zip([linsys, psysA, psysC], axes):
    results.append(odesys.integrate(tend, iv, pars, **integrate_kw))
    plot_res(results[-1], ax)


    

In [ ]:

    

linsys.autonomous_exprs, psysA.autonomous_exprs, psysC.autonomous_exprs


    

In [ ]:

    

linsys.autonomous_interface, psysA.autonomous_interface, psysC.autonomous_interface


    

In [ ]:

    

psysAr = PartiallySolvedSystem.from_linear_invariants(linsys, preferred=[0], roots=[1000*linsys.dep[0] - linsys.dep[2]])
psysAr.autonomous_exprs, psysAr.autonomous_interface


    

In [ ]:

    

resAroot = psysAr.integrate(tend, iv, pars, return_on_root=True, **integrate_kw)
resAswitch = resAroot.copy()
resAswitch.extend_by_integration(tend, odesys=psysC, **integrate_kw)
fig, axes = plt.subplots(1, 2, figsize=(16, 4))
axes[1].axvline(resAswitch.xout[resAswitch.info['root_indices'][0]], ls='--', alpha=.5)
for res, ax in zip([resAroot, resAswitch], axes):
    plot_res(res, ax)