Find classical solution of transport equation with FDM

Question

We know the classical solution of transport equation is determined by one initial (boundary?) condition, for example, the solution of

$$\frac{\partial u(t,x)}{\partial t}+\frac{\partial u(t,x)}{\partial x}=0$$ $$u(0,x)=e^{-x^2}$$

is $$u(t,x)= e^{-(x-t)^2}$$

My question is, how to find this solution numerically with finite difference method (FDM)? One initial condition (i.c.) is apparently insufficient. Do we need a special difference scheme, or a special boundary condition (b.c.) that amounts to "nothing" in $x$ direction?

"But why do you need that? This problem can be easily solved analytically!" Yeah, I thought so, too, until I came across this boundary value problem of Navier-Stokes equation recently:

http://www.featflow.de/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark2_re100.html

As one can see, only one b.c. is given for $p$, so there seems to be a family of partial differential equation that "suffers" the same issue, and that's the reason I decided to look into it.

Answer 1

This equation is a pure transport equation. It defines the conservation of the quantity $u$ all over its path, i.e. $$\frac{du}{dt}=0 \qquad \forall t,x$$ and thus whatever the initial condition is, it is transported with constant velocity.

When the equation is solved in a bounded domain, the concept of boundary condition arise. It is necessary to provide what happens with the incomings at the boundary. Mathematically the operator $\partial/\partial x$ always comes with a boundary condition when is bounded. In an unbounded domain this BC is implicitly defined in the initial condition, but in a bounded one this BC must be provided. This BC is imposed in the inflow boundary, that in your case is always on the left side of your domain.

Wrt. the NS equations... the pressure does not need any BC but to be defined in one point if only Dirichlet BC are provided for the velocity. If one Do nothing BC (free stress, usually in outflow boundary) is imposed, the pressure is completely defined implicitly in the normal stress: $$\vec{n}^T\sigma\vec{n}|_\Gamma=-p+2\vec{n}^T\frac{\partial \vec{v}}{\partial n}=0$$

EDIT

Let us rewrite your problem to be able to compare analytical vs numerical results. Let us define a bounded region in which the equation

$$\frac{\partial u}{\partial t}+v\frac{\partial u}{\partial x}=0$$ is valid, namely, for $t_0>0$ and $x\in[x_1,x_2]$ with an initial condition of :$u(x,t_0)=u_0(x)$ Therefore, due to the fact of the boundedness of the domain an inflow BC is required. Because the transport velocity $v$ is constant and positive the BC will be for the extremum in the coordinate $x_1$. Let this BC be $u(x_1,t)=g(t)$.

For the sake of simplicity let us take $v=1$.

Now the analytical solution is not what you said, but $u_0(x)$ moving right and a new curve drawn by the function $g(t)$.

With ANY discretisation, say FDM, FEM, FVM... you will have a discrete system as follows: $$\frac{du}{dt}+Du+BC=0$$ Where $u=[u_2,...,u_N]^T$, and $BC=[g(t)/\Delta x,0,...,0]$ and $D(i,i)=1/\Delta x,D(i,i-1)=-1/\Delta x$ for $i\geq 2$ and $D(1,1)=1$. It is clear that $u_1=g(t)$.

For the case of NSE, one has the following system: $$div(\vec{v})=0$$ $$\frac{\partial \vec{v}}{\partial t}+(\vec{v}^T\vec{grad})\vec{v}=div(\widetilde{\sigma})$$ with $\widetilde{\sigma}=-p\mathbb{I} +\tfrac{1}{2}[\vec{grad}\,\vec{v}+(\vec{grad}\,\vec{v})^T] $

This can be spatially discretised, following any method, as follows: $$D\vec{v}=0$$ $$ \frac{d \vec{v}}{dt}+\vec{H}(\vec{v})=L\vec{V}-Gp+\vec{BC}$$

In order to have a stable discretisation the rank of the divergence matrix $D$ must be less than the number of velocity unknowns. Otherwise $\vec{v}$ would be uniquely determined as the kernel of $D$. This induces the usage of different meshes for pressure and velocity, i.e. you cannot use whatever mesh you want to solve the problem due to stability issues.

The operators $G$ and $L$ are the gradient and Laplacian respectively. The vector $\vec{H}$ corresponds to the nonlinear terms in NSE. The system, once discretised in time may be arranged as follows: $$\left[\matrix{A &G\\ D&0}\right]\left[\matrix{\vec{v}^{n+1}\\ p^{n+1}}\right]=\left[\matrix{\vec{R}(\vec{v}^n,\vec{v}^{n-1})\\ 0}\right]+\vec{BC}$$ One example of matrix $A$ coincides with: $$A=\left[\mathbb{I}-c\frac{\Delta t}{Re}L\right]/\Delta t$$ Where c is a constant and depends on the selected temporal scheme.

