In [ ]:

    

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


    

In [ ]:

    

def mk_odesys(power=1):
    def rhs(t, y, p):
        k = t**power
        return [-k*y[0], k*y[0]]
    
    def analytic(tout, init_y):
        y0ref = init_y[0]*np.exp(-tout**(power+1)/(power+1))
        return np.array([y0ref, init_y[0] - y0ref + init_y[1]]).T
    
    odesys = SymbolicSys.from_callback(rhs, 2, 0, linear_invariants=[[1, 1]])
    return odesys, analytic


    

In [ ]:

    

odesys1, analytic1 = mk_odesys(1)
odesys1.dep, odesys1.exprs, odesys1.get_jac()


    

In [ ]:

    

def integrate_and_plot(odesys, analytic, axes=None, tout=3, init_y=(1,0)):
    result = odesys.integrate(tout, init_y, integrator='cvode')
    if axes is None:
        _fig, axes = plt.subplots(1, 3, figsize=(14, 4))
    result.plot(ax=axes[0])
    result.plot(ax=axes[1], y=result.yout - analytic(result.xout, result.yout[0, :]))
    result.plot_invariant_violations(ax=axes[2])
    print({k: v for k, v in result.info.items() if not k.startswith('internal')})


    

In [ ]:

    

integrate_and_plot(odesys1, analytic1)


    

In [ ]:

    

asys1 = odesys1.as_autonomous()
asys1.dep, asys1.exprs, asys1.get_jac()


    

In [ ]:

    

integrate_and_plot(asys1, analytic1)


    

In [ ]:

    

def compare(power, **kwargs):
    odesys, analytic = mk_odesys(power)
    autsys = odesys.as_autonomous()
    fig, axes = plt.subplots(2, 3, figsize=(14, 8))
    integrate_and_plot(odesys, analytic, axes=axes[0], **kwargs)
    integrate_and_plot(autsys, analytic, axes=axes[1], **kwargs)


    

In [ ]:

    

compare(2)


    

In [ ]:

    

compare(3, init_y=(3,1))


    

In [ ]:

    

compare(10, tout=1.5)