Inverse advection-diffusion problem, solving for a drift coefficient with experimental data?

Question

I am investigating a physical process where I believe the 1-D advection-diffusion equation: \begin{equation} \frac{\partial u}{\partial t} = -\frac{\partial}{\partial x}[\mu(x,t) u(x,t)] + \frac{\partial}{\partial x}[D(x,t)\frac{\partial}{\partial x}u(x,t)] \end{equation}

may be applicable. I have a set of experimental data of the density, $u_{exp}(x,t)$, across $x$ for each time $t$. To begin my investigation of this I want to assume that $D(x,t)$ is a constant, $D(x,t)$=$D_0$, and $\mu(x,t) = \mu_0 + \tilde{\mu}(x,t) $ such that I can approach the problem numerically as:

\begin{equation} \frac{\partial u}{\partial t} = -\frac{\partial}{\partial x}[\mu(x,t)] u(x,t) -\mu(x,t)\frac{\partial}{\partial x}[u(x,t)] + D_0\frac{\partial^2} {\partial x^2}[u(x,t)] \end{equation}

\begin{equation} \frac{\partial u}{\partial t} = -\frac{\partial}{\partial x}[\mu_0 + \tilde{\mu}(x,t)] u(x,t) -[\mu_0 + \tilde{\mu}(x,t)]\frac{\partial}{\partial x}[u(x,t)] + D_0\frac{\partial^2} {\partial x^2}[u(x,t)] \end{equation}

Rearranging and expanding, I can pick out the terms such that this resembles a typical advection-diffusion equation (with constant diffusion and drift coefficients) with an additional term like so:

\begin{equation} \frac{\partial u}{\partial t} = \bigg[ -\mu_0\frac{\partial}{\partial x}[u(x,t)] + D_0\frac{\partial^2}{\partial x^2}[u(x,t)] \bigg]_{linear-adv.-diff. \;equation} \; -\bigg[ \frac{\partial}{\partial x}[\tilde{\mu}(x,t)] u(x,t) +\tilde{\mu}(x,t)\frac{\partial}{\partial x}[u(x,t)] \bigg]_{additional\;term} \end{equation}

Now I'll define:

\begin{equation} f_0(x,t) = -\mu_0\frac{\partial}{\partial x}[u(x,t)] + D_0\frac{\partial^2}{\partial x^2}[u(x,t)] \end{equation}

To rewrite as:

\begin{equation} \frac{\partial u}{\partial t} = f_0(x,t) -\frac{\partial}{\partial x}[\tilde{\mu}(x,t)] u(x,t) -\tilde{\mu}(x,t)\frac{\partial}{\partial x}[u(x,t)] \end{equation}

\begin{equation} \frac{\partial u}{\partial t} = f_0(x,t) -\tilde{\mu}_x(x,t) u(x,t) -\tilde{\mu}(x,t)u_x(x,t) \end{equation}

Here I now want to begin supplementing my equation with my experimental data, and predictions from a numerical simulation of the linear-adv.-diff. equation, to get this show on the road. I substitute in a numerical approximation for $f_0(x,t')$, $f_{num}(x,t')$, at time $t'$ coming from a realization with the initial profile $u_{exp}(x,0)$, and coefficients $u_0$ and $D_0$. I thus write the error for my simulation, using this numerical prediction relative to the experimental $\mu_{exp}$ as:

