Quantifying Uncertainty with the PDF Method

In the last century, scientists and mathematicians have discovered many ways to combine what we know with what we don’t. We no longer believe that the physical world can be explained and predicted with purely deterministic equations, as Newton and Laplace did. It turns out that we can’t get rid of randomness. Most physical phenomena cannot be measured, modeled and/or predicted exactly. They can only be studied probabilistically.

In some cases, probability is a fundamental part of the physics like in quantum mechanics. In other cases, the physical state we are observing is very complex that it only makes sense study it statistically. For example, it makes more sense to model the concentration of contaminants in an aquifer rather than track every single molecule in the fluid. Not only is it impossible to track every particle, but it’s also useless. Sometimes, the motion is chaotic even when the system looks simple, like in the case of a double pendulum. In fact, there is a transition from seemingly simple and predictable motion to completely random processes. This is the subject matter of nonlinear dynamics and complex systems.

In engineering applications, uncertainty is a practical problem. In most applications, we only have access to a very limited number of measurements. In addition, the measurements themselves have a finite accuracy. Some problems are too large for any accurate measurements to be even possible. How would you predict the weather tomorrow based on what you see today? How would you predict the propagation of pressure waves in the Earth’s crust around an earthquake?

There is no way to know everything about a real world physical system. But it’s also not enough to just study the “average” behavior. For example, robots are much better at navigating their environment when they’re uncertain about their sensors measurements. So what’s the solution?

Many engineering methods are used to account for this kind of uncertainty. In this articles I will talk mainly about the PDF method. This method was popularized in modeling turbulent flows and combustion. But it can be applied to any set of differential equations where the parameters, initial and/or boundary conditions are uncertain.

Problem definition

We start from a physical problem with a known set of differential equations for a space-time dependent variable $u(x, t)$ given by

$\displaystyle{ \frac{\partial u}{\partial t} + \mathcal{L}_x^\mu u= 0 }$

where $\mathcal L_x^\mu$ is a differential operator with respect to $x$ (e.g. $\partial ^2/\partial x^2$ ), and $\mu$ is a physical parameter (e.g. viscosity, porosity).

To solve the problem, we need initial conditions (IC), $u(x, 0) = u_0(x)$ , and/or boundary conditions (BC), $u(x_b, t) = u_b(x_b, t)$ , depending on the properties of the differential equation. The uncertainty can be introduced through the ICs, BCs and/or parameters, by defining a corresponding probability density function (PDF). The PDF of an uncertainty initial condition is given by $f_{u_0}(U; x)$ , the PDF of an uncertainty boundary condition by $f_{u_b}(U; x_b, t)$ , and the PDF of a random parameter by $f_{\mu}(M)$ (assuming the parameter is not space-time dependent, which can be relaxed).

The uncertainty in the IC, BC and/or parameters is “propagated” to the function $u(x, t)$ we are interested to solve; making it also a random field. In this setup, instead of solving for $u(x, t)$ , we are interested in finding $f_u(U; x, t)$ , or, more generally, the joint PDF $f_{u \mu}(U, M; x, t)$ .

We can compute a good approximation of $f_{u \mu}(U, M; x, t)$ numerically using Monte Carlo methods, but at a very high computational cost. The purpose of the PDF method is to systematically transform the original set of differential equations in $u(x, t)$ into a single one for $f_{u \mu}(U, M; x, t)$ . In short, we’re trying to make this transformation

Deterministic

Probabilistic

PDE/ODE:

$\displaystyle{\frac{\partial u(x, t)}{\partial t} = \mathcal{L}_x^\mu[u(x,t); \mu]}$

$\displaystyle{\frac{\partial f_{u\mu}(U, M; x, t)}{\partial t} = \mathcal{M}_{xU}^\mu[f_{u\mu}(U, M; x, t)]}$

B.C:

$u(x_b, t) = u_b(x_b, t)$

$f_{u\mu}(U, M; x_b, t) = \hat{g}_{u\mu}(U, M; x_b, t)$

I.C:

$u(x, 0) = u_0(x)$

$f_{u\mu}(U, M; x, 0) = \hat{h}_{u\mu}(U, M, x)$

Note that $\hat{g}_{u\mu}(U, M; x_b, t)$ and $\hat{h}_{u\mu}(U, M, x)$ are known boundary and initial conditions of the probabilistic problem, where the mean is usually their deterministic counterpart. Our purpose is to find the new differential operator $\mathcal{M}_{xU}^\mu$ (now also a function of $U$ ), and subsequently solve the equation for the joint PDF. This is sometimes possible analytically, but can more generally be done numerically.

A simple example

The best way to introduce the method is through a simple example. Let’s start with a simple ODE for a time-dependent variable $u(t)$ ,

(1) $\begin{equation*} \frac{du}{dt} = ku \end{equation*}$

with parameter $k$ and initial condition $u(0) = u_0$ . Let’s assume that only the initial condition is uncertain, with PDF $f_{u_0}(U)$ .

The PDF method starts by defining a generalized function we call the raw PDF, given by

$\displaystyle{ \pi_u = \delta( u(t) - U ) }$

where $\delta(\cdot)$ is a delta function and $U$ the sample space variable; that is, $U$ is a sample from the random function $u(t)$ . The nice thing about this strange mathematical object is that its ensemble average is nothing but the PDF. That is

$\displaystyle{ \langle \pi_u \rangle = \int_{-\infty}^{+\infty} f_u(u', t) \delta(u' - U) du' = f_u(U, t) }$

The goal is to use this property to transform our ODE of $u$ into a PDE of $f_u(U, t)$ . For the simple example at hand, here are the steps

Multiply both sides of the equation by $\partial \pi_u / \partial u$

$\displaystyle{ \frac{\partial \pi_u}{\partial u} \frac{du}{dt} = \frac{\partial \pi_u}{\partial u} k u }$

Note the chain rule on the left hand side, and use the property $\partial \pi_u/\partial u = -\partial\pi_u/\partial U$ (explained later), and bring $u$ into the derivative because it’s independent of $U$

$\displaystyle{ \frac{\partial \pi_u}{\partial t} = - k \frac{\partial }{\partial U} ( \pi_u u ) }$

Using the variable switching property $\delta(u(t) - U) u = \delta(u(t) - U) U$ characteristic of delta functions. Then take an ensemble average

$\displaystyle{ \frac{\partial f_u}{\partial t} + k U \frac{\partial f_u }{\partial U} + kf_u = 0 }$

This is a linear hyperbolic equation that can be solved analytically using the method of characteristics. If we assume no boundaries in this problem (thus no boundary conditions), and recall that the initial condition is $f_u(U, 0) = f_{u_0}(U)$ , the solution is given by

$\displaystyle{ f_u(U, t) = f_{u_0}(U e^{-kt}) e^{-kt} }$

Remember that the solution of the deterministic problem is $u_0 e^{-kt}$ . For a $k>0$ the solution will die out and for $k<0$ it will explode. What should happen to the distribution? Here’s are plots of the solution, with different values of $k$ , with $u_0 =1$ , and with a Gaussian initial distribution.

Notice that the have the same distribution at $t=0$ (with mean and mode 1.0). For $k<0$ , all initial conditions die out towards $U=0$ ; thus increasing the density around $U=0$ . The opposite effect happens for positive $k$ .

Does it always work?

In this simple example, we are twice lucky. First, we can derive a PDF equation, and second we can solve it analytically. In many cases, when the PDF equation is complicated, the second step is done numerically. This is usually not a big deal because PDF equations are linear (although high dimensional). But more importantly, there is a large class of differential equations (like a diffusion equation) for which the first step doesn’t work. Only a few terms in the equations can be turned into functions of the joint PDF; while the rest needs to be approximated. There are many ways to make those approximations, and they usually depend on the problem at hand. It would be great if we could have a general method for making those approximations. Can machine learning help?

In a recent paper, we propose a data-driven method to discover the PDF equation from a collection of single realizations (Monte Carlo simulations). This is a promising direction for data-driven discovery of differential equations in general. In the future, I will post more on how machine learning is changing the way we do physics.

Problem definition

A simple example

Does it always work?

Leave a Reply Cancel Reply