Units in ODEs

| categories: odes, units | tags: | View Comments

We reconsider a simple ODE but this time with units. We will use the quantities package again.

Here is the ODE, \(\frac{dCa}{dt} = -k Ca\) with \(C_A(0) = 1.0\) mol/L and \(k = 0.23\) 1/s. Compute the concentration after 5 s.

import quantities as u

k = 0.23 / u.s
Ca0 = 1 * u.mol / u.L

def dCadt(Ca, t):
    return -k * Ca

import numpy as np
from scipy.integrate import odeint

tspan = np.linspace(0, 5) * u.s

sol = odeint(dCadt, Ca0, tspan)

print sol[-1]
[ 0.31663678]

No surprise, the units are lost. Now we start wrapping odeint. We wrap everything, and then test two examples including a single ODE, and a coupled set of ODEs with mixed units.

import quantities as u
import matplotlib.pyplot as plt

import numpy as np
from scipy.integrate import odeint as _odeint

def odeint(func, y0, t, args=(),
           Dfun=None, col_deriv=0, full_output=0,
           ml=None, mu=None, rtol=None, atol=None,
           tcrit=None, h0=0.0, hmax=0.0, hmin=0.0,
           ixpr=0, mxstep=0, mxhnil=0, mxordn=12,
           mxords=5, printmessg=0):

    def wrapped_func(Y0, T, *args):
        # put units on T if they are on the original t
        # check for units so we don't put them on twice
        if not hasattr(T, 'units') and hasattr(t, 'units'):
            T = T * t.units
        # now for the dependent variable units. Y0 may be a scalar or
        # a list or an array. we want to check each element of y0 for
        # units, and add them to the corresponding element of Y0 if we
        # need to.
        try:
            uY0 = [x for x in Y0] # a list copy of contents of Y0
            # this works if y0 is an iterable, eg. a list or array
            for i, yi in enumerate(y0):
                if not hasattr(uY0[i],'units') and hasattr(yi, 'units'):
               
                    uY0[i] = uY0[i] * yi.units
                
        except TypeError:
            # we have a scalar
            if not hasattr(Y0, 'units') and hasattr(y0, 'units'):
                uY0 = Y0 * y0.units
       
        val = func(uY0, t, *args)

        try:
            return np.array([float(x) for x in val])
        except TypeError:
            return float(val)
    
    if full_output:
        y, infodict = _odeint(wrapped_func, y0, t, args,
                              Dfun, col_deriv, full_output,
                              ml, mu, rtol, atol,
                              tcrit, h0, hmax, hmin,
                              ixpr, mxstep, mxhnil, mxordn,
                              mxords, printmessg)
    else:
        y = _odeint(wrapped_func, y0, t, args,
                    Dfun, col_deriv, full_output,
                    ml, mu, rtol, atol,
                    tcrit, h0, hmax, hmin,
                    ixpr, mxstep, mxhnil, mxordn,
                    mxords, printmessg)

    # now we need to put units onto the solution units should be the
    # same as y0. We cannot put mixed units in an array, so, we return a list
    m,n = y.shape # y is an ndarray, so it has a shape
    if n > 1: # more than one equation, we need a list
        uY = [0 for yi in range(n)]
        
        for i, yi in enumerate(y0):
            if not hasattr(uY[i],'units') and hasattr(yi, 'units'):
                uY[i] = y[:,i] * yi.units
            else:
                uY[i] = y[:,i]
                
    else:
        uY = y * y0.units

    y = uY


    if full_output:
        return y, infodict
    else:
        return y

##################################################################
# test a single ODE
k = 0.23 / u.s
Ca0 = 1 * u.mol / u.L

def dCadt(Ca, t):
    return -k * Ca

tspan = np.linspace(0, 5) * u.s
sol = odeint(dCadt, Ca0, tspan)

print sol[-1]

plt.plot(tspan, sol)
plt.xlabel('Time ({0})'.format(tspan.dimensionality.latex))
plt.ylabel('$C_A$ ({0})'.format(sol.dimensionality.latex))
plt.savefig('images/ode-units-ca.png')

##################################################################
# test coupled ODEs
lbmol = 453.59237*u.mol

kprime = 0.0266 * lbmol / u.hr / u.lb
Fa0 = 1.08 * lbmol / u.hr
alpha = 0.0166 / u.lb
epsilon = -0.15

def dFdW(F, W, alpha0):
    X, y = F
    dXdW = kprime / Fa0 * (1.0 - X)/(1.0 + epsilon * X) * y
    dydW = - alpha0 * (1.0 + epsilon * X) / (2.0 * y)
    return [dXdW, dydW]

X0 = 0.0 * u.dimensionless
y0 = 1.0

# initial conditions
F0 = [X0, y0] # one without units, one with units, both are dimensionless

wspan = np.linspace(0,60) * u.lb

sol = odeint(dFdW, F0, wspan, args=(alpha,))
X, y = sol

print 'Test 2'
print X[-1]
print y[-1]

plt.figure()
plt.plot(wspan, X, wspan, y)
plt.legend(['X','$P/P_0$'])
plt.xlabel('Catalyst weight ({0})'.format(wspan.dimensionality.latex))
plt.savefig('images/ode-coupled-units-pdrpo.png')
[ 0.31663678] mol/L
Test 2
0.665569578156 dimensionless
0.263300470681

That is not too bad. This is another example of a function you would want to save in a module for reuse. There is one bad feature of the wrapped odeint function, and that is that it changes the solution for coupled ODEs from an ndarray to a list. That is necessary because you apparently cannot have mixed units in an ndarray. It is fine, however, to have a list of mixed units. This is not a huge problem, but it changes the syntax for plotting results for the wrapped odeint function compared to the unwrapped function without units.

Copyright (C) 2013 by John Kitchin. See the License for information about copying.

org-mode source

blog comments powered by Disqus