At the Section 4.2 of this paper (which is very well known in the computational fluid dynamic community), the author claims that it is enough, for the compact finite difference formulation in eq. 4.2.3 (as applied to problem 4.2.1) to be conservative, that the matrix $B$ (note, the derivation matrix is $A^{-1}B$) have all columns but first and last which sum to zero, but gives no proof or reference for this.

Well I can't get this obvious proof.

I've searched on the web for other papers (citing or not this one) and found some results (see Section 2.1 and Appendix A) and proofs for cases which are much simpler.


The author first considers the model equation

$$\partial_t f + \partial_x F = 0, \quad x \in [a,b] \subseteq \mathbb R, t \in \mathbb R^+$$

for the quantity $f(x,t)$, which has flux $F(x,t) = \hat{F}(f(x,t))$. By integrating over $[a,b]$, it is shown that $f$ is conserved, in that the total $f$ in the domain $[a,b]$ changes in time only due to the flux $F$ at the boundaries.

$$\frac{d}{dt}\int^b_a f(x,t) dx = - (F(b,t) - F(a,t))\tag{1}\label{integral}$$

(All this is clear to me so far.)

At this point the author seeks a discrete analog for Eq. \eqref{integral}.

What I would do is the following.

To keep it simple and get the basic idea, I would consider a linear version of the model equation, i.e. the linear advection equation with $F = f$, $$\partial_t f + \partial_x f = 0$$ then I would discretize it in space,

$$ \frac{\text d}{\text dt} \textbf f + \textbf D \textbf f = 0 \tag{2}\label{discretized}$$

where $\textbf f$ contains nodal values of $f$ and $\textbf D$ is a derivation matrix. A reasonable discrete analog for the integral $\int^b_a f(x,t) dx$ is, in my opinion, the sum of all nodal values of $\textbf f$, i.e. $\textbf 1^{\text T} \textbf f$, being $\textbf 1$ the unitary column, so I'd multiply Eq. \eqref{discretized} by $\textbf 1^{\text T}$ (and take it in the time derivative), obtaining

$$\frac{\text d}{\text dt} \textbf 1^{\text T} \textbf f = - \textbf 1^{\text T}\textbf D \textbf f \tag{3}\label{guilty}$$

At this point I'd conclude that, for the LHS of Eq. \eqref{guilty} to depend only on first and last elements of $\textbf f$ in the way as Eq. \eqref{integral}, the row vector $\textbf 1^{\text T}\textbf D$ should be $(-1,0,0,\dots,0,0,1)$, i.e. all columns of $\mathbf D$ should sum up to zero, except the first and last (which should sum up to $-1$ and $1$ respectively).

Turning point

The paper in question is about compact finite difference schemes, in which case the matrix $\mathbf D$ can be written as $\mathbf A^{-1} \mathbf B$.

The author claims that

In order to satisfy the global conservation constraint it is sufficient to require that the columns $2$ through $N-1$ of the matrix $\mathbf B$ sum exactly to zero. This ensures that only the boundary nodes contribute to the boundary fluxes.

which is clearly not consistent with what I wrote just after Eq. \eqref{guilty}.

So either the author is wrong (unlikely, since the paper has been cited more than five thousands times since 1992), or I misunderstood what he means by global conservation constraint as applied to the space-discrete equation or, in other words, multiplying Eq. \eqref{discretized} by $\mathbf 1^\text{T}$ is not the intended discrete analog of integrating $f$ over the domain.


The equation: $$\partial_t f + \partial_xF = 0$$ with $f_{x=0} = f_0(t)$ (which is known) and appropiate initial conditions can be discretised as follows: $$\partial_t\vec{f}+D\vec{F}=\vec{b}\tag{*}$$ Where $\vec{f} = (f_1,f_2,...,f_{N+1})^T$ and $D$ is the discrete operator which replaces $\partial_x$, and is equal to: $$D = \frac{1}{h}\left[\begin{matrix}1 & & & &\\ -1 &1&&&\\ & -1&1&&\\ & &-1&1&&\\ & & & \ddots& \ddots \end{matrix}\right] $$ and $\vec{b} = (F(f_0(t))/h,0,...,0)^T$ the boundary condition (inflow at first node).

Conservation implies that if we sum all fluxes of equation $(*)$ multiplying by a vector $\vec{1}=(1,1,...)$ we must obtain the total variation of fluxes $f_T = \sum_if_i$ $$\partial_t f_T+\vec{1}^TD\vec{F}=F(f_0(t))/h$$ since $\vec{1}^TD = (0,0,...,1)$ $$\partial_t f_T = F(f_0(t))/h-\vec{F}(L)/h$$ which is the conservation requirement.

Since this applies for both cases, in which node $0$ or $N+1$ are inflow boundaries (matrix $D$ changes accordingly) we can say that all the sum of all the columns must add to $0$, as shown ($\vec{1}^TD$ must be equal to $(1,0,0,...)$ or $(0,0,....,1)$ if the inflow is at the last or first node respectively).

  • $\begingroup$ I follow your reasoning, but it's not clear to me to what extent it is true that this property (zero sum of columns) of the $D = A^{-1}B$ matrix is implied by the matrix $B$ having that property. $\endgroup$ – Enrico Maria De Angelis May 10 '18 at 11:30

Matrix $\mathbf B$ having columns summing up to $(-1,0,0,...,0,0,+1)$ does not imply that the same happens for $\mathbf A^{-1}\mathbf B$.

Indeed, if $\mathbf A^{-1}\mathbf B$ has that property, than it is correct to use $\mathbf 1^{\text{T}}$ to discretize the integral over the domain. If the scheme used is such that $\mathbf B$ has that property, than the row vector $\textbf 1^{\text{T}} \textbf A$ should be used to discretize the integral.

This appears clearer in the light of R. Knikker's 2009 paper on the Journal of Computational Physics.


Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.