The term $\vec{BC}$ only has BC for the velocity. This can be easily seen applying FEM. For FDM pressure is usually calculated at interior nodes. Keep in mind that the pressure is a Lagrange multiplier that forces the incompressibility condition, and plays a mathematical rather a physical role. Thats the reason of its lack of BC.

Particularisation for finite difference scheme

As I previously commented, the mesh must be different for pressure and for velocity (Inf-sup condition). You can use an square element in which there is only a pressure node in the middle, 2 nodes for vertical component of velocity in the upper and lower part, and finally 2 more for the horizontal component of the velocity in the right and left part of the element. This mesh is known as a staggered grid.

This has a clear advatage: the face fluxes are calculated exactly in order to assure the incompressibility condition.

As you can see, due to the position of the pressure nodes, it is not required any boundary condition.

Answer 2

Let me try to answer your question concerning the transport equation from a different perspective. Consider the following initial-value problem. We search for a function $v(x,t)$ such that \begin{align*} v_t + v_x &= 0 \\ v(x, 0) &= f(x) \end{align*} for $t \ge 0$ and $x \in \mathbb{R}$. The solution of this PDE can be approximated using the upwind scheme, which can be derived as follows.

Let $\Delta t > 0$ and $\Delta x > 0$ be two small real numbers. We can use that $$ v_t(x,t) \approx \frac{v(x, t+\Delta t) - v(x, t)}{\Delta t} \quad \text{and} \quad v_x(x,t) \approx \frac{v(x, t) - v(x - \Delta x, t)}{\Delta x} $$ to find an approximation to the solution. Replacing the partial derivatives in the original equation, we find that $$ \frac{v(x, t+\Delta t) - v(x, t)}{\Delta t} + \frac{v(x, t) - v(x - \Delta x, t)}{\Delta x} = 0 $$ approximates the PDE. Solving this equation for $v(x, t+\Delta t)$ gives $$ v(x, t+\Delta t) = v(x, t) + \frac{\Delta t}{\Delta x}(v(x, t) - v(x - \Delta x, t)) \,. $$

Now, to turn this approximation into an algorithm, consider the grid $$ \mathcal{G} := \{ (k\, \Delta x, n\, \Delta t) : k \in \mathbb{Z}, n \in \mathbb{N}_0 \} \,, $$ and denote $u_k^n$ as the approximation to the value of $v(k\, \Delta x, n\, \Delta t)$. Then, \begin{align*} u_{k}^0 &= f(k\,\Delta x) \\ u_{k}^{n+1} &= u_k^n + \frac{\Delta t}{\Delta x}(u_k^n - u_{k-1}^n) \end{align*} is the upwind scheme for solving the transport equation. Which is a convergent scheme in case $\Delta t$ is not too large in comparison to $\Delta x$.

The problem is, that the grid $\mathcal{G}$ has infinitely many points and thus is not suitable for a computer implementation. However, if we take a closer look, we will realize that we actually only need a finite subset of $\mathcal{G}$.

Assume we would like to compute the value $u_k^n$ using the upwind scheme. We see that to compute $u_k^n$, we need to compute the values of $u_{k-1}^{n-1}$ and $u_k^{n-1}$, first. Using the same reasoning, we see that to compute these two values, we need the values of $u_{k-2}^{n-2}$, $u_{k-1}^{n-2}$, and $u_{k}^{n-2}$. We can repeat this steps until we finally reach $n = 0$, where we can use the values given by the function $f$. Hence, to compute the value of $u_k^n$, we need the values given in the following diagramm: $$ \newcommand{\iddots}{\kern3mu\raise1mu{.}\kern3mu\raise6mu{.}\kern3mu\raise12mu{.}} \begin{matrix} &&&& u_k^n \\ &&& u_{k-1}^{n-1} & u_{k}^{n-1} \\ && u_{k-2}^{n-2} & u_{k-1}^{n-2} & u_{k}^{n-2} \\ & \iddots & \vdots & \vdots & \vdots \\ u_{k-n}^0 & \dots & u_{k-2}^{0} & u_{k-1}^{0} & u_{k}^{0} \end{matrix} $$ The triangle shown above is the so called numerical domain of dependence of the point $u_k^n$.

What you can see is that to compute the point $u_k^n$, you just need to compute the values in this triangular subset of the grid. Hence, no additional boundary conditions are needed. In case you want to compute more than one point, your domain of dependence will be trapezoidal, as can be seen easily. Note, that for other PDEs these shapes will usually be different.

The reason for this behavior is that the transport equation is a hyperbolic PDE, and thus the information of the initial values are only transported with finite speed, in contrast to, e.g., the heat equation where the speed of transport of information is infinite.

Note, that for the scheme to actually converge to the true solution, the domain of dependence of the numerical scheme has to include the domain of dependence of the PDE, otherwise the scheme will diverge. This is the so called Courant–Friedrichs–Lewy (CFL) condition, which explains the restriction on the step size $\Delta t$, mentioned earlier.