In [ ]:

    

from itertools import product, repeat
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import RadioButtons
from pyodesys import ODESys
from pyodesys.integrators import (RK4_example_integrator, EulerForward_example_integrator,
                                  EulerBackward_example_integrator, Trapezoidal_example_integrator,
                                  Midpoint_example_integrator)
from pyodesys.tests.bateman import bateman_full
%matplotlib inline


    

In [ ]:

    

choice = RadioButtons(
    options=['decay', 'sine'],
    description='Case:',
    disabled=False
)
choice


    

In [ ]:

    

if choice.value == 'decay':
    odesys = ODESys(lambda x, y, p: [
        (0 if i == 0 else p[i-1]*y[i-1]) -
        (0 if i == len(y) - 1 else p[i]*y[i])
         for i in range(len(y))
    ], lambda x, y, p: [[(-p[ri] if ri == ci and ri < len(y) - 1 else (p[ci] if ri == ci + 1 else 0)) for ci in range(len(y))] for ri in range(len(y))])
    x0, xend, y0, p, n0 = 0, 1, [5, 2, .3, .1, 0], [13, 7, 3, 1], 24

    def ref(x):
        return np.array(bateman_full(y0, p+[0], x, exp=np.exp)).T
elif choice.value == 'sine':
    A = 3
    k = 14
    x0, xend, y0, p, n0 = 0, 1, [0, A*k], [k], 128
    odesys = ODESys(lambda x, y, p: [y[1], -p[0]**2*y[0]],
                   lambda x, y, p: [[0, 1], [-p[0]**2, 0]])
    def ref(x):
        return np.array([A*np.sin(k*x), A*k*np.cos(k*x)]).T
else:
    raise ValueError("Unknown choice")


    

In [ ]:

    

xout1 = np.linspace(x0, xend, n0 + 1)
res1 = odesys.integrate(xout1, y0, p, first_step=xend/n0, integrator=EulerForward_example_integrator)
yref1 = ref(res1.xout)
fig, axes = plt.subplots(1, 2, figsize=(14, 4))
res1.plot(ls=('-',), ax=axes[0])
xplt = np.linspace(x0, xend)
res1.plot(x=xplt, y=ref(xplt), ls=(':',), ax=axes[0])
res1.plot(y=res1.yout - yref1, ax=axes[1])


    

In [ ]:

    

def plot_order_of_convergence(Integrator, nruns=6):
    res, nruns = [], 7
    for i in range(nruns):
        n = n0 * 2**i
        x = np.linspace(0, 1, n)
        res.append(odesys.integrate(x, y0, p, integrator=Integrator))

    vals = {}
    for ix, iy in product(range(n0), range(len(y0))):
        vals[ix, iy] = [res[i].yout[ix*(2**i), iy] for i in range(nruns)]

    logn = np.log(n0 * 2**np.arange(nruns))
    fig, axes = plt.subplots(1, len(y0), figsize=(14, 4))
    for iy, ax in enumerate(axes):
        for ix in range(1, n0):
            lnabserr = np.log(np.abs(np.array(vals[ix, iy]) - yref1[ix, iy]))
            ax.plot(logn, lnabserr, '.-', c=(ix/n0, 0, 0))
        popt = np.polyfit(logn, lnabserr, 1)
        ax.plot(logn[[0,-1]], np.polyval(popt, logn[[0,-1]]), label='y=%.2f x %+.2f' % tuple(popt))
        ax.legend()


    

In [ ]:

    

plot_order_of_convergence(EulerForward_example_integrator)


    

In [ ]:

    

plot_order_of_convergence(EulerBackward_example_integrator)


    

In [ ]:

    

plot_order_of_convergence(Midpoint_example_integrator)


    

In [ ]:

    

plot_order_of_convergence(Trapezoidal_example_integrator)


    

In [ ]:

    

plot_order_of_convergence(RK4_example_integrator)