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 [ ]:

    

A0, spring_const = 3, 5


    

In [ ]:

    

analytic_xmin = spring_const**-.5 * np.pi
def analytic(x):
    return np.array([+A0*np.cos(spring_const**.5*x),
                     -A0*np.sin(spring_const**.5*x)*spring_const**.5]).T


    

In [ ]:

    

integrate_kw = dict(integrator='cvode', method='adams', nsteps=2000)


    

In [ ]:

    

odesys = SymbolicSys.from_callback(dydt, 2, 1)
args = 17, [A0, 0], [spring_const]
res = odesys.integrate(*args, **integrate_kw)
res.plot()
print({k: v for k, v in res.info.items() if not k.startswith('internal_')})


    

In [ ]:

    

%load_ext autoreload
%autoreload 2


    

In [ ]:

    

from pyodesys.convergence import integrate_tolerance_series, fit_factory


    

In [ ]:

    

def plot_tolerance_series(tols, relative=False, **kwargs):
    result0, results, extra = integrate_tolerance_series(odesys, tols, None, *args, **integrate_kw, **kwargs)
    fig, axes = plt.subplots(2, 2, sharex=True, figsize=(14, 4), gridspec_kw=dict(hspace=0))
    ytru = analytic(result0.xout)
    errtru = np.abs(result0.yout - ytru)
    errest = extra['errest']
    if relative:
        errtru /= ytru
        errest /= result0.yout
    result0.plot(yerr=errest, ax=axes[0, 0])
    result0.plot(y=ytru, ax=axes[1, 0])
    result0.plot(y=errest, ax=axes[0, 1])
    result0.plot(y=errtru, ax=axes[1, 1])
    for ax in axes[:, 1]:
        ally = np.concatenate((errtru, extra['errest'])).flat
        if relative:
            ymax = np.max(np.sort(ally)[:-(1 + len(ally)//100)])
        else:
            ymax = np.max(ally)
        ax.set_ylim([0, 1.1*ymax])


    

In [ ]:

    

plot_tolerance_series(np.logspace(-6, -4, 7))


    

In [ ]:

    

plot_tolerance_series(np.logspace(-6, -4, 27))


    

In [ ]:

    

plot_tolerance_series(np.logspace(-10, -7, 50))


    

In [ ]:

    

plot_tolerance_series(np.logspace(-10, -8, 6))


    

In [ ]:

    

plot_tolerance_series(np.logspace(-10, -8, 12))


    

In [ ]:

    

plot_tolerance_series(np.logspace(-10, -8, 24))


    

In [ ]:

    

plot_tolerance_series(np.logspace(-10, -8, 6), fit=fit_factory())


    

In [ ]:

    

plot_tolerance_series(np.logspace(-10, -8, 10), relative=True)


    

In [ ]:

    

from pyodesys.plotting import plot_result


    

In [ ]:

    

plot_tolerance_series(np.logspace(-2.5, -1.5, 12), fit=fit_factory(3))


    

In [ ]: