## Uncertainty in the solution of an ODE

Posted July 14, 2013 at 01:36 PM | categories: ode, uncertainty | tags: | View Comments

Updated October 18, 2013 at 03:55 PM

Our objective in this post is to examine the effects of uncertainty in parameters that define an ODE on the integrated solution of the ODE. My favorite method for numerical uncertainty analysis is Monte Carlo simulation because it is easy to code and usually easy to understand. We take that approach first.

The problem to solve is to estimate the conversion in a constant volume batch reactor with a second order reaction \(A \rightarrow B\), and the rate law: \(-r_A = k C_A^2\), after one hour of reaction. There is 5% uncertainty in the rate constant \(k=0.001\) and in the initial concentration \(C_{A0}=1\).

The relevant differential equation is:

\(\frac{dX}{dt} = -r_A /C_{A0}\).

We have to assume that 5% uncertainty refers to a normal distribution of error that has a standard deviation of 5% of the mean value.

from scipy.integrate import odeint import numpy as np N = 1000 K = np.random.normal(0.001, 0.05*0.001, N) CA0 = np.random.normal(1, 0.05*1, N) X = [] # to store answer in for k, Ca0 in zip(K, CA0): # define ODE def ode(X, t): ra = -k * (Ca0 * (1 - X))**2 return -ra / Ca0 X0 = 0 tspan = np.linspace(0,3600) sol = odeint(ode, X0, tspan) X += [sol[-1][0]] s = 'Final conversion at one hour is {0:1.3f} +- {1:1.3f} (1 sigma)' print s.format(np.average(X), np.std(X))

Final conversion at one hour is 0.782 +- 0.013 (1 sigma)

See, it is not too difficulty to write. It is however, a little on the expensive side to run, since we typically need 1e3-1e6 samples to get the statistics reasonable. Let us try the uncertainties package too. For this we have to wrap a function that takes uncertainties and returns a single float number.

from scipy.integrate import odeint import numpy as np import uncertainties as u k = u.ufloat(0.001, 0.05*0.001) Ca0 = u.ufloat(1.0, 0.05) @u.wrap def func(k, Ca0): # define the ODE def ode(X, t): ra = -k * (Ca0 * (1 - X))**2 return -ra / Ca0 X0 = 0 # initial condition tspan = np.linspace(0, 3600) # integrate it sol = odeint(ode, X0, tspan) return sol[-1][0] result = func(k, Ca0) s = 'Final conversion at one hour is {0}(1 sigma)' print s.format(result)

Final conversion at one hour is 0.783+/-0.012(1 sigma)

This is about the same amount of code as the Monte Carlo approach, but it runs much faster, and gets approximately the same results. You have to remember the wrapping technique, since the uncertainties package does not run natively with the odeint function.

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