\begin{equation} \bigg[ \frac{\partial u_{exp}(x,t)}{\partial t} - f_{num}(x,t) = -\tilde{\mu}_x(x,t) u_{exp}(x,t) -\tilde{\mu}(x,t)u_{x,exp}(x,t) \bigg]\bigg|_{t=t'} \end{equation}

\begin{equation} \bigg[ - \epsilon(x,t) = \tilde{\mu}_x(x,t) u_{exp}(x,t) +\tilde{\mu}(x,t)u_{x,exp}(x,t) \bigg]\bigg|_{t=t'} \end{equation}

where $\epsilon$ denotes the error in the prediction across $x$ at time $t'$. Next we recast this into the form of a familiar first-order-differential equation at time $t'$:

\begin{equation} \bigg[ \tilde{\mu}_x(x,t) +\tilde{\mu}(x,t) \frac{u_{x,exp}(x,t)}{u_{exp}(x,t)} = - \frac{\epsilon(x,t)}{u_{exp}(x,t)} \bigg]\bigg|_{t=t'} \end{equation}

Thus, where $p(x,t) = \frac{u_{x,exp}(x,t)}{u_{exp}(x,t)}$, and $q(x,t) = - \frac{\epsilon(x,t)}{u_{exp}(x,t)}$:

\begin{equation} \bigg[ \tilde{\mu}_x(x,t) + p(x,t) \tilde{\mu}(x,t) = q(x,t) \bigg]\bigg|_{t=t'} \end{equation}

This can then be solved by the use of an integrating factor, as detailed: http://tutorial.math.lamar.edu/Classes/DE/Linear.aspx

I suppose my question is: is there a more straight-forward method to solving for the expression $\mu(x,t)$ than first assuming a linear portion of the equation and then rederiving the $\mu(x,t)$ profile as I have done here? Also, is what I've done here coherent? I mean, is it permissible even to try to solve for $\mu(x,t)$ via the error in the linear-adv.-diff. approximation? A motivating image for this is shown below:

(Top) Experimental data. (Bottom) Linear-Adv.-Diff equation. Time is colored, with red being $t=0$ and purple being $t=t_{final}$

Answer 1

There is another approach called the adjoint method, which is commonly used in inverse problems for PDE and which is quite easy to generalize to other problems. This is going to be long.

Your observations give you a field $u_{\text{exp}}$; you'd like to find a value of $\mu$ for which

$\frac{1}{2}\iint(u - u_\text{exp})^2dx\hspace{2pt}dt$

is a minimum, where $u$ satisfies the advection-diffusion equation with $\mu$ as the drift coefficient. That's a constraint, which we can enforce by the introduction of a Lagrange multiplier $\lambda$:

$\iint\lambda\left\{\frac{\partial u}{\partial t} + \frac{\partial}{\partial x}(\mu u) - \frac{\partial}{\partial x}\left(D\frac{\partial u}{\partial x}\right)\right\}dx\hspace{2pt}dt$

So, altogether now, we seek 3 fields $u$, $\mu$, $\lambda$ which are an extremum of the functional

$J[u, \mu, \lambda] = \frac{1}{2}\iint(u - u_\text{exp})^2dx\hspace{2pt}dt + \iint\lambda\left\{\frac{\partial u}{\partial t} + \frac{\partial}{\partial x}(\mu u) - \frac{\partial}{\partial x}\left(D\frac{\partial u}{\partial x}\right)\right\}dx\hspace{2pt}dt$.

That's a constrained optimization problem; the first term is what we're optimizing, and the second is the constraint.

Using the usual technique of the calculus of variations, you can compute the gradient $\nabla J$ of the functional $J$ to a small perturbation $\delta u$, $\delta \mu$, $\delta\lambda$, which gives you a set of PDEs to solve:

$\frac{\partial u}{\partial t} + \frac{\partial}{\partial x}(\mu u) - \frac{\partial}{\partial x}\left(D\frac{\partial u}{\partial x}\right) = 0$

$\frac{\partial\lambda}{\partial t} + \mu\frac{\partial\lambda}{\partial x} + \frac{\partial}{\partial x}\left(D\frac{\partial\lambda}{\partial x}\right) = u - u_{\text{exp}}$

The first one is (unsurprisingly) your original PDE. The second PDE looks odd; the diffusion coefficient is negative. However, no one said you had to solve it from $t = 0$ to $t = T$; you can go backwards and then everything's fine. This second system is called the adjoint equation. Note that if you had a perfect guess for $\mu$ which made $u$ exactly equal to the experimental data, the right-hand side of the adjoint equation would be 0 and $\lambda$ would also be zero.

Having solved both of these PDE, the change in $J$ with respect to a small perturbation $\delta\mu$ is

$\nabla_\mu J\delta\mu = -\iint u\frac{\partial\lambda}{\partial x}\delta\mu\hspace{3pt} dx\hspace{2pt}dt$.

But, we know that $J$ decreases fastest in the direction opposite its gradient, so if you wanted to improve $J$, you could set

$\mu' = \mu + \alpha u\frac{\partial\lambda}{\partial x}$

and search for a value of the scalar $\alpha$ that works best. Or, you could use something more sophisticated like a quasi-Newton method, but that's just a detail.

In summary:

1. Pick some guess for $\mu$
2. Solve the forward PDE for $u$
3. Solve the adjoint PDE for $\lambda$
4. Compute the gradient $\nabla J$
5. Update $\mu$ and repeat if the error isn't good enough

I can offer more detail if necessary, there's quite a bit left out. I don't know of a good introductory reference for the subject that isn't too general to be useful.

Generally, you want to avoid numerically differentiating experimental data. If it has some random red noise errors on top of the actual signal and you differentiate those, you get a bunch of nonsense out; you could apply a low-pass filter to it first to avoid that. In the adjoint method, you can take advantage of all the smoothing properties of parabolic PDE to do that for you.