In [ ]:

    

import numpy as np
import matplotlib.pyplot as plt
from pyodesys.symbolic import SymbolicSys
%matplotlib inline


    

In [ ]:

    

def dydt(t, y, params, backend):
    return [y[1], -y[0]*params[0]]  # y'' = -k y => y = C1 + sin(k**0.5 * x + C2)


    

In [ ]:

    

odesys = SymbolicSys.from_callback(dydt, 2, 1)
odesys.nonlinear_invariants = [odesys.params[0]*odesys.dep[0]**2/2 + odesys.dep[1]**2/2]
A0, spring_const = 3, 5
res = odesys.integrate(17, [A0, 0], [spring_const], integrator='cvode', method='adams', nsteps=2000)
res.plot()
print({k: v for k, v in res.info.items() if not k.startswith('internal_')})


    

In [ ]:

    

res.plot(deriv=True)


    

In [ ]:

    

analytic_xmin = spring_const**-.5 * np.pi
def analytic(x):
    return A0*np.cos(spring_const**.5*x)


    

In [ ]:

    

fig, axes = plt.subplots(1, 2, figsize=(12, 3))
lo, hi = 1.27, 1.52
xplt = np.linspace(lo, hi)
axes[0].plot(xplt, analytic(xplt), ls=':')
res.plot(ax=axes[0], between=(lo, hi), indices=(0,))
axes[0].axvline(analytic_xmin, 0, 1, linestyle='--') #, transform=axes[1].get_xaxis_transform())
axes[0].axvline(res.xout[np.argmin(res.yout[res.xout<2, 0])], 0, 1, c='k') # transform=axes[1].get_xaxis_transform(), 
x_select, y_select = res.between(lo, hi)
axes[1].plot(x_select, y_select[:, 0] - analytic(x_select))
plt.tight_layout()


    

In [ ]:

    

fig, axes = plt.subplots(2, 2, figsize=(12, 6))
x_at = np.linspace(lo, hi)

for ri, use_deriv in enumerate([False, True]):
    y_at = np.array([res.at(_, use_deriv=use_deriv)[0][0] for _ in x_at])
    found_xmin = x_at[np.argmin(y_at[x_at<2])]
    axes[ri][0].plot(x_at, y_at, 'k')
    axes[ri][0].plot(x_at, analytic(x_at), ls=':')
    axes[ri][0].axvline(analytic_xmin, 0, 1, linestyle='--')
    axes[ri][0].axvline(found_xmin, 0, 1, c='k')
    axes[ri][0].set_title('min_loc_err: %.5g' % (found_xmin - analytic_xmin))
    axes[ri][1].plot(x_at, y_at - analytic(x_at))


    

In [ ]:

    

_ = res.plot_invariant_violations()


    

In [ ]:

    

import matplotlib.pyplot as plt
odesys.nonlinear_invariant_names = ['energy']
_ = res.plot_invariant_violations()
_ = plt.legend(loc='best')