# A NEW CLASS OF FINITE DIFFERENCE METHODS: THE ZIGZAG SCHEMES

Lorenzo Poggioni, Didier Clamond, Yves D'angelo

► To cite this version:

Lorenzo Poggioni, Didier Clamond, Yves D'angelo. A NEW CLASS OF FINITE DIFFERENCE METHODS: THE ZIGZAG SCHEMES. 2025. hal-05075534

HAL Id: hal-05075534

https://hal.science/hal-05075534v1

Preprint submitted on 22 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A NEW CLASS OF FINITE DIFFERENCE METHODS: THE ZIGZAG SCHEMES

LORENZO POGGIONI, DIDIER CLAMOND & YVES D'ANGELO

**ABSTRACT.** We introduce a novel class of finite difference approximations, termed *zigzag schemes*, that employ a hybrid stencil that is neither symmetrical, nor fully one-sided. These zigzag schemes often enjoy more permissive stability constraints and see their coefficients vanish as the order tends to infinity. This property permits the formulation of higher order schemes. An explicit formula is given for both collocated and staggered grids for an arbitrary order and a closed-form expression for the infinite-order scheme is also provided. A linear stability analysis indicates that the zigzag scheme offer a broader range of conditional stability compared to the centred and upwind schemes, sometimes being the only stable scheme. Additionally, the asymmetrical structure of the stencil of zigzag schemes prevents some issues such as the formation of “ghost solutions”. Moreover, implementing zigzag schemes is relatively easy when a code using classical finite differences is available, that is an important feature for well-tested legacy codes. Overall, zigzag schemes provide a compelling alternative for finite differences methods by enabling faster and more stable numerical simulations without sacrificing accuracy or ease of use.

## 1. INTRODUCTION

Introduced by Taylor [23] and Euler [4] in the eighteenth century, finite difference methods remain a cornerstone in the numerical approximation of partial differences equations. These methods are used in a wide array of domains, including the study of immersed boundary problems [5, 15, 18, 19] and the study of hyperbolic conservation laws [9, 12, 20]. Additionally, the local nature of finite difference schemes make them an adapted tool regarding high performance computing, for both CPU [1, 22] and GPU [8, 11, 24]. Traditional schemes, most notably the centred and upwind schemes, are well known for their effectiveness and their simplicity. However, these schemes are often characterised by tight restrictions regarding their stability. These restrictions also tend to worsen as the order of the schemes increases. Thus, hybrid difference schemes have been proposed for convection-diffusion problems [21]. The idea is to use locally a centred or an upwind/downwind schemes according to the characteristics of the problem, and thus to take advantage of the good properties of both schemes. An alternative approach is described in the present paper.

In this work, we introduce a novel finite difference method, termed *zigzag schemes*, designed to relax the aforementioned restrictions while retaining the simplicity of implementation enjoyed by the centred and upwind schemes. This allows the zigzag schemes to be implementable into already well established legacy codes based on finite difference methods such as the direct numerical simulation (DNS) solver H-Allegro [3]. Furthermore, the performance of this new general method does not rely on properties inherent to a specific class of equation. Contrary to the upwind schemes, the coefficients of the zigzag schemes vanish as the order tends to infinity, which allows the use of higher-order schemes. These zigzag schemes are characterised

<footnote>Date: May 22, 2025.</footnote>
by stencils that are neither fully centred nor fully one-sided, but instead represent a hybrid of both. We present an explicit formula for the coefficients of these schemes for arbitrary orders, as well as a closed-form expression for the infinite-order case.

The paper is organised as follows. Sections 2 and 3 present, respectively, the classical centred and onde-sided (upwind) schemes. In section 4, we introduce the new collocated and staggered zigzag schemes. Section 5 is devoted to a linear stability analysis of the zigzag schemes, as well as a comparison with the centred and upwind schemes showcasing situations where the zigzag schemes are the only stable schemes. In the section 6, we highlight a situation where zigzag schemes are able to avoid problems inherent to symmetrical schemes (such as the formation of “ghost solutions”) without being limited to lower orders.

## 2. CLASSICAL CENTRED FINITE DIFFERENCES

### 2.1. Schemes definition.
Consider a function $f$ of one real variable $x$, discretised with a constant step size $\Delta x$. The discrete abscissa are noted as $x_i \stackrel{\text{def}}{=} i\Delta x$ and the discrete value of $f$ at $x_i$ is noted $f_i = f(x_i)$ for $i \in \{1, 2, ..., N\}$. A 2N-th order centred finite approximation of $f'$ and $f''$ can be written [6]

$$f'_{i} = \sum_{j=1}^{N} C_{N}^{j} \frac{f_{i+j} - f_{i-j}}{2j\Delta x}, \qquad (1)$$

$$f''_{i} = \sum_{j=1}^{N} C_{N}^{j} \frac{f_{i+j} - 2f_{i} + f_{i-j}}{(j \Delta x)^2}, \qquad (2)$$

where

$$C_{N}^{j} \stackrel{\text{def}}{=} \frac{2(-1)^{j+1}(N!)^2}{(N+j)!(N-j)!} \qquad (3)$$

Note that as $N \to \infty$, the centred coefficients become [7]

$$C_{\infty}^{j} = 2(-1)^{j+1}, \qquad (4)$$

an infinite-order scheme can then be considered. If the interval of definition of $f$ is discretised using $P \ge 2$ equally spaced points, then necessarily $1 \le N \le \lfloor(P-1)/2\rfloor$ ($\lfloor\cdot\rfloor$ the rounding towards $-\infty$), and appropriate boundary conditions must be taken into account to deduce values at virtual nodes outside the definition domain. Note that the information contained in (1) and (2) is sufficient to approximate derivatives of any order, since higher-order derivatives can be expressed as combinations of first- and second-order derivatives.

### 2.2. Staggered variant.
Instead of using nodes at integer positions $x = j\Delta x$, one may advantageously consider nodes at intermediate positions. The generic staggered centred scheme of order 2N then is [6]:

$$f'_{i} = \sum_{j=1}^{N} \check{C}_{N}^{j} \frac{f_{i+j-1/2} - f_{i-j+1/2}}{(2j-1)\Delta x}, \qquad (5)$$

where the coefficients $\check{C}_N^j$ are defined by

$$\check{C}_{N}^{j} \stackrel{\text{def}}{=} \frac{(-1)^{j+1} 4^{1-2N} (2N)!^2}{(2j-1) N!^2 (N+j-1)! (N-j)!}, \qquad \check{C}_{\infty}^{j} = \frac{4 (-1)^{j+1}}{(2j-1) \pi}. \qquad (6)$$
Note that the staggered centred scheme coefficients decay faster than the coefficients of the collocated centred scheme. The use of such staggered schemes leads to a substantial gain in accuracy over the centred ones with identical $N$ and $\Delta x$ [6].

2.3. Relation to Lanczos' $\sigma$-factors. For a pure frequency of wavenumber $k$, i.e. $f = e^{ikx}$ (with $i^2 = -1$ not to be confused with an index $i$), we have

$$ \left[ \left( e^{ikx} \right)' \right]_i = i k e^{ikx} \sum_{j=1}^{N} C_N^j \frac{\sin(kx_j)}{kx_j} = \sigma_c(k, N) i k e^{ikx}, \quad (7) $$

where

$$ \sigma_c(k, N) \stackrel{\text{def}}{=} \sum_{j=1}^{N} C_N^j \text{sinc} \left( \frac{j\pi k}{k_{\text{max}}} \right) \quad \text{with} \quad \left\{ \begin{array}{ll} k_{\text{max}} & \stackrel{\text{def}}{=} \pi/\Delta x \\ \text{sinc} & \stackrel{\text{def}}{=} \sin(x)/x \\ x_j & = j\Delta x = j\pi/k_{\text{max}} \end{array} \right. \quad (8) $$

Here, $\sigma_c(k, N)$ is called the centred $\sigma$-factor of order $2N$ and $\sigma_c(k, 1)$ is the classical $\sigma$-factor of Lanczos [16]. Note that $\sigma_c$ satisfies the properties: $\sigma_c(k, N) \in \mathbb{R}$, $\sigma_c(k, N) = \sigma_c(-k, N)$, $\sigma_c(0, N) = 1$ and $\sigma_c(k_{\text{max}}, N) = 0$.

For a pure frequency, higher-order derivatives can concisely be defined as

$$ \left[ \left( e^{ikx} \right)^{2n} \right]_i = (ik)^{2n} e^{ikx_i} \sum_{j=1}^{N} C_N^j \text{sinc} \left( \frac{j \pi k}{2 k_{\text{max}}} \right)^{2n}, \quad (9) $$

$$ \left[ \left( e^{ikx} \right)^{(2n+1)} \right]_i = (ik)^{2n+1} e^{ikx_i} \sum_{j=1}^{N} C_N^j \cos\left(\frac{j \pi k}{2 k_{\text{max}}}\right) \text{sinc}\left(\frac{j \pi k}{2 k_{\text{max}}}\right)^{2n+1}, \quad (10) $$

the elementary proof being left to the reader.

2.4. Practical computation of high orders. In practice, the coefficients $C_N^j$ are difficult to handle numerically for the (not so) large $N$ because of the factorials that appear within the formulas, leading to large round-off errors. To avoid both inaccuracies and integer overflows, certain computational techniques should be used.

One trick is to rely on the logarithm of the $\Gamma$-function. Using this approach, one can compute $C_N^j$ for higher values of $N$ faster and more accurately using the optimised Gammaln function [2] and exploiting the relation

$$ \frac{(N!)^2}{(N+j)!(N-j)!} = \exp(2 \log[\Gamma(N+1)] - \log[\Gamma(N+j+1)] - \log[\Gamma(N-j+1)]). $$

However, this approach merely delays the issue by a few orders. Indeed, even using this technique, we found that it was sometimes impossible to go beyond a few hundred in order without overflow. To overcome this limitation, we found that an accurate and robust way to compute $C_N^j$ for any $N \in \mathbb{N}$ is to exploit the relation

$$ \frac{(N!)^2}{(N+j)!(N-j)!} = \prod_{\ell=1}^{j} \left(1 - \frac{j}{N+\ell}\right) = \exp\left[\sum_{\ell=1}^{j} \log 1_p\left(\frac{-j}{N+\ell}\right)\right], \quad (11) $$

where $\log 1_p(x) \stackrel{\text{def}}{=} \log(1+x)$ is specifically designed [10] to be accurate to machine precision as $1+x \approx 1$. To our knowledge, this technique never appeared before in the literature (for computed the finite difference weights). It allows to compute finite differences coefficients of
<table><caption>TABLE 1. Computational time (seconds) of the centred staggered σ-factor using different approaches. For direct computation (left column), we compute $N!$ using recursive multiplication $N(N-1)(N-2)...1$: this approach is limited to fairly low values of $N$. For large $N$, the log1p approach (last column, Eq. (11)) compares to Fornberg's approach, while allowing a direct computation.</caption><thead><tr><th>order N</th><th>direct computation</th><th>Γ approach</th><th>Fornberg</th><th>log1p</th></tr></thead><tbody><tr><td>10</td><td>0.0079</td><td>0.0081</td><td>0.0088</td><td>0.0083</td></tr><tr><td>80</td><td>0.2171</td><td>0.2145</td><td>0.2180</td><td>0.2187</td></tr><tr><td>100</td><td>overflow</td><td>0.2815</td><td>0.2819</td><td>0.2938</td></tr><tr><td>250</td><td>overflow</td><td>0.7886</td><td>0.7923</td><td>0.8014</td></tr><tr><td>300</td><td>overflow</td><td>overflow</td><td>0.9658</td><td>0.9624</td></tr></tbody></table>

large orders accurately within a reasonable time, as illustrated in Table 1 for the computation of the staggered coefficients. Note that this concise formula performs similarly to Fornberg's efficient algorithm [7] (c.f. Table 1).

### 3. CLASSICAL FORWARD/BACKWARD FINITE DIFFERENCES

It is well known that, in certain situations, centred schemes can be unstable. An alternative is to consider stencils that are not symmetric. Fully non-centred (i.e. one-sided) schemes can then be considered. Usual $N$-th order forward/backward schemes for $f'$ and $f''$ approximations can be written as

$$f'_{i} = \sum_{j=1}^{N} D_{N}^{j} \frac{f_{i \pm j} - f_{i}}{\pm j \Delta x}, \qquad f''_{i} = \sum_{j=1}^{N} D_{N}^{j} \frac{f_{i \pm 2j} - 2f_{i \pm j} + f_{i}}{(j \Delta x)^2}, \qquad (12)$$

where

$$D_{N}^{j} \stackrel{\text{def}}{=} \frac{(-1)^{j+1} N!}{j! (N-j)!}. \qquad (13)$$

Thus, for a pure frequency $f = e^{ikx}$, this yields

$$\left[ \left( e^{ikx} \right)' \right]_{i} = ik e^{ikx_i} \sum_{j=1}^{N} D_{N}^{j} \frac{\exp(ikx_{\pm j}) - 1}{ikx_{\pm j}} \qquad (14)$$

$$= ik e^{ikx_i} \sum_{j=1}^{N} D_{N}^{j} \left[ \frac{\sin(kx_{\pm j})}{kx_{\pm j}} + i \frac{1 - \cos(kx_{\pm j})}{kx_{\pm j}} \right]. \qquad (15)$$

The latter relation suggests the introduction of a forward/backward σ-factor of order $N$ which, keeping the notation introduced in (8), is defined by

$$\sigma_{\text{F/B}}(k, N) \stackrel{\text{def}}{=} \sum_{j=1}^{N} D_{N}^{j} \left[ \text{sinc} \left( \frac{j\pi k}{k_{\text{max}}} \right) \pm i \frac{1 - \cos(j\pi k/k_{\text{max}})}{j\pi k/k_{\text{max}}} \right]. \qquad (16)$$

This σ-factor verifies $\sigma_{\text{F/B}}(k, N) \in \mathbb{C}$, $\sigma_{\text{F/B}}(-k, N) = \overline{\sigma_{\text{F/B}}(k, N)}$, $\sigma_{\text{F/B}}(0, N) = 1$, $\sigma_{\text{F/B}}(2k_{\text{max}}, N) = 0$ and $\text{Re}\{\sigma_{\text{F/B}}(k_{\text{max}}, N)\} = 0$, ($\overline{\cdot}$) denoting the complex conjugate of ($\cdot$)).
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Forward-first zigzag scheme stencil used to compute approximation at the j-th node. Dark points are used, while hollow points are not.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 2. Backward-first zigzag scheme stencil used to compute approximation at the j-th node.</figcaption></figure>

One can notice that, as $N \to \infty$, the coefficients $D_N^j$ do not converge since

$$|D_N^j| \sim N^j / j! \quad \text{as} \quad N \to \infty. \qquad (17)$$

Consequently, this approach is restricted to relatively low orders $N$. To overcome these limitations, we propose a novel non-centred scheme.

## 4. ZIGZAG FINITE DIFFERENCES

Rather than solely relying on forward or backward schemes, we consider here a hybrid approach combining both. More specifically, high-order non-centred schemes can be constructed alternating forward and backward differences, starting with either a forward or backward difference. We call these hybrid schemes *forward/backward-first zigzag schemes*.

### 4.1. Definition of zigzag schemes for regular grids.
The figures 1 and 2 illustrate how the stencils of such methods look like. A Nth order zigzag scheme approximating $f'$ and $f''$ is defined by

$$f'_{i} = \sum_{j=1}^{N} Z_{N}^{j} \frac{f_{i \mp j(-1)^{j}} - f_{i}}{\mp(-1)^{j} j \Delta x}, \qquad (18)$$

$$f''_{i} = \sum_{j=1}^{N} Z_{N}^{j} \frac{f_{i \mp 2j(-1)^{j}} - 2f_{i \mp j(-1)^{j}} + f_{i}}{(j \Delta x)^{2}}, \qquad (19)$$

where the coefficients $Z_N^j$ are given analytically by the formulae (20)–(25) (see Table 3 for the numerical values of the first ones).

The coefficients $Z_N^j$ were obtained using a computer algebra system, low-order ones are given in Table 3. From this table, the general formula (20)–(25) for $Z_N^j$ was inferred and
<table><caption>FIGURE 3. First coefficients $Z_N^j$ of the zigzag schemes.</caption><thead><tr><th>N<br>j</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th></tr></thead><tbody><tr><th>1</th><td>1</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>2</th><td><sup>2</sup>&frasl;<sub>3</sub></td><td><sup>1</sup>&frasl;<sub>3</sub></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>3</th><td>1</td><td><sup>1</sup>&frasl;<sub>5</sub></td><td>-<sup>1</sup>&frasl;<sub>5</sub></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>4</th><td><sup>4</sup>&frasl;<sub>5</sub></td><td><sup>2</sup>&frasl;<sub>5</sub></td><td>-<sup>4</sup>&frasl;<sub>35</sub></td><td>-<sup>3</sup>&frasl;<sub>35</sub></td><td></td><td></td><td></td><td></td></tr><tr><th>5</th><td>1</td><td><sup>2</sup>&frasl;<sub>7</sub></td><td>-<sup>2</sup>&frasl;<sub>7</sub></td><td>-<sup>1</sup>&frasl;<sub>21</sub></td><td><sup>1</sup>&frasl;<sub>21</sub></td><td></td><td></td><td></td></tr><tr><th>6</th><td><sup>6</sup>&frasl;<sub>7</sub></td><td><sup>3</sup>&frasl;<sub>7</sub></td><td>-<sup>4</sup>&frasl;<sub>21</sub></td><td>-<sup>1</sup>&frasl;<sub>7</sub></td><td><sup>2</sup>&frasl;<sub>77</sub></td><td><sup>5</sup>&frasl;<sub>231</sub></td><td></td><td></td></tr><tr><th>7</th><td>1</td><td><sup>1</sup>&frasl;<sub>3</sub></td><td>-<sup>1</sup>&frasl;<sub>3</sub></td><td>-<sup>1</sup>&frasl;<sub>11</sub></td><td><sup>1</sup>&frasl;<sub>11</sub></td><td><sup>5</sup>&frasl;<sub>429</sub></td><td>-<sup>5</sup>&frasl;<sub>429</sub></td><td></td></tr><tr><th>8</th><td><sup>8</sup>&frasl;<sub>9</sub></td><td><sup>4</sup>&frasl;<sub>9</sub></td><td>-<sup>8</sup>&frasl;<sub>33</sub></td><td>-<sup>2</sup>&frasl;<sub>11</sub></td><td><sup>8</sup>&frasl;<sub>143</sub></td><td><sup>20</sup>&frasl;<sub>429</sub></td><td>-<sup>8</sup>&frasl;<sub>1287</sub></td><td>-<sup>7</sup>&frasl;<sub>1287</sub></td></tr></tbody></table>

validated via symbolic computations for orders up to five thousands<fnref n="1" />

$$Z_N^1 = \frac{2N!}{(N+1)!} \left\lfloor \frac{N+1}{2} \right\rfloor, \qquad (20)$$

$$Z_N^2 = \frac{2N!}{(N+2)!} \left\lfloor \frac{N}{2} \right\rfloor \left\lfloor \frac{N+1}{2} \right\rfloor, \qquad (21)$$

$$Z_N^j = (-1)^{\lfloor j/2 \rfloor} \frac{2N + j - 2 - (j+2)(-1)^N}{(j/2 - 1/2)!(N+j)!} \frac{(j-2)!N!}{(j/2 - 3/2)!} \left\lfloor \frac{N}{2} \right\rfloor \left\lfloor \frac{N+1}{2} \right\rfloor \qquad (22)$$

$$\times \prod_{\ell=1}^{j/2-3/2} \left\lfloor \frac{N-2\ell-1}{2} \right\rfloor \left\lfloor \frac{N+2\ell+1}{2} \right\rfloor \quad \ell \text{ odd} \ge 3, \qquad (23)$$

$$Z_N^j = (-1)^{1+j/2} \frac{j!N!}{(j/2)!^2 (N+j)!} \left\lfloor \frac{N}{2} \right\rfloor \left\lfloor \frac{N+1}{2} \right\rfloor \left\lfloor \frac{N-3}{2} \right\rfloor \left\lfloor \frac{N+3}{2} \right\rfloor \qquad (24)$$

$$\times \prod_{\ell=1}^{j/2-2} \left\lfloor \frac{N-2\ell-2}{2} \right\rfloor \left\lfloor \frac{N+2\ell+4}{2} \right\rfloor \quad \ell \text{ even} \ge 4, \qquad (25)$$

where $\lfloor \cdot \rfloor$ and $\lceil \cdot \rceil$ are, respectively, the rounding towards $-\infty$ and $+\infty$. Note that these schemes are not centred and, unlike the classical forward/backward schemes, the coefficients $Z_N^j$ converge as $N \to \infty$, i.e.,

$$Z_{\infty}^{2\ell} = -Z_{\infty}^{2\ell+1} = \frac{(-1)^{\ell+1} (2\ell)!}{4^{\ell} (\ell!)^2}. \qquad (26)$$

Furthermore, unlike their classical centred counterparts (4), the $\infty$-order zigzag coefficients $Z_{\infty}^j$ decay towards zero as $j \to \infty$ since $(2\ell)!/4^{\ell}(\ell!)^2 \sim 1/\sqrt{\pi\ell}$ as $\ell \to \infty$. A visual summary of the evolution of the coefficients' magnitude from the forward, centred and zigzag schemes is available in Appendix A.

<footnote n="1">¹Despite our efforts, we did not manage to prove the formula rigorously, nor to simplify it.</footnote>
4.2. **Zigzag σ-factors.** Considering a pure frequency $f = e^{ikx}$, we have

$$ \begin{aligned} \left( e^{ikx} \right)'_i &= ik e^{ikx} \sum_{j=1}^{N} Z_N^j \frac{\exp(ik x_{\mp j(-1)^j}) - 1}{ik x_{\mp j(-1)^j}} \\ &= ik e^{ikx} \sum_{j=1}^{N} Z_N^j \left[ \frac{\sin(k x_{\mp j(-1)^j})}{k x_{\mp j(-1)^j}} + i \frac{1 - \cos(k x_{\mp j(-1)^j})}{k x_{\mp j(-1)^j}} \right], \end{aligned} \quad (27) $$

which, once again using the notations introduced in (8), suggests the zigzag σ-factor

$$ \sigma_z(k, N) \stackrel{\text{def}}{=} \sum_{j=1}^{N} Z_N^j \left[ \text{sinc} \left( \frac{j \pi k}{k_{\text{max}}} \right) \mp (-1)^j i \frac{1 - \cos (j \pi k / k_{\text{max}})}{j \pi k / k_{\text{max}}} \right], \quad (28) $$

where the upper (resp. lower) sign corresponds to the forward-first (resp. backward-first) scheme.

A closed-form of the $\infty$-order zigzag σ-factor can be obtained. For the forward-first (the backward-first being the complex conjugate of the forward-first), we have

$$ \begin{aligned} \sigma_z(k, \infty) &= \sum_{\ell=1}^{\infty} \frac{4(-1)^{\ell+1}(2\ell-2)! \exp((2\ell-1)i\pi k / k_{\text{max}}) - 1}{4^{\ell}(\ell-1)!^2 (2\ell-1)i\pi k / k_{\text{max}}} \\ &\quad - \sum_{\ell=1}^{\infty} \frac{(-1)^{\ell+1}(2\ell)! \exp(-2\ell i\pi k / k_{\text{max}}) - 1}{4^{\ell}(\ell!)^2 2\ell i\pi k / k_{\text{max}}} \\ &= \frac{k_{\text{max}}}{i\pi k} \log \left( \frac{\exp(i\pi k / k_{\text{max}}) + \sqrt{1 + \exp(2i\pi k / k_{\text{max}})}}{1 + \sqrt{1 + \exp(-2i\pi k / k_{\text{max}})}} \right), \end{aligned} \quad (29) $$

where we have exploited the two relations

$$ \sum_{\ell=1}^{\infty} \frac{4(-1)^{\ell+1}(2\ell-2)!}{4^{\ell}(\ell!)^2} = \arcsinh(X) = \log \left( X + \sqrt{1+X^2} \right), \quad (30) $$

$$ \sum_{\ell=1}^{\infty} \frac{(-1)^{\ell+1}(2\ell)!}{4^{\ell}(\ell!)^2} \frac{X^{2\ell}}{2\ell} = \log \left( \frac{1 + \sqrt{1 + X^2}}{2} \right). \quad (31) $$

4.3. **Zigzag schemes for staggered grids.** Zigzag schemes can be also considered for staggered grids, the staggered variant of these zigzag schemes being obviously

$$ f'_i = \sum_{j=1}^{N} \tilde{Z}_N^j \frac{f_{i \mp (j-1/2)(-1)^j} - f_i}{\mp (-1)^j (j-1/2) \Delta x}, \quad (32) $$

where the coefficients $\tilde{Z}_N^j$ are defined by

$$ \tilde{Z}_N^j \stackrel{\text{def}}{=} \frac{(-1)^{\lceil 1+j/2 \rceil} (2N-1)!}{(2j-1) 8^{N-1} (N-1)!} \left( \left\lfloor \frac{N+j-1}{2} \right\rfloor! \left\lfloor \frac{N-j}{2} \right\rfloor! \right)^{-1}. \quad (33) $$

Again, this formula was inferred from numerically computed coefficients, then verified using a computer algebra system up to order one thousand.
<figure><img src="image_4.png" /><figcaption>FIGURE 4. Stability region as a function of $\Lambda = -i\kappa\pi\lambda\sigma(\kappa)$ in the complex plane. The darker area is where the scheme is stable, and the dotted line represent the critical values where $|G| = 1$.</figcaption></figure>

## 5. STABILITY ANALYSIS ON LINEAR ADVECTION EQUATION

In order to compare the performances of the zigzag schemes to the classical centred and forward/backward ones, we conduct a stability analysis<fnref n="2" /> for the reference transport equation

$$u_t + c u_x = 0, \quad u(0, x) = u_0(x), \quad (t, x) \in \mathbb{R}_+ \times \mathcal{T}, \qquad (34)$$

where $\mathcal{T}$ is a periodic, one-dimensional domain and $c \in \mathbb{R}$ is a wave celerity. The exact solution of this equation is $u(t, x) = u_0(x - ct)$. Looking for solutions of the form $u(t, x) = U(t)e^{ikx}$ (with $k \in \mathbb{R}$) and using a finite-difference scheme for the $x$-derivative, approximation of Eq. (34) yields the semi-discrete equation

$$\frac{dU}{dt} + ikc\sigma(k)U = 0, \qquad (35)$$

where $\sigma$ denotes the $\sigma$-factor of the scheme used to approach the spatial derivative. Note that the scheme is centred, $\sigma(k) \in \mathbb{R}$. The temporal variable $t$ remains to be discretised and we investigate below several possibilities.

### 5.1. Euler explicit scheme.
We start with the first-order explicit Euler scheme $dU(t)/dt \approx [U(t + \Delta t) - U(t)]/\Delta t$ which, after substitution in (35), leads to the fully discrete equation

$$\frac{U(t + \Delta t) - U(t)}{\Delta t} + ikc\sigma(k)U(t) = 0. \qquad (36)$$

With $U = a \exp(i\omega t)$ ($a \in \mathbb{R}$), one obtains $e^{i\omega\Delta t} = 1 - ik\Delta t c\sigma(k)$. Introducing the amplification factor $G \stackrel{\text{def}}{=} U(t + \Delta t)/U(t) = e^{i\omega\Delta t}$, the von Neumann stability condition [14] $||G||_\infty \le 1$ shows that the solution is bounded if and only if (for any wavenumber $k \in \mathbb{R}$)

$$|1 - i\Delta t c\sigma(k)| \le 1. \qquad (37)$$

<footnote n="2">²A comparison of the numerical diffusion introduced by these schemes is also available in the Annex B.</footnote>
Let be the Nyquist wavenumber $k_{\max} \stackrel{\text{def}}{=} \pi/\Delta x$, the scaled wavenumber $\kappa = k/k_{\max}$ and the “CFL-like” stability-number $\lambda \stackrel{\text{def}}{=} c\Delta t/\Delta x$. Then (37) is rewritten

$$|1 - i\pi\lambda\kappa\sigma(\kappa)| \le 1, \quad \forall |\kappa| \le 1, \qquad (38)$$

this condition allowing the computation of the stability region of the general explicit Euler scheme of order 1 (Figure 4). When centred schemes are considered, $\sigma$ is real and the condition (38) yields

$$(\pi\lambda\kappa\sigma)^2 \le 0, \qquad (39)$$

which is never satisfied. In other words, this means that centred schemes are unconditionally unstable, a very well-known result [17].

5.1.1. Forward finite differences scheme of order 1. With the spatial derivative approximated by a first-order forward finite difference scheme, the corresponding $\sigma$-factor is

$$\sigma(\kappa) = \frac{e^{i\pi\kappa} - 1}{i\pi\kappa}. \qquad (40)$$

Substituting (40) into the condition (38), we obtain the stability condition

$$\forall \kappa \in [-1, 1], \quad |1 + \lambda - \lambda e^{i\pi\kappa}| \le 1 \quad \Rightarrow \quad \lambda(1 + \lambda) \le 0 \quad \Rightarrow \quad -1 \le \lambda \le 0. \qquad (41)$$

In order to satisfy the above condition, the wave celerity $c$ needs to be negative (when the time increases). This means that the scheme must be “upwind” [17]. Condition (41) also implies that $|\Delta x/\Delta t| \ge c$, meaning that the “numerical speed” must be greater than the physical speed $c$. (If the celerity of the wave is positive, then a backward scheme must be considered instead.) The corresponding stability region is given in appendix C.

5.1.2. Forward finite differences scheme of order 2. For this scheme, we have

$$\sigma(\kappa) = i \frac{e^{2i\pi\kappa} - 4e^{i\pi\kappa} + 3}{2\pi\kappa}, \qquad (42)$$

and the substitution of (42) into (38) yields

$$\left| 1 + \frac{\lambda}{2} \left( e^{2i\pi\kappa} - 4e^{i\pi\kappa} + 3 \right) \right| \le 1. \qquad (43)$$

After doing computations similar to those done in (41), we show that (43) implies $\lambda = 0$. Hence, we can conclude that the second-order purely forward scheme (with first-order explicit Euler temporal scheme) is stable only for some wavenumbers, but never for all wavenumbers. Therefore, it is unconditionally unstable.

5.1.3. Zigzag finite differences scheme of order 2. Consider now the zigzag forward-first finite difference scheme of order 2 (recall that the first-order zigzag scheme coincides with the first-order purely forward scheme); this scheme can be defined via its $\sigma$-factor

$$\sigma(\kappa) = i \frac{e^{-2i\pi\kappa} - 4e^{i\pi\kappa} + 3}{6\pi\kappa}. \qquad (44)$$

Substituting $\sigma$ into (38) gives the condition

$$\left| 1 + \frac{\lambda}{6} \left( e^{-2i\pi\kappa} - 4e^{i\pi\kappa} + 3 \right) \right| \le 1. \qquad (45)$$
$$ \begin{array}{c|c} 0 & \\ \frac{1}{2} & \frac{1}{2} \\ \hline 0 & 1 \end{array} $$

<figure><img src="image_5.png" /><figcaption>FIGURE 5. Stability region in function of $\Lambda = -i\kappa\pi\lambda\sigma(\kappa)$ in the complex plane. The darker area contains the values of $\Lambda$ for which the scheme is stable. The dotted line represent the critical values where $|G| = 1$.</figcaption></figure>

After some computations we can show that this also implies $\lambda = 0$. Hence, as for the pure forward scheme, the zigzag scheme of order two (with explicit Euler) is unconditionally unstable (c.f. Appendix C) due to the lack of stability for low wavenumbers.

5.2. **Second-order Runge-Kutta scheme.** To increase the order of the time scheme, we consider the second-order Runge-Kutta method given by the Butcher's tableau in Table 2. Doing so, the semi-discrete equation (35) becomes the fully discrete equation

$$ U(t + \Delta t) = \left[ 1 + \Delta t A + \frac{1}{2}(\Delta t A)^2 \right] U(t), \quad A \stackrel{\text{def}}{=} -ik c \sigma(k), \qquad (46) $$

with the von Neumann stability condition

$$ \left| 1 + \Delta t A + \frac{1}{2}(\Delta t A)^2 \right| \le 1 \quad \forall k. \qquad (47) $$

Using the notations introduced in section 1, the condition (47) becomes

$$ \left| 1 - i\pi\kappa\lambda\sigma - \frac{1}{2}(\pi\kappa\lambda\sigma)^2 \right| \le 1 \quad \forall |\kappa| \le 1, \qquad (48) $$

the corresponding stability region being depicted in figure 5. For centred schemes, $\sigma$ is real and this condition yields

$$ (\pi\lambda\kappa\sigma)^4 \le 0. \qquad (49) $$
It cannot be satisfied meaning that centred schemes are unconditionally unstable, like the explicit Euler scheme.

5.2.1. Forward finite differences scheme of order 1. The σ-factor corresponding to this case, i.e.

$$ \sigma(\kappa) = \frac{e^{i\pi\kappa} - 1}{i\pi\kappa}, \qquad (50) $$

substituted into (48) gives the condition

$$ \left| 1 - \lambda (e^{i\pi\kappa} - 1) + \frac{1}{2}\lambda^2 (e^{i\pi\kappa} - 1)^2 \right| \le 1 \quad \Rightarrow \quad -1 \le \lambda^3 + 2\lambda^2 + 2\lambda \le 0. \qquad (51) $$

$$ \Rightarrow \quad \lambda \in [-1, 0]. \qquad (52) $$

Hence, the forward finite differences scheme of order one is conditionally stable, requiring the stability-number λ to be in [-1, 0]. We note that, like in the case of the explicit Euler scheme, the scheme must be “upwind” to be stable.

5.2.2. Forward finite differences scheme of order 2. In this case, the σ-factor

$$ \sigma(\kappa) = i \frac{e^{2i\pi\kappa} - 4e^{i\pi\kappa} + 3}{2\pi\kappa}, \qquad (53) $$

substituted into the von Neumann stability condition (48) yields

$$ \left| 1 + \frac{\lambda}{2} (e^{2i\pi\kappa} - 4e^{i\pi\kappa} + 3) + \frac{\lambda^2}{8} (e^{2i\pi\kappa} - 4e^{i\pi\kappa} + 3)^2 \right| \le 1. \qquad (54) $$

A study of the condition (54) shows that $\lambda \in [-\frac{1}{2}, 0]$. Thus, the second-order forward finite differences scheme is conditionally stable. Similarly, for the first-order scheme, the strictest constraint on the stability-number is imposed by the Nyquist wavenumber (see Annex D). Note that the stability strip is half as wide as that of the first-order.

5.2.3. Forward finite differences scheme of order 3. For this scheme, the σ-factor is

$$ \sigma(\kappa) = -\frac{i}{6\pi\kappa} (2e^{3i\pi\kappa} - 9e^{2i\pi\kappa} + 18e^{i\pi\kappa} - 11), \qquad (55) $$

which gives

$$ \left| 1 - \frac{\lambda}{6} (2e^{3i\pi\kappa} - 9e^{2i\pi\kappa} + 18e^{i\pi\kappa} - 11) + \frac{\lambda^2}{72} (2e^{3i\pi\kappa} - 9e^{2i\pi\kappa} + 18e^{i\pi\kappa} - 11)^2 \right| \le 1. \qquad (56) $$

A study of this condition reveals that $\lambda = 0$. This scheme is thus unstable. Moreover, numerical observations seem to show that upwind schemes of orders 3 and above are unstable.

5.2.4. Zigzag finite differences scheme of order 2. The σ-factor corresponding to this scheme is

$$ \sigma(\kappa) = i \frac{e^{-2i\pi\kappa} - 4e^{i\pi\kappa} + 3}{6\pi\kappa}, \qquad (57) $$

which once injected in (48) yields

$$ \left| 1 + \frac{\lambda}{6} (e^{-2i\pi\kappa} - 4e^{4i\pi\kappa} + 3) + \frac{\lambda^2}{72} (e^{-2i\pi\kappa} - 4e^{4i\pi\kappa} + 3)^2 \right| \le 1. \qquad (58) $$
This condition implies $\lambda P(\lambda) \le 0$, with $P$ defined by

$$ P(\lambda) \stackrel{\text{def}}{=} \lambda^3 \left[ 64 \cos(\kappa\pi)^4 + 80 \cos(\kappa\pi)^3 + 105 \cos(\kappa\pi)^2 + 50 \cos(\kappa\pi) + 25 \right] \\ + \lambda^2 \left[ -96 \cos(\kappa\pi)^3 + 36 \cos(\kappa\pi)^2 + 60 \right] \\ + \lambda \left[ \cos(\kappa\pi)^2 - 144 \cos(\kappa\pi) + 72 \right] + 216. $$

The polynomial $P$ possesses a global extremum equal to $324\lambda^3 + 216$ at $\kappa = 0$ for any $\lambda$. This extremum is null for the critical $\lambda_c = -\sqrt[3]{2} (\sqrt[3]{3})^2/3 \approx -0.8736$. It corresponds to a minimum of $P$ if and only if $\lambda \ge \lambda_c$, and it is a maximum if and only if $\lambda \le \lambda_c$. Hence, this scheme is conditionally stable, under the condition $\lambda_c \le \lambda \le 0$.

The condition on the stability-number for this scheme is less restrictive than that of the forward-scheme for the same order. Indeed, the stability strip $[-1/2, 0]$ is included in $[\lambda_c, 0]$. Note that this time, the most restrictive stability condition is imposed by the null wavenumber.

5.2.5. *Zigzag finite differences scheme of order 3.* The $\sigma$-factor of this scheme is

$$ \sigma(\kappa) = - \frac{2e^{3i\pi\kappa} + 3e^{-2i\pi\kappa} - 30e^{i\pi\kappa} + 25}{30i\pi\kappa} \qquad (59) $$

Which after substitution in the von Neumann condition (48) gives

$$ \left| 1 + \frac{\lambda}{30} \left( 2e^{3i\kappa\pi} + 3e^{-2i\kappa\pi} - 30e^{i\kappa\pi} + 25 \right) + \frac{\lambda^2}{1800} \left( 2e^{3i\kappa\pi} + 3e^{-2i\kappa\pi} - 30e^{i\kappa\pi} + 25 \right)^2 \right| \le 1. $$

A similar procedure to the one followed for the condition (58) gives $-\frac{15}{14} \le \lambda \le 0$. This scheme remains conditionally stable and it is the first one allowing $|\lambda| > 1$ thus far. Note that for the same order, the zigzag-scheme is conditionally stable whilst the equivalent forward-scheme is not.

5.2.6. *Zigzag finite differences scheme of order 4.* For this scheme,

$$ \sigma(\kappa) = \frac{336 e^{i\kappa\pi} - 245 - 84 e^{-2i\kappa\pi} - 16 e^{3i\kappa\pi} + 9 e^{-4i\kappa\pi}}{420i\kappa\pi}. \qquad (60) $$

Substituting this $\sigma$-factor in (48) yields $\lambda = 0$. This scheme is thus unconditionally unstable. Numerically, we observe that the third-order seems to be the highest order allowing the stability of the zigzag-schemes while using (RK2).

5.3. **Summary on the stability of higher order schemes.** For orders larger than three, the computations of the stability regions become highly tedious. Hence, for these higher orders, we proceeded by computing the stability-number numerically with an accuracy of $\pm 10^{-4}$. In the case of staggered zigzag schemes, we extrapolated the results for the infinite order from the order 300. A summary of all the critical stability-number values can be found in the tables 6,7,8,9, 10,11 and 12 and the corresponding stability regions are in appendix C to I. Note that tables 7, 10 and 11, highlight situations where the zigzag schemes are the only ones to be conditionally stable. Furthermore, these tables also reveal that zigzag schemes consistently exhibit greater stability for even orders. This alternating pattern is expected, as the stencils of zigzag schemes are asymmetrical by design. Consequently, one side of the stencil
<table><caption>FIGURE 6. Highest stability-number $|\lambda| = |c\Delta t/\Delta x|$ allowed for the forward Euler integration scheme in time with spatial schemes of orders $N = 1$ and $N = 2$.</caption><thead><tr><th>scheme</th><th>N = 1</th><th>N &ge; 2</th></tr></thead><tbody><tr><td>centred</td><td>0</td><td>0</td></tr><tr><td>centred staggered</td><td>0</td><td>0</td></tr><tr><td>forward</td><td>1</td><td>0</td></tr><tr><td>zigzag</td><td>1</td><td>0</td></tr><tr><td>zigzag staggered</td><td>1</td><td>0</td></tr></tbody></table>

<table><caption>FIGURE 7. Highest stability-number $|\lambda| = |c\Delta t/\Delta x|$ allowed for RK2 time integrator coupled with spatial schemes of orders $N = 1, 2, 3$ and orders $N \ge 4$. For a given order $N$ (i.e. for a given column), the zigzag schemes enjoy the best stability constraints.</caption><thead><tr><th>scheme</th><th>N = 1</th><th>N = 2</th><th>N = 3</th><th>N &ge; 4</th></tr></thead><tbody><tr><td>centred</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>centred staggered</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>forward</td><td>1</td><td>1/2</td><td>0</td><td>0</td></tr><tr><td>zigzag</td><td>1</td><td>√2(√3)²/3 ≈ 0.8736</td><td>15/14 ≈ 1.0714</td><td>0</td></tr><tr><td>zigzag staggered</td><td>1.5436</td><td>1.2599</td><td>1.7099</td><td>0</td></tr></tbody></table>

<table><caption>FIGURE 8. Highest stability-number $|\lambda| = |c\Delta t/\Delta x|$ allowed for RK3 time integrator with spatial schemes of orders $N = 1$ to 6, and $N = \infty$ (extrapolated from order 300 in the case of the zigzag staggered). For a given order $N$ (i.e. for a given column), the zigzag schemes enjoy the best stability constraints.</caption><thead><tr><th>scheme</th><th>N = 1</th><th>N = 2</th><th>N = 3</th><th>N = 4</th><th>N = 5</th><th>N = 6</th><th>N = &infin;</th></tr></thead><tbody><tr><td>centred</td><td></td><td>1.7320</td><td></td><td>1.2622</td><td></td><td>1.0920</td><td>0.5533</td></tr><tr><td>centred staggered</td><td></td><td>1.7320</td><td></td><td>1.4845</td><td></td><td>1.3949</td><td>1.1282</td></tr><tr><td>forward</td><td>1.2563</td><td>0.6280</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>zigzag</td><td>1.2563</td><td>1.8845</td><td>1.3461</td><td>1.6490</td><td>1.3741</td><td>1.5727</td><td>1.1037</td></tr><tr><td>zigzag staggered</td><td>1.6791</td><td>2.5187</td><td>1.9374</td><td>2.3199</td><td>2.0112</td><td>2.2569</td><td>2.1407</td></tr></tbody></table>

always incorporates more points than the other. However, the enhanced stability observed at even orders is unexpected. Indeed, at even orders, the backward side of the stencil will contain more points and thus, should contribute to a greater extent. This is counter-intuitive, given the well-known result that backward schemes are unconditionally unstable for the transport equation with negative velocity.

## 6. PREVENTION OF GHOST SOLUTIONS

In the recent work, Inglard et al. [13] showcased the formation of ghost solutions when solving the advection equation (34) using centred schemes with von Neumann boundary conditions. The term “ghost solution” designates the anomalous formation of a periodic
<table><caption>FIGURE 9. Highest stability-number $|\lambda| = |c\Delta t/\Delta x|$ allowed for a RK4 time integrator with spatial schemes of orders $N = 1$ to 6, and $N = \infty$ (extrapolated from order 300 in the case of the zigzag staggered).</caption><thead><tr><th>scheme</th><th>N = 1</th><th>N = 2</th><th>N = 3</th><th>N = 4</th><th>N = 5</th><th>N = 6</th><th>N = ∞</th></tr></thead><tbody><tr><td>centred</td><td></td><td>2.8284</td><td></td><td>2.0612</td><td></td><td>1.7834</td><td>0.9035</td></tr><tr><td>centred staggered</td><td></td><td>2.8284</td><td></td><td>2.4243</td><td></td><td>2.2779</td><td>1.8425</td></tr><tr><td>forward</td><td>1.3926</td><td>0.6963</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>zigzag</td><td>1.3926</td><td>2.0889</td><td>1.4921</td><td>1.8278</td><td>1.5232</td><td>1.7433</td><td>1.5800</td></tr><tr><td>zigzag staggered</td><td>1.9122</td><td>2.8683</td><td>2.2064</td><td>2.6419</td><td>2.2904</td><td>2.5701</td><td>2.4378</td></tr></tbody></table>

<table><caption>FIGURE 10. Highest stability-number $|\lambda| = |c\Delta t/\Delta x|$ allowed for a RK5 time integrator with spatial schemes of orders $N = 1$ to 7. For a given order $N$ (i.e. for a given column), the zigzag schemes enjoy the best stability constraints.</caption><thead><tr><th>scheme</th><th>N = 1</th><th>N = 2</th><th>N = 3</th><th>N = 4</th><th>N = 5</th><th>N = 6</th><th>N = 7</th></tr></thead><tbody><tr><td>centred</td><td></td><td>0</td><td></td><td>0</td><td></td><td>0</td><td></td></tr><tr><td>centred staggered</td><td></td><td>0</td><td></td><td>0</td><td></td><td>0</td><td></td></tr><tr><td>forward</td><td>1.6085</td><td>0.8042</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>zigzag</td><td>1.6085</td><td>2.4127</td><td>1.7234</td><td>2.0100</td><td>1.7592</td><td>0.4542</td><td>0</td></tr><tr><td>zigzag staggered</td><td>2.2646</td><td>3.3969</td><td>2.6129</td><td>3.1287</td><td>2.7123</td><td>0.7234</td><td>0.8121</td></tr></tbody></table>

<table><caption>FIGURE 11. Highest stability-number $|\lambda| = |c\Delta t/\Delta x|$ allowed for RK6 time integrator with spatial schemes of orders $N = 1$ to 7. For a given order $N$ (i.e. for a given column), the zigzag schemes enjoy the best stability constraints.</caption><thead><tr><th>scheme</th><th>N = 1</th><th>N = 2</th><th>N = 3</th><th>N = 4</th><th>N = 5</th><th>N = 6</th><th>N = 7</th></tr></thead><tbody><tr><td>centred</td><td></td><td>0</td><td></td><td>0</td><td></td><td>0</td><td></td></tr><tr><td>centred staggered</td><td></td><td>0</td><td></td><td>0</td><td></td><td>0</td><td></td></tr><tr><td>forward</td><td>1.7767</td><td>0.88824</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>zigzag</td><td>1.7767</td><td>2.6650</td><td>1.9035</td><td>2.3319</td><td>1.9431</td><td>2.2240</td><td>0</td></tr><tr><td>zigzag staggered</td><td>2.4915</td><td>3.7373</td><td>2.8749</td><td>3.4423</td><td>2.9842</td><td>3.3488</td><td>3.0350</td></tr></tbody></table>

<table><caption>FIGURE 12. Highest stability-number $|\lambda| = |c\Delta t/\Delta x|$ allowed for a RK7 time integrator with spatial schemes orders $N = 1$ to 6, and $N = \infty$ (extrapolated from order 300 in the case of the zigzag staggered). For a given order $N$ (i.e. for a given column), the zigzag schemes enjoy the best stability constraints.</caption><thead><tr><th>scheme</th><th>N = 1</th><th>N = 2</th><th>N = 3</th><th>N = 4</th><th>N = 5</th><th>N = 6</th><th>N = ∞</th></tr></thead><tbody><tr><td>centred</td><td></td><td>1.7644</td><td></td><td>1.2857</td><td></td><td>1.1124</td><td>0.5638</td></tr><tr><td>centred staggered</td><td></td><td>1.7644</td><td></td><td>1.5122</td><td></td><td>1.4209</td><td>1.1493</td></tr><tr><td>forward</td><td>1.9770</td><td>0.9884</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>zigzag</td><td>1.9770</td><td>2.9656</td><td>2.1182</td><td>2.5949</td><td>2.1624</td><td>2.3282</td><td>1.1198</td></tr><tr><td>zigzag staggered</td><td>2.7440</td><td>4.1160</td><td>3.1661</td><td>3.791</td><td>3.2865</td><td>3.6881</td><td>2.5618</td></tr></tbody></table>
solution, despite imposing von Neumann conditions on the boundaries of the domain (see Fig. 13). Under the same conditions formulated in [13], the replacement of the centred scheme by a one-sided scheme, or by a zigzag scheme completely prevents the formation of ghost solutions. However, as mentioned in section 3, the coefficients of the one-sided scheme diverge as the order tends to infinity. This property limits the use of one-sided schemes to low orders. The zigzags schemes do not suffer such a restriction. In other words, the zigzag schemes allow the use of high-order schemes using von Neumann boundary conditions, without forming ghost solutions. A comparison of the centred and the zigzag schemes with von Neumann conditions is provided in Fig. 13.

## 7. CONCLUSION

In this work, we introduced a novel class of finite differences approximations, referred to as zigzag schemes. These schemes employ hybrid stencils that are neither fully centred nor fully one-sided, leading to encouraging results. Specifically, as summarised in tables 6 to 12, zigzag schemes exhibit not only a broader range of conditional stability but also, in some cases, a significantly higher critical stability-number, sometimes exceeding that of standard schemes by a factor of two or more. Moreover, the coefficients $Z_N^j$ of the zigzag schemes vanish as the order tends to infinity (see Eq. (26)). The zigzag schemes thus allow the use of higher order schemes. Overall, zigzag schemes offer a compelling alternative for finite difference simulations. They enable faster computations at comparable accuracy levels, while retaining the simplicity of implementation characteristic of conventional schemes. Furthermore, a staggered variant of the zigzag schemes was presented in section 4.1, alongside a closed-form expression for the stencil coefficients at infinite order.

Moreover, the asymmetrical structure of the zigzag schemes prevents the formation of some issues arising when using symmetrical (centred) schemes. This was illustrated in section 6 where we showed that zigzags schemes are able to prevent the formation of anomalous “ghost solutions”.

Lastly, we remark that the zigzag schemes are defined by a signature following the pattern `---+---+---+(or --+-+--+)`. This choice was made because it seemed to be the most “natural” and “easier” way to blend centred and one-sided schemes together. However, other schemes employing other (potentially random) signatures could be investigated.

## REFERENCES

[1] S. Abide, M. S. Binous, and B. Zeghmati. An efficient parallel high-order compact scheme for the 3D incompressible Navier-Stokes equations. *Int. J. Comput. Fluid Dyn.*, 31(4-5): 214-229, 2017. ISSN 1061-8562,1029-0257. doi: 10.1080/10618562.2017.1326592. URL https://doi.org/10.1080/10618562.2017.1326592.

[2] M. Abramowitz and I. A. Stegun. *Handbook of Mathematical Functions*. Dover, 1965.

[3] CORIA-CFD. H-allegro. high-order finite difference code that solves the unsteady compressible reacting navier-stokes equations system on structured cartesian meshes. URL https://www.coria-cfd.fr/index.php/H-Allegro.

[4] L. Euler. *Institutiones calculi integralis*. Acad. Imper. scientiarum, 1768.

[5] E. A. Fadlun, R. Verzicco, P. Orlandi, and J. Mohd-Yusof. Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations. *J. Comput.*
<figure><img src="image_6.png" /><figcaption>FIGURE 13. Numerical simulations of the transport equation (34) using the second order implicit centred scheme of [13] (left column) and the second order implicit zigzag backward-first scheme (right column) at times $t = 0, 0.5, 1, 1.5, 2$ and 4 using 1000 points to discretise the space.</figcaption></figure>
*Phys.*, 161(1):35–60, 2000. ISSN 0021-9991,1090-2716. doi: 10.1006/jcph.2000.6484.
URL https://doi.org/10.1006/jcph.2000.6484.

[6] B. Fornberg. High-order finite differences and the pseudospectral method on staggered grids. *SIAM J. Numer. Anal.*, 27(4):904–918, 1990.

[7] B. Fornberg. *A Practical Guide to Pseudospectral Methods*. Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, 1998. ISBN 9780521645645.

[8] M. Giles, E. László, I. Reguly, J. Appleyard, and J. Demouth. Gpu implementation of finite difference solvers. In *2014 Seventh Workshop on High Performance Computational Finance*, pages 1–8, 2014. doi: 10.1109/WHPCF.2014.10.

[9] E. Godlewski and P.-A. Raviart. *Numerical approximation of hyperbolic systems of conservation laws*, volume 118 of *Applied Mathematical Sciences*. Springer-Verlag, New York, 1996. ISBN 0-387-94529-6. doi: 10.1007/978-1-4612-0713-9. URL https://doi.org/10.1007/978-1-4612-0713-9.

[10] D. Goldberg. What every computer scientist should know about floating-point arithmetic. *ACM Comp. Surv.*, 23(1):5–48, 1991.

[11] B. Hamilton and C. Webb. Room acoustics modelling using gpu-accelerated finite difference and finite volume methods on a face-centered cubic grid. 01 2013.

[12] A. Harten. High resolution schemes for hyperbolic conservation laws. *J. Comput. Phys.*, 49(3):357–393, 1983. ISSN 0021-9991,1090-2716. doi: 10.1016/0021-9991(83)90136-5.
URL https://doi.org/10.1016/0021-9991(83)90136-5.

[13] M. Inglard, F. Lagoutière, and H. H. Rugh. Ghost solutions with centered schemes for one-dimensional transport equations with Neumann boundary conditions. *Annales de la Faculté des sciences de Toulouse : Mathématiques*, Ser. 6, 29(4):927–950, 2020. doi: 10.5802/afst.1650.

[14] E. Isaacson and H. B. Keller. *Analysis of Numerical Methods*. Dover, 1994.

[15] M.-C. Lai and C. S. Peskin. An immersed boundary method with formal second-order accuracy and reduced numerical viscosity. *J. Comput. Phys.*, 160(2):705–719, 2000. ISSN 0021-9991,1090-2716. doi: 10.1006/jcph.2000.6483. URL https://doi.org/10.1006/jcph.2000.6483.

[16] C. Lanczos. *Applied Analysis*. Dover Books on Mathematics. Dover Publications, 1988. ISBN 9780486656564.

[17] R. J. LeVeque. *Finite Volume Methods for Hyperbolic Problems*. Cambridge Texts in Applied Mathematics. Cambridge University Press, 2002.

[18] R. J. LeVeque and Z. L. Li. The immersed interface method for elliptic equations with discontinuous coefficients and singular sources. *SIAM J. Numer. Anal.*, 31(4):1019–1044, 1994. ISSN 0036-1429. doi: 10.1137/0731054. URL https://doi.org/10.1137/0731054.

[19] C. S. Peskin. The immersed boundary method. *Acta Numer.*, 11:479–517, 2002. ISSN 0962-4929,1474-0508. doi: 10.1017/S0962492902000077. URL https://doi.org/10.1017/S0962492902000077.

[20] G. A. Sod. A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws. *J. Comput. Phys.*, 27(1):1–31, 1978. ISSN 0021-9991,1090-2716. doi: 10.1016/0021-9991(78)90023-2. URL https://doi.org/10.1016/0021-9991(78)90023-2.

[21] D. B. Spalding. A novel finite-difference formulation for differential expression involving both first and second derivatives. *Int. J. Numer. Methods Eng.*, 4:551–559, 1972.
[22] L. Szustak, K. Halbiniak, A. Kulawik, J. Wrobel, and P. Gepner. Toward parallel modeling of solidification based on the generalized finite difference method using Intel Xeon Phi. In *Parallel processing and applied mathematics. Part I*, volume 9573 of *Lecture Notes in Comput. Sci.*, pages 411–422. Springer, [Cham], 2016. ISBN 978-3-319-32149-3; 978-3-319-32148-6. doi: 10.1007/978-3-319-32149-3\_39. URL https://doi.org/10.1007/978-3-319-32149-3_39.

[23] B. Taylor. *Methodus incrementorum directa inversa. Auctore Brook Taylor, LL. D. and Regiae Societatis Secretario*. Typis Pearsonianis Prostant apud Gul. Innys ad Insignia Principis in Coemeterio Paulino MDCCXV, Londini, 1715. URL https://books.google.com/books?id=iXN1xgEACAAJ.

[24] C.-C. Ye, P.-J.-Y. Zhang, Z.-H. Wan, R. Yan, and D.-J. Sun. Accelerating cfd simulation with high order finite difference method on curvilinear coordinates for modern gpu clusters. *Advances in Aerodynamics*, 4(1), Feb. 2022. ISSN 2524-6992. doi: 10.1186/s42774-021-00098-3. URL http://dx.doi.org/10.1186/s42774-021-00098-3.
<figure><img src="image_7.png" /><figcaption>(A) Magnitude of the coefficients $D_N^j$ for the forward schemes (c.f. equation (13)).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(B) Magnitude of the coefficients $C_N^j$ for the centred schemes (c.f. equation (3)).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(C) Magnitude of the coefficients $Z_N^j$ for the forward-first zigzag schemes (c.f. equations (20)–(25)).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(D) Magnitude of the coefficients $Z_N^j$ for the backward-first zigzag scheme (c.f. equations (20)–(25)).</figcaption></figure>

FIGURE 14. Magnitude of the coefficients used by the centred, one-sided and zigzag schemes.

## APPENDIX A. MAGNITUDE OF THE FINITE DIFFERENCE COEFFICIENTS

To make it easier to compare the values of the various coefficients involved, we provide four figures (Fig. 14a to 14d) illustrating the magnitude of the finite difference coefficients for the forward, centred and zigzag schemes.
## APPENDIX B. NUMERICAL DIFFUSION OF THE ZIGZAG SCHEMES

The centred, upwind and zigzag finite differences schemes all introduce some numerical diffusion. If we compare the rate at which these schemes dissipate, we notice that the energy loss from the zigzag schemes is lower than the one-sided schemes and greater than the centred schemes. This middle-ground is not surprising since the zigzag schemes can be seen as a mix of one-sided and centred schemes. To illustrate this observation with a concrete example, we consider the solution of the transport equation (34) using a Runge-Kutta scheme of order 3 and spatial schemes of order 2. We solve this equation for $(t, x) \in [0, 15] \times [-20, 20]$ using a periodised error function as the initial condition, with the parameters $\Delta t = 0.05$ and $\Delta x = 0.01$. For the centred, zigzag and upwind schemes we obtain an energy loss of respectively $2.18 \cdot 10^{-9}$, $3.71 \cdot 10^{-8}$ and $1.07 \cdot 10^{-7}$.
## APPENDIX C. STABILITY REGIONS FOR THE FORWARD EULER METHOD

We provide Fig. 15 to 19 the stability regions in the $\lambda - \kappa$ plane ($\lambda$ and $\kappa$ respectively being the stability-number and the scaled wavenumber) for the forward Euler (FE) method in time, for spatial orders $N = 1$ to 5, respectively (left to right) for centred (c), centred staggered (cs), forward (f), zigzag forward-first (zff) and zigzag forward-first staggered (zffs) schemes. The greyed out areas represent the couples $(\lambda, \kappa)$ for which the scheme is stable (i.e. the amplification factor $|G| \le 1$), while the dotted contours correspond to the critical case $|G| = 1$. A given scheme is conditionally stable if and only if there exist a $\lambda_c \in \mathbb{R}$ such that the line of equation $\lambda = \lambda_c$ is included in the grey area.

<figure><img src="image_11.png" /><figcaption>FIGURE 15. FE; left to right : stability regions for **Order 1** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 16. FE; left to right : stability regions for **Order 2** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 17. FE; left to right : stability regions for **Order 3** c, cs, f, zff and zffs.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>FIGURE 18. FE; left to right : stability regions for **Order 4** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIGURE 19. FE; left to right : stability regions for **Order 5** c, cs, f, zff and zffs.</figcaption></figure>
## APPENDIX D. STABILITY REGIONS FOR RK2

Fig. 20 to 25 show the stability regions in the $\lambda—\kappa$ plane ($\lambda$ and $\kappa$ respectively being the stability-number and the scaled wavenumber) for RK2 time integrator, for spatial orders $N = 1$ to 6, respectively (left to right) for centred (c), centred staggered (cs), forward (f), zigzag forward-first (zff) and zigzag forward-first staggered (zffs) schemes. The greyed out areas represent the couples $(\lambda, \kappa)$ for which the scheme is stable (i.e. the amplification factor $|G| \le 1$), while the dotted contours correspond to the critical case $|G| = 1$. A given scheme is conditionally stable if and only if there exist a $\lambda_c \in \mathbb{R}$ such that the line of equation $\lambda = \lambda_c$ is included in the grey area.

<figure><img src="image_16.png" /><figcaption>FIGURE 20. RK2; left to right : stability regions for **Order 1** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIGURE 21. RK2; left to right : stability regions for **Order 2** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIGURE 22. RK2; left to right : stability regions for **Order 3** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>FIGURE 23. RK2; left to right : stability regions for **Order 4** c, cs, f, zff and zffs.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>FIGURE 24. RK2; left to right : stability regions for **Order 5** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>FIGURE 25. RK2; left to right : stability regions for **Order 6** c, cs, f, zff and zffs.</figcaption></figure>
## APPENDIX E. STABILITY REGIONS FOR RK3

Fig. 26 to 31 show the stability regions in the $\lambda-\kappa$ plane ($\lambda$ and $\kappa$ respectively being the stability-number and the scaled wavenumber) for RK3 time integrator, for spatial orders $N=1$ to $6$, respectively (left to right) for centred (c), centred staggered (cs), forward (f), zigzag forward-first (zff) and zigzag forward-first staggered (zffs) schemes. The greyed out areas represent the couples $(\lambda, \kappa)$ for which the scheme is stable (i.e. the amplification factor $|G| \le 1$), while the dotted contours correspond to the critical case $|G| = 1$. A given scheme is conditionally stable if and only if there exist a $\lambda_c \in \mathbb{R}$ such that the line of equation $\lambda = \lambda_c$ is included in the grey area.

<figure><img src="image_22.png" /><figcaption>FIGURE 26. RK3; left to right : stability regions for **Order 1** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>FIGURE 27. RK3; left to right : stability regions for **Order 2** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>FIGURE 28. RK3; left to right : stability regions for **Order 3** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>FIGURE 29. RK3; left to right : stability regions for **Order 4** c, cs, f, zff and zffs.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>FIGURE 30. RK3; left to right : stability regions for **Order 5** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>FIGURE 31. RK3; left to right : stability regions for **Order 6** c, cs, f, zff and zffs.</figcaption></figure>
APPENDIX F. STABILITY REGIONS FOR RK4

Fig. 32 to 37 show the stability regions in the $\lambda-\kappa$ plane ($\lambda$ and $\kappa$ respectively being the stability-number and the scaled wavenumber) for RK4 time integrator, for spatial orders $N=1$ to $6$, respectively (left to right) for centred (c), centred staggered (cs), forward (f), zigzag forward-first (zff) and zigzag forward-first staggered (zffs) schemes. The greyed out areas represent the couples $(\lambda, \kappa)$ for which the scheme is stable (i.e. the amplification factor $|G| \le 1$), while the dotted contours correspond to the critical case $|G| = 1$. A given scheme is conditionally stable if and only if there exist a $\lambda_c \in \mathbb{R}$ such that the line of equation $\lambda = \lambda_c$ is included in the grey area.

<figure><img src="image_28.png" /><figcaption>FIGURE 32. RK4; left to right : stability regions for **Order 1** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>FIGURE 33. RK4; left to right : stability regions for **Order 2** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>FIGURE 34. RK4; left to right : stability regions for **Order 3** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>FIGURE 35. RK4; left to right : stability regions for **Order 4** c, cs, f, zff and zffs.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>FIGURE 36. RK4; left to right : stability regions for **Order 5** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>FIGURE 37. RK4; left to right : stability regions for **Order 6** c, cs, f, zff and zffs.</figcaption></figure>
## APPENDIX G. STABILITY REGIONS FOR RK5

Fig. 38 to 43 show the stability regions in the $\lambda-\kappa$ plane ($\lambda$ and $\kappa$ respectively being the stability-number and the scaled wavenumber) for RK5 time integrator, for spatial orders $N=1$ to 6, respectively (left to right) for centred (c), centred staggered (cs), forward (f), zigzag forward-first (zff) and zigzag forward-first staggered (zffs) schemes. The greyed out areas represent the couples $(\lambda, \kappa)$ for which the scheme is stable (i.e. the amplification factor $|G| \le 1$), while the dotted contours correspond to the critical case $|G| = 1$. A given scheme is conditionally stable if and only if there exist a $\lambda_c \in \mathbb{R}$ such that the line of equation $\lambda = \lambda_c$ is included in the grey area.

<figure><img src="image_34.png" /><figcaption>FIGURE 38. RK5; left to right : stability regions for **Order 1** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>FIGURE 39. RK5; left to right : stability regions for **Order 2** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>FIGURE 40. RK5; left to right : stability regions for **Order 3** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>FIGURE 41. RK5; left to right : stability regions for **Order 4** c, cs, f, zff and zffs.</figcaption></figure>
<figure><img src="image_38.png" /><figcaption>FIGURE 42. RK5; left to right : stability regions for **Order 5** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>FIGURE 43. RK5; left to right : stability regions for **Order 6** c, cs, f, zff and zffs.</figcaption></figure>
## APPENDIX H. STABILITY REGIONS FOR RK6

Fig. 44 to 49 show the stability regions in the $\lambda—\kappa$ plane ($\lambda$ and $\kappa$ respectively being the stability-number and the scaled wavenumber) for RK6 time integrator, for spatial orders $N = 1$ to 6, respectively (left to right) for centred (c), centred staggered (cs), forward (f), zigzag forward-first (zff) and zigzag forward-first staggered (zffs) schemes. The greyed out areas represent the couples $(\lambda, \kappa)$ for which the scheme is stable (i.e. the amplification factor $|G| \le 1$), while the dotted contours correspond to the critical case $|G| = 1$. A given scheme is conditionally stable if and only if there exist a $\lambda_c \in \mathbb{R}$ such that the line of equation $\lambda = \lambda_c$ is included in the grey area.

<figure><img src="image_40.png" /><figcaption>FIGURE 44. RK6; left to right : stability regions for **Order 1** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>FIGURE 45. RK6; left to right : stability regions for **Order 2** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>FIGURE 46. RK6; left to right : stability regions for **Order 3** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>FIGURE 47. RK6; left to right : stability regions for **Order 4** c, cs, f, zff and zffs.</figcaption></figure>
<figure><img src="image_44.png" /><figcaption>FIGURE 48. RK6; left to right : stability regions for **Order 5** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>FIGURE 49. RK6; left to right : stability regions for **Order 6** c, cs, f, zff and zffs.</figcaption></figure>
# APPENDIX I. STABILITY REGIONS FOR RK7

Fig. 50 to 55 show the stability regions in the $\lambda-\kappa$ plane ($\lambda$ and $\kappa$ respectively being the stability-number and the scaled wavenumber) for RK7 time integrator, for spatial orders $N=1$ to 6, respectively (left to right) for centred (c), centred staggered (cs), forward (f), zigzag forward-first (zff) and zigzag forward-first staggered (zffs) schemes. The greyed out areas represent the couples $(\lambda, \kappa)$ for which the scheme is stable (i.e. the amplification factor $|G| \le 1$), while the dotted contours correspond to the critical case $|G| = 1$. A given scheme is conditionally stable if and only if there exist a $\lambda_c \in \mathbb{R}$ such that the line of equation $\lambda = \lambda_c$ is included in the grey area.

<figure><img src="image_46.png" /><figcaption>FIGURE 50. RK7; left to right : stability regions for **Order 1** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_47.png" /><figcaption>FIGURE 51. RK7; left to right : stability regions for **Order 2** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_48.png" /><figcaption>FIGURE 52. RK7; left to right : stability regions for **Order 3** c, cs, f, zff and zffs.</figcaption></figure>

<figure><img src="image_49.png" /><figcaption>FIGURE 53. RK7; left to right : stability regions for **Order 4** c, cs, f, zff and zffs.</figcaption></figure>
<figure><img src="image_50.png" /></figure>

FIGURE 54. RK7; left to right : stability regions for **Order 5** c, cs, f, zff and zffs.

<figure><img src="image_51.png" /></figure>

FIGURE 55. RK7; left to right : stability regions for **Order 6** c, cs, f, zff and zffs.