In [ ]:

    

%load_ext autoreload
%autoreload 2


    

In [ ]:

    

import numpy as np
import matplotlib.pyplot as plt
from pyodesys.symbolic import SymbolicSys, get_logexp, symmetricsys, PartiallySolvedSystem
import sympy as sp
sp.init_printing()
%matplotlib inline


    

In [ ]:

    

l3sys = SymbolicSys.from_callback(lambda t, y, p: [-p[0]*y[0], p[0]*y[0] - p[1]*y[1], p[1]*y[1]],
                                   3, 2, linear_invariants=[[1, 1, 1]])
l3sys.exprs, l3sys.linear_invariants, l3sys.nonlinear_invariants or []


    

In [ ]:

    

logexp = get_logexp()
LogLogSys = symmetricsys(logexp, logexp)
def expr_proc(exprs):
    return [sp.powsimp(expr.expand(), force=True).simplify() for expr in exprs]

t3sys = LogLogSys.from_other(l3sys, exprs_process_cb=expr_proc)
t3sys.exprs, t3sys.linear_invariants or [], t3sys.nonlinear_invariants


    

In [ ]:

    

l2sys = PartiallySolvedSystem.from_linear_invariants(l3sys)
l2sys.exprs, l2sys.linear_invariants or [], l2sys.nonlinear_invariants or []


    

In [ ]:

    

t2sys = LogLogSys.from_other(l2sys, exprs_process_cb=expr_proc)
t2sys.exprs, l2sys.linear_invariants or [], l2sys.nonlinear_invariants or []


    

In [ ]:

    

y0, y1, y2 = l3sys.dep
p0, p1 = l3sys.params
x = l3sys.indep
exp = l3sys.be.exp
assert l3sys.exprs[0] + p0*y0 == 0
assert l3sys.exprs[1] - p0*y0 + p1*y1 == 0
assert l3sys.exprs[2] - p1*y1 == 0
assert l3sys.linear_invariants.tolist() == [[1, 1, 1]]

assert t3sys.exprs[0] + p0*exp(x) == 0
assert t3sys.exprs[1] + p1*exp(x) - p0*exp(x+y0-y1) == 0
assert t3sys.exprs[2] - p1*exp(x+y1-y2) == 0
assert t3sys.nonlinear_invariants[0] - exp(y0) - exp(y1) - exp(y2) == 0

i0, i1, i2 = l2sys.init_dep
assert l2sys.exprs[0] - p0*(i0 + i1 + i2 - y1 - y2) + p1*y1 == 0
assert l2sys.exprs[1] - p1*y1 == 0

assert t2sys.exprs[0] - exp(x-y1)*(p0*(i0 + i1 + i2 - exp(y1) - exp(y2)) - p1*exp(y1))


    

In [ ]:

    

l2sys.init_dep


    

In [ ]:

    

params = [2.05, 2.4]
results = [odesys.integrate(
    [1e-8, 3.5], [1, 1e-16, 1e-16], params, integrator='cvode', nsteps=2000,
    atol=1e-9, rtol=1e-9) for odesys in [l3sys, l2sys, t3sys, t2sys]]
fig, axes = plt.subplots(1, 4, figsize=(16, 4))
for res, ax in zip(results, axes):
    res.plot(ax=ax)


    

In [ ]:

    

from pyodesys.tests.bateman import bateman_full
for res in results:
    ref = bateman_full(res.yout[0, :], params + [0], res.xout - res.xout[0], exp=np.exp)
    assert np.allclose(res.yout, np.array(ref).T)


    

In [ ]:

    

logexp16 = get_logexp(1, 10**-16)
LogLog16 = symmetricsys(logexp16, logexp16)
t3sys16 = LogLog16.from_other(l3sys, lower_bounds=[0]*l3sys.ny, upper_bounds=[1+1e-8]*l3sys.ny)
assert np.all(t3sys16.lower_bounds < -30) and np.all(t3sys16.upper_bounds > 0)


    

In [ ]: