* Uncertainty in the solution of an ODE
:PROPERTIES:
:categories: ODE, uncertainty
:date: 2013/07/14 13:36:36
:updated: 2013/10/18 15:55:01
:END:
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.
#+BEGIN_SRC python
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))
#+END_SRC
#+RESULTS:
: 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.
#+BEGIN_SRC python
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)
#+END_SRC
#+RESULTS:
: 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.