Solving a transformed system of first order ordinary differential equations

In this notebook we explore how we can use different convenient subclasses of SymbolicSys to reformulate a problem (for other stability/accuracy/performance characteristics)



In [ ]:

    
from __future__ import print_function, division, absolute_import
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
from pyodesys import OdeSys
from pyodesys.symbolic import SymbolicSys, symmetricsys
sp.init_printing()
%matplotlib inline
print(sp.__version__)

We will consider a chain of three coupled decays:



In [ ]:

    
f = lambda t, x, k: [-k[0]*x[0], k[0]*x[0] - k[1]*x[1], k[1]*x[1] - k[2]*x[2], k[2]*x[2] - k[3]*x[3]]
k = [7, 3, 2, 0]
y0 = [1, 0, 0, 0]
tend = 1.0

The above equation has an analytic solution, evaluated at t = 1 it is:



In [ ]:

    
ref = [9.11881965554516244e-04, 8.55315762040415040e-02, 3.07983556726319885e-01, 6.05572985104084083e-01]
def rmsd(x):
    diff = np.array(ref[4-len(x):]) - x
    return np.sum(diff**2)**0.5

We can integrate this using e.g. lsoda from scipy:



In [ ]:

    
odesys = OdeSys(f)
res = odesys.integrate(np.linspace(0, tend), y0, params=k)
res.plot(names='abcd')
plt.legend()
{k: v for k, v in res.info.items() if not k.startswith('internal')}, rmsd(res.yout[-1, :])

For a stiff system (requiring an implicit method) we would have needed to define the jacobian too, then using symbolic manipulation can be very useful (less error prone and less user code needed to be written):



In [ ]:

    
odesys = SymbolicSys.from_callback(f, 4, 4, names='abcd')
odesys.exprs, odesys.get_jac(), odesys.get_dfdx()

Let's define a convienience function for plotting both with linear and logartithmic scale:



In [ ]:

    
def integrate_and_plot(odesys, integrator, tout, y0, k, interpolate=False, stiffness=False, **kwargs):
    plt.figure(figsize=(14,5))
    xout, yout, info = odesys.integrate(tout, y0, k, integrator=integrator, **kwargs)
    plt.subplot(1, 3 if stiffness else 2, 1)
    odesys.plot_result(interpolate=interpolate)
    plt.legend(loc='best')
    
    plt.subplot(1, 3 if stiffness else 2, 2)
    plt.gca().set_xscale('log')
    plt.gca().set_yscale('log')
    odesys.plot_result(interpolate=interpolate)
    plt.legend(loc='best')
    
    if stiffness:
        ratios = odesys.stiffness()
        plt.subplot(1, 3, 3)
        plt.yscale('linear')
        plt.plot(odesys._internal[0], ratios)
    info.pop('internal_xout')
    info.pop('internal_yout')
    return len(xout), info, rmsd(yout[-1, :])

Let's use the vode integrator this time and an implicit algorithm:



In [ ]:

    
integrate_and_plot(odesys, 'scipy', np.linspace(0, tend), y0, k, first_step=1e-14, name='vode', method='bdf')

We see that the jacobian was evaluated twice. The final error is slightly higher (which might be expected since the system is not particularly stiff)

SymbolicSys has provided us with callbacks if we manually want to evaluate f or its jacobian:



In [ ]:

    
print(odesys.f_cb(0, y0, k))
print()
print(odesys.j_cb(0, y0, k))

We can use the cvode integrator (through the use of pycvodes)



In [ ]:

    
integrate_and_plot(odesys, 'cvode', tend, y0, k, atol=1e-8, rtol=1e-8, method='bdf')

Since we are using symbolic manipulation it is very easy to perform a variable transformation. We will look at how this system behaves in logarithmic space (we need $y>0$ so we'll add a $\delta$ much smaller than our original linear absolute tolerance)

First we need to define our transformation, we use the helper function symmetricsys for this:



In [ ]:

    
y0aug = [1, 1e-20, 1e-20, 1e-20]
logexp = sp.log, sp.exp
LogLogSys = symmetricsys(
    logexp, logexp, exprs_process_cb=lambda exprs: [
    sp.powsimp(expr.expand(), force=True) for expr in exprs])

Alright, so our newly defined LogLogSys class can now take the same input and give back transformed expressions:



In [ ]:

    
tsys = LogLogSys.from_callback(f, 4, 4, names=True)
tsys.exprs, tsys.get_jac(), tsys.get_dfdx()

Now let us integrate the transformed system:



In [ ]:

    
integrate_and_plot(tsys, 'cvode', [1e-12, tend], y0aug, k, atol=1e-7, rtol=1e-7, first_step=1e-4, nsteps=5000)



In [ ]:

    
integrate_and_plot(tsys, 'odeint', [1e-12, tend], y0aug, k, atol=1e-8, rtol=1e-8, first_step=1e-4, nsteps=50000)

Substantially more work was needed to solve this transformed system, also the accuracy suffered. So this transformation was not of any help for this particular problem, choice of initial conditions and length of integration. There may be situations where it is useful though.

Stiffness



In [ ]:

    
def f2(t, x, p):
    y0 = p[0]*sp.exp(-p[1]*t)
    return [p[1]*y0 - p[2]*x[0]] + [
        p[i+2]*x[i] - p[i+3]*x[i+1] for i in range(len(x) - 1)
    ]
tsys2 = LogLogSys.from_callback(f2, 3, 5, names=True)
tsys2.exprs, tsys2.get_jac()



In [ ]:

    
integrate_and_plot(tsys2, 'cvode', [1e-12, tend], y0aug[1:], [y0aug[0], 7, 3, 2, 1],
                   atol=1e-7, rtol=1e-7, first_step=1e-4, stiffness=True, nsteps=5000)



In [ ]:

    
integrate_and_plot(odesys, 'cvode', [1e-12, tend], y0aug, [7, 3, 2, 1], atol=1e-7, rtol=1e-7, first_step=1e-4, stiffness=True)

Using PartiallySolvedSys

Sometimes we can solve a system paritally by integrating some of the dependent variables. The ODE system then needs to be reformulated in terms of the new anayltic function. PartiallySolvedSystem helps us with this task:



In [ ]:

    
from pyodesys.symbolic import PartiallySolvedSystem
psys = PartiallySolvedSystem(odesys, lambda x0, y0, p0: {odesys.dep[0]: y0[0]*sp.exp(-p0[0]*(odesys.indep-x0))})
psys.exprs, psys.get_jac()



In [ ]:

    
integrate_and_plot(psys, 'cvode', [1e-12, tend], y0aug, [7, 3, 2, 1], atol=1e-7, rtol=1e-7, first_step=1e-4, stiffness=True)

Using fixed RK4 stepper

pyodesys provides an example class of an integrator (RK4_example_integrator) in pyodesys.integrators. It can be used as a model integrator when designing custom steppers.



In [ ]:

    
from pyodesys.integrators import RK4_example_integrator



In [ ]:

    
integrate_and_plot(odesys, RK4_example_integrator, [1e-12, tend], y0aug, k, atol=1e-8, rtol=1e-8, first_step=1e-1)



In [ ]:

    
integrate_and_plot(psys, RK4_example_integrator, [1e-12, tend], y0aug, [7, 3, 2, 1], atol=1e-8, rtol=1e-8, first_step=1e-1)



In [ ]: