finite difference formula and its stencil

Question

I would like to prove the following finite-difference formula for functions $u$ and the flux $f$:

$$\frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t}+\frac{-f\left(u_{i+2}^{n+1}\right)+8 f\left(u_{i+1}^{n+1}\right)-8 f\left(u_{i-1}^{n+1}\right)+f\left(u_{i-2}^{n+1}\right)}{12 \Delta x}=0$$

This pair $(u,f)$ and it's discretized form as written above, should be conservative, ie. suppose mass is conserved. Then for example

$$\begin{array}{l}{\text {change in total mass in }[a, b] \text { in time interval }\left[t_{1}, t_{2}\right]=\text { net mass passing through}} \\ {\text {boundaries of }[a, b] \text { in time interval }\left[t_{1}, t_{2}\right] .}\end{array}$$

i.e.,

$$\int_{a}^{b}\left[\mathbf{u}\left(x, t_{2}\right)-\mathbf{u}\left(x, t_{1}\right)\right] d x=-\int_{t_{1}}^{t_{2}}[\mathbf{f}(b, t)-\mathbf{f}(a, t)] d t$$

Now for each cell during each time interval we have

$$\int_{x_{i-1 / 2}}^{x_{i+1 / 2}}\left[u\left(x, t^{n+1}\right)-u\left(x, t^{n}\right)\right] d x=-\int_{t^{n}}^{t^{n+1}}\left[f\left(u\left(x_{i+1 / 2}, t\right)\right)-f\left(u\left(x_{i-1 / 2}, t\right)\right)\right] d t\,\,\,\,\,(*)$$

$\begin{array}{l}{\text {where we suppose that space is divided into cells }\left[x_{i-1 / 2}, x_{i+1 / 2}\right], \text { where } x=x_{i+1 / 2} \text { is called }} \\ {\text { a cell edge. Also, suppose that time is divided into time intervals }\left[t^{n}, t^{n+1}\right], \text { where } t=t^{n}} \\ {\text { is called a time level. }}\end{array}$

Then the numerical version of this conservative form is , by using trapezoidal method to approximate the integral, as

$$\overline{u}_{i}^{n+1}=\overline{u}_{i}^{n}-\lambda\left(\hat{f}_{i+1 / 2}^{n}-\hat{f}_{i-1 / 2}^{n}\right)$$

where $$\overline{u}_{i}^{n} \approx \frac{1}{\Delta x} \int_{x_{i-1 / 2}}^{x_{i+1 / 2}} u\left(x, t^{n}\right) d x$$

$$\hat{f}_{i+1 / 2}^{n} \approx \frac{1}{\Delta t} \int_{t^{n}}^{t^{n+1}} f\left(u\left(x_{i+1 / 2}, t\right)\right) d t$$

and $$\lambda=\frac{\Delta t}{\Delta x}$$

and where $\Delta x=x_{i+1 / 2}-x_{i-1 / 2}$ and $\Delta t=t^{n+1}-t^{n}$. In the above, an overbar indicates spatial cell-integral averages, the hat means time-integtral averages

So yeah.. I am slightly convinced, given the long-windedness of the notation above, that I have complicated my task by asking the full question.. instead maybe it should have been more prudent (and sufficient) to merely ask:

about showing the relation

$$\frac{-f\left(u_{i+2}^{n+1}\right)+8 f\left(u_{i+1}^{n+1}\right)-8 f\left(u_{i-1}^{n+1}\right)+f\left(u_{i-2}^{n+1}\right)}{12 \Delta x}=0$$

might hold.. but I am not sure this would alter the problem, so I asked the full question

Answer 1

The finite difference scheme at hand can be derived by standard Taylor expansions. For the time discretization we expand around time $t_n$ and impose first order of accuracy: $$ u(t_n+\Delta t) = u(t_n) + \Delta t \dot{u}(t_n) + \mathcal{O}(\Delta t^2) \quad \Rightarrow \quad \dot{u}(t_n) = \frac{u(t_n + \Delta t) - u(t_n)}{\Delta t} + \mathcal{O}(\Delta t). $$ Here, $\dot{u}$ denotes the first derivative in time of $u$.

For the space discretization you do a similar thing: Expand $u(x_j \pm \Delta x)$ and $u(x_j \pm 2\Delta x)$ around $x_j$, combine them to solve for $u'(x_j)$ (first derivative in space) and impose $4^{\text{th}}$ order of accuracy. This will lead you to the coefficients of your spatial scheme. I'll leave the details to you.

It appears that your definition of conservation is relevant for problems posed on a finite domain whereas what you want to show is conservation for a scheme that cannot be applied without modification near domain boundaries. However, the scheme can be applied in a periodic domain, so I will answer the question with such a domain in mind.

In vector form, the stencil can be written as $$ \vec{u}^{n+1} - \vec{u}^n = \lambda D \vec{f}^n, $$ where $\vec{u}^n = (u_1^n ,u_2^n, \dots, u_{m-1}^n, u_m^n)^\top$ and $u_j^n \approx u(x_j, t_n)$ and $\vec{f}^n$ is defined similarly. Here, the matrix $D$ is given by $$ D = \frac{1}{12} \begin{pmatrix} 0 & 8 & -1 & 0 & \dots & 0 & 1 & -8 \\ -8 & 0 & 8 & -1 & 0 & \dots & 0 & 1 \\ 1 & -8 & 0 & 8 & -1 & 0 & \dots & 0 \\ & & & \ddots & \ddots & & & \\ 0 & \dots & 0 & 1 & -8 & 0 & 8 & -1 \\ -1 & 0 & \dots & 0 & 1 & -8 & 0 & 8 \\ 8 & -1 & 0 & \dots & 0 & 1 & -8 & 0 \end{pmatrix}. $$ Note that $D$ is skew-symmetric, i.e. $D^\top = -D$.

Showing that the scheme is conservative means, as you have pointed out, that some approximation of the integral of the solution over the spatial domain is conserved. A suitable approximation for the integral of the solution on a periodic domain is simply $$ \int u(x,t_n) dx \approx \Delta x \vec{1}^\top \vec{u}^n = \Delta x \sum_{j=1}^m u_j^n, $$ where the vector $\vec{1}^\top = (1,1,\dots, 1,1)$. The trapezoidal rule is not going to work here because it assumes that the domain has boundaries.

Multiplying the scheme from the left by $\Delta x \vec{1}^\top$ gives $$ \Delta x \vec{1}^\top \vec{u}^{n+1} - \Delta x \vec{1}^\top \vec{u}^n = \Delta t \vec{1}^\top D \vec{f}^n = -\Delta t (D \vec{1})^\top \vec{f}^n = 0, $$ i.e. $$ \Delta x \vec{1}^\top \vec{u}^{n+1} = \Delta x \vec{1}^\top \vec{u}^n $$ or $$ \Delta x \sum_{j=1}^m u_j^{n+1} = \Delta x \sum_{j=1}^m u_j^n. $$ In this calculation, we used the fact that $D$ is skew-symmetric in the second equality, and that $D \vec{1} = \vec{0}$ (the columns of $D$ sum to zero) in the third equality. If $u$ represents mass, then the final relation tells us that the total mass in the domain does not change from one time step to another, i.e. the scheme is conservative.