The Kitchin Research Group

As we have seen, handling units with third party functions is fragile, and often requires additional code to wrap the function to handle the units. An alternative approach that avoids the wrapping is to rescale the equations so they are dimensionless. Then, we should be able to use all the standard external functions without modification. We obtain the final solutions by rescaling back to the answers we want.

Before doing the examples, let us consider how the quantities package handles dimensionless numbers.

import quantities as u

a = 5 * u.m
L = 10 * u.m # characteristic length

print a/L
print type(a/L)

0.5 dimensionless
<class 'quantities.quantity.Quantity'>

As you can see, the dimensionless number is scaled properly, and is listed as dimensionless. The result is still an instance of a quantities object though. That is not likely to be a problem.

Now, we consider using fsolve with dimensionless equations. Our goal is to solve \(C_A = C_{A0} \exp(-k t)\) for the time required to reach a desired \(C_A\). We let \(X = Ca / Ca0\) and \(\tau = t * k\), which leads to \(X = \exp{-\tau}\) in dimensionless terms.

import quantities as u
import numpy as np
from scipy.optimize import fsolve

CA0 = 1 * u.mol / u.L
CA = 0.01 * u.mol / u.L  # desired exit concentration
k = 1.0 / u.s

# we need new dimensionless variables
# let X = Ca / Ca0
# so, Ca = Ca0 * X

# let tau = t * k
# so t = tau / k

X = CA / CA0 # desired exit dimensionless concentration

def func(tau):
    return X - np.exp(-tau)

tauguess = 2

print func(tauguess) # confirm we have a dimensionless function

tau_sol, = fsolve(func, tauguess)
t = tau_sol / k
print t

-0.125335283237 dimensionless
4.60517018599 s

Now consider the ODE \(\frac{dCa}{dt} = -k Ca\). We let \(X = Ca/Ca0\), so \(Ca0 dX = dCa\). Let \(\tau = t * k\) which in this case is dimensionless. That means \(d\tau = k dt\). Substitution of these new variables leads to:

\(Ca0*k \frac{dX}{d\tau} = -k Ca0 X \)

or equivalently: \(\frac{dX}{d\tau} = -X \)

import quantities as u

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

# Let X = Ca/Ca0  -> Ca = Ca0 * X  dCa = dX/Ca0
# let tau = t * k -> dt = 1/k dtau


def dXdtau(X, tau):
    return -X

import numpy as np
from scipy.integrate import odeint

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

X0 = 1
X_sol = odeint(dXdtau, X0, tauspan)

print 'Ca at t = {0} = {1}'.format(tspan[-1], X_sol.flatten()[-1] * Ca0)

Ca at t = 5.0 s = 0.316636777351 mol/L

That is pretty much it. Using dimensionless quantities simplifies the need to write wrapper code, although it does increase the effort to rederive your equations (with corresponding increased opportunities to make mistakes). Using units to confirm your dimensionless derivation reduces those opportunities.

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

org-mode source