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






\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)]}


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)


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.

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes:

<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>