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 SymbolicSys, 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)

def plot_res(res, ax):
    res.plot(xscale='log', yscale='log', ax=ax, title_info=2)  # , info_vlines_kw=True)
    ax.set_xlim([1e-12, tend])
    ax.set_ylim([1e-30, 1e9])
    ax.legend(loc='best')

def integrate_and_plot_systems(systems):
    results = []
    fig, axes = plt.subplots(1, 3, figsize=(16, 4))
    for odesys, ax in zip(systems, axes):
        results.append(odesys.integrate(tend, iv, pars, **integrate_kw))
        plot_res(results[-1], ax)
    return results

odes = [linsys, psysA, psysC]


    

In [ ]:

    

results = integrate_and_plot_systems(odes)


    

In [ ]:

    

_ = integrate_and_plot_systems([SymbolicSys.from_other(odesys, lower_bounds=[0]*odesys.ny) for odesys in odes])


    

In [ ]: