The Kitchin Research Group

In the previous section we examined an analytical approach to error propagation, and a simulation based approach. There is another approach to error propagation, using the uncertainties module (https://pypi.python.org/pypi/uncertainties/). You have to install this package, e.g. pip install uncertainties. After that, the module provides new classes of numbers and functions that incorporate uncertainty and propagate the uncertainty through the functions. In the examples that follow, we repeat the calculations from the previous section using the uncertainties module.

Addition and subtraction

import uncertainties as u

A = u.ufloat((2.5, 0.4))
B = u.ufloat((4.1, 0.3))
print A + B
print A - B

>>> >>> >>> 6.6+/-0.5
-1.6+/-0.5

Multiplication and division

F = u.ufloat((25, 1))
x = u.ufloat((6.4, 0.4))

t = F * x
print t

d = F / x
print d

>>> >>> >>> 160.0+/-11.8726576637
>>> >>> 3.90625+/-0.289859806243

Exponentiation

t = u.ufloat((2.03, 0.0203))
print t**5

from uncertainties.umath import sqrt
A = u.ufloat((16.07, 0.06))
print sqrt(A)
# print np.sqrt(A) # this does not work

from uncertainties import unumpy as unp
print unp.sqrt(A)

34.4730881243+/-1.72365440621
>>> >>> >>> >>> 4.00874045057+/-0.00748364738749
... >>> >>> 4.00874045057+/-0.00748364738749

Note in the last example, we had to either import a function from uncertainties.umath or import a special version of numpy that handles uncertainty. This may be a limitation of teh uncertainties package as not all functions in arbitrary modules can be covered. Note, however, that you can wrap a function to make it handle uncertainty like this.

import numpy as np

wrapped_sqrt = u.wrap(np.sqrt)
print wrapped_sqrt(A)

>>> >>> 4.00874045057+/-0.00748364738774

Propagation of errors in an integral

import numpy as np
import uncertainties as u

x = np.array([u.ufloat((1, 0.01)), 
              u.ufloat((2, 0.1)),
              u.ufloat((3, 0.1))])

y = 2 * x

print np.trapz(x, y)

>>> >>> ... ... >>> >>> >>> >>> 8.0+/-0.600333240792

Chain rule in error propagation

v0 = u.ufloat((1.2, 0.02))
a = u.ufloat((3.0, 0.3))
t = u.ufloat((12.0, 0.12))

v = v0 + a * t
print v

>>> >>> >>> >>> 37.2+/-3.61801050303

A real example? This is what I would setup for a real working example. We try to compute the exit concentration from a CSTR. The idea is to wrap the “external” fsolve function using the uncertainties.wrap function, which handles the units. Unfortunately, it does not work, and it is not clear why. But see the following discussion for a fix.

from scipy.optimize import fsolve

Fa0 = u.ufloat((5.0, 0.05))
v0 = u.ufloat((10., 0.1))

V = u.ufloat((66000.0, 100))  # reactor volume L^3
k = u.ufloat((3.0, 0.2))      # rate constant L/mol/h

def func(Ca):
    "Mole balance for a CSTR. Solve this equation for func(Ca)=0"
    Fa = v0 * Ca     # exit molar flow of A
    ra = -k * Ca**2  # rate of reaction of A L/mol/h
    return Fa0 - Fa + V * ra

# CA guess that that 90 % is reacted away
CA_guess = 0.1 * Fa0 / v0

wrapped_fsolve = u.wrap(fsolve)
CA_sol = wrapped_fsolve(func, CA_guess)

print 'The exit concentration is {0} mol/L'.format(CA_sol)

>>> >>> >>> >>> >>> >>> >>> ... ... ... ... ... >>> ... >>> >>> >>> <function fsolve at 0x148f25f0>
>>> >>> The exit concentration is NotImplemented mol/L

I got a note from the author of the uncertainties package explaining the cryptic error above, and a solution for it. The error arises because fsolve does not know how to deal with uncertainties. The idea is to create a function that returns a float, when everything is given as a float. Then, we wrap the fsolve call, and finally wrap the wrapped fsolve call!

• Step 1. Write the function to solve with arguments for all unitted quantities. This function may be called with uncertainties, or with floats.
• Step 2. Wrap the call to fsolve in a function that takes all the parameters as arguments, and that returns the solution.
• Step 3. Use uncertainties.wrap to wrap the function in Step 2 to get the answer with uncertainties.

Here is the code that does work:

import uncertainties as u
from scipy.optimize import fsolve

Fa0 = u.ufloat((5.0, 0.05))
v0 = u.ufloat((10., 0.1))

V = u.ufloat((66000.0, 100.0))  # reactor volume L^3
k = u.ufloat((3.0, 0.2))      # rate constant L/mol/h

# Step 1
def func(Ca, v0, k, Fa0, V):
    "Mole balance for a CSTR. Solve this equation for func(Ca)=0"
    Fa = v0 * Ca     # exit molar flow of A
    ra = -k * Ca**2  # rate of reaction of A L/mol/h
    return Fa0 - Fa + V * ra

# Step 2
def Ca_solve(v0, k, Fa0, V): 
    'wrap fsolve to pass parameters as float or units'
    # this line is a little fragile. You must put [0] at the end or
    # you get the NotImplemented result
    sol = fsolve(func, 0.1 * Fa0 / v0, args=(v0, k, Fa0, V))[0]
    return sol

# Step 3
print u.wrap(Ca_solve)(v0, k, Fa0, V)

0.005+/-0.000167764327667

It would take some practice to get used to this, but the payoff is that you have an “automatic” error propagation method.

Being ever the skeptic, let us compare the result above to the Monte Carlo approach to error estimation below.

import numpy as np
from scipy.optimize import fsolve

N = 10000
Fa0 = np.random.normal(5, 0.05, (1, N))
v0 = np.random.normal(10.0, 0.1, (1, N))
V =  np.random.normal(66000, 100, (1,N))
k = np.random.normal(3.0, 0.2, (1, N))

SOL = np.zeros((1, N))

for i in range(N):
    def func(Ca):
        return Fa0[0,i] - v0[0,i] * Ca + V[0,i] * (-k[0,i] * Ca**2)
    SOL[0,i] = fsolve(func, 0.1 * Fa0[0,i] / v0[0,i])[0]

print 'Ca(exit) = {0}+/-{1}'.format(np.mean(SOL), np.std(SOL))

Ca(exit) = 0.00500829453185+/-0.000169103578901

I am pretty content those are the same!

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

org-mode source