The Kitchin Research Group

Suppose we have an equation \(y = e^{a y}\) that we want to solve, where \(a\) is a constant with some uncertainty. What is the uncertainty in the solution \(y\)?

Finding a solution is not difficult. The uncertainty in the solution, however, is not easy, since we do not have an explicit function to propagate errors through. Let us examine the solution first.

import numpy as np
from scipy.optimize import fsolve

a = 0.20

def f(y):
    return y - np.exp(a * y)

sol, = fsolve(f, 1)
print sol

1.2958555091

A way to estimate the uncertainty is by Monte Carlo simulation. We solve the equation many times, using values sampled from the uncertainty distribution. Here we assume that the \(a\) parameter is normally distributed with an average of 0.2 and a std deviation of 0.02. We solve the equation 10000 times for different values of \(a\) sampled according to the normal distribution. That gives us a distribution of solutions that we can do statistical analysis of to get the average and std deviation.

import numpy as np
from scipy.optimize import fsolve
N = 10000

A = np.random.normal(0.2, 0.02, size=N)

sol = np.zeros(A.shape)

for i, a in enumerate(A):
    s, = fsolve(lambda y:y - np.exp(a * y), 1)
    sol[i] = s

ybar = np.mean(sol)
s_y = np.std(sol)

print ybar, s_y, s_y / ybar

import matplotlib.pyplot as plt
count, bins, ignored = plt.hist(sol)
plt.savefig('images/implicit-uncertainty.png')

1.29887470397 0.0465110111613 0.0358086973433

We get approximately the same answer, and you can see here the distribution of solution values is not quite normal. We compute the standard deviation anyway, and find the standard deviation is about 3.6%. It would be nice to have some analytical method to estimate this uncertainty. So far I have not figured that out.

This method could have relevance in estimating the uncertainty in the friction factor for turbulent flow (\(Re > 2100\)). In that case we have the implicit equation \(\frac{1}{\sqrt{f_F}}=4.0 \log(Re \sqrt{f_F})-0.4\). Uncertainties in the Re number would lead to uncertainties in the friction factor. Whether those uncertainties are larger than the uncertainties from the original correlation would require some investigation.

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

org-mode source