Compact Approximate Taylor methods for systems of conservation laws

TL;DR

This paper introduces the Compact Approximate Taylor (CAT) methods, high-order schemes for conservation laws that are stable, accurate, and capable of handling discontinuities with shock-capturing techniques, showing promising results in various test cases.

Contribution

The paper develops a new family of high-order, linearly stable CAT methods based on centered stencils, extending Lax-Wendroff methods and incorporating shock-capturing techniques.

Findings

CAT methods achieve order 2p accuracy with CFL<1 stability.

WENO-CAT methods effectively prevent oscillations near discontinuities.

CAT methods perform well on test cases like Burgers and Euler equations.

Abstract

A new family of high order methods for systems of conservation laws are introduced: the Compact Approximate Taylor (CAT) methods. These methods are based on centered (2p + 1)-point stencils where p is an arbitrary integer. We prove that the order of accuracy is 2p and that CAT methods are an extension of high-order Lax-Wendroff methods for linear problems. Due to this, they are linearly L2-stable under a CFL<1 condition. In order to prevent the spurious oscillations that appear close to discontinuities two shock-capturing techniques have been considered: a fux-limiter technique (FL-CAT methods) and WENO reconstruction for the frst time derivative (WENO-CAT methods). We follow the WENO-Lax Wendroff Approximate Taylor method of Zorio, Baeza and Mullet (2017) in the second approach. A number of test cases are considered to compare these methods with other WENO-based schemes: the linear transport equation, Burgers equation, and the 1D compressible Euler system are considered. Although CAT methods present an extra computational cost due to the local character, this extra cost is compensated by the fact that they still give good solutions with CFL values close to 1.

Figures

Tables

Table 1. Table 1: Order of the formula ( 3 ).

	1	2	3	4
1	2	4	6	8
2	2	4	6	8
3		2	4	6
4		2	4	6
5			2	4
6			2	4
7				2
8				2

Table 2. Table 2: Average rate time and flops to increase from LW-CAT2 to LW-CAT2p, with p = 2 , 3 , 4 , 5 𝑝 2 3 4 5 p=2,3,4,5 for the scalar transport equation with initial conditions ( 62 ) and t = 1 𝑡 1 t=1 s.

	LW-CAT2		LW-CAT4		LW-CAT6		LW-CAT8		LW-CAT10
time	1		2.98		7.72		18.87		42.66
flops	1		1.61		2.51		3.69		5.16

Table 3. Table 3: Linear transport equation with initial condition ( 63 ), C ​ F ​ L = 0.5 𝐶 𝐹 𝐿 0.5 CFL=0.5 and t = 1 𝑡 1 t=1 s: L 1 superscript 𝐿 1 L^{1} errors and accuracy order for LW-CAT2p, p = 1 , 2 , 3 𝑝 1 2 3 p=1,2,3 .

	LW-CAT2		LW-CAT4		LW-CAT6
$\Delta x$	Error $\parallel \cdot {\parallel }_1$	Order $\parallel \cdot {\parallel }_1$	Error $\parallel \cdot {\parallel }_1$	Order $\parallel \cdot {\parallel }_1$	Error $\parallel \cdot {\parallel }_1$	Order $\parallel \cdot {\parallel }_1$
0.1053	3.68e-02		1.40e-02		7.88e-03
0.0526	6.84e-03	2.43	3.50e-05	8.64	4.25e-08	7.50
0.0263	1.70e-03	2.00	2.19e-06	4.00	6.49e-10	6.03
0.0132	4.27e-04	2.00	1.36e-07	4.00	9.89e-12	6.04
0.0066	1.06e-04	2.00	8.55e-09	4.00	1.53e-13	6.01
0.0033	2.66e-05	2.00	5.34e-10	4.00	2.64e-15	5.96

Table 4. Table 4: Linear transport equation with initial condition ( 63 ), C ​ F ​ L = 0.5 𝐶 𝐹 𝐿 0.5 CFL=0.5 and t = 1 𝑡 1 t=1 s: L 1 superscript 𝐿 1 L^{1} errors and accuracy order for WENO5-RK3 and WENO5-LWA5.

	WENO5-RK3		WENO-LWA5
$\Delta x$	Error $\parallel \cdot {\parallel }_1$	Order $\parallel \cdot {\parallel }_1$	Error $\parallel \cdot {\parallel }_1$	Order $\parallel \cdot {\parallel }_1$
0.1053	2.03e-03		5.44e-05
0.0526	6.06e-05	5.06	1.65e-06	5.04
0.0263	1.87e-06	5.02	5.04e-08	5.04
0.0132	5.83e-08	5.00	1.51e-09	5.05
0.0066	1.82e-09	5.00	4.41e-11	5.10
0.0033	5.65e-11	5.01	1.15e-12	5.25

Table 5. Table 5: Burgers equation with initial condition ( 62 ), C ​ F ​ L = 0.5 𝐶 𝐹 𝐿 0.5 CFL=0.5 and t = 2 𝑡 2 t=2 s: computational time in seconds

Nodes	LW-CAT4	FL-CAT4		WENO5-CAT4		WENO5-LWA5
80	0.469	0.4702		0.5508		0.1998
160	0.8090	0.9668		0.9700		0.4616
320	1.9838	1.9908		2.0064		1.0346

Table 6. Table 6: Burgers equation with initial condition ( 62 ), C ​ F ​ L = 0.9 𝐶 𝐹 𝐿 0.9 CFL=0.9 and t = 2 𝑡 2 t=2 s: computational time in seconds

Nodes	FL-CAT4		WENO5-CAT4
80	0.2003		0.2078
160	0.4024		0.5001
320	1.0494		1.045

Table 7. Table 7: Burgers equation with initial condition ( 63 ), C ​ F ​ L = 0.5 𝐶 𝐹 𝐿 0.5 CFL=0.5 and t = 0.5 𝑡 0.5 t=0.5 s: L 1 superscript 𝐿 1 L^{1} errors and accuracy order for LW-CAT2p, p = 1 ,  2 ,  3 𝑝 1 2 3 p=1,\,2,\,3 .

	LW-CAT2		LW-CAT4		LW-CAT6
$\Delta x$	Error $\parallel \cdot {\parallel }_1$	Order $\parallel \cdot {\parallel }_1$	Error $\parallel \cdot {\parallel }_1$	Order $\parallel \cdot {\parallel }_1$	Error $\parallel \cdot {\parallel }_1$	Order $\parallel \cdot {\parallel }_1$
0.1053	7.94e-03		9.01e-04		2.09e-04
0.0526	2.08e-03	1.93	6.13e-05	3.88	4.27e-06	5.62
0.0263	5.22e-04	1.99	3.89e-06	3.98	7.49e-08	5.83
0.0132	1.29e-04	2.01	2.44e-07	4.00	1.20e-09	5.96
0.0066	3.08e-05	2.00	1.51e-08	4.00	1.87e-11	6.00
0.0033	6.16e-06	2.00	8.76e-10	4.00	2.84e-13	6.00

Table 8. Table 8: 1D Euler equations with initial condition ( 77 ), C ​ F ​ L = 0.5 𝐶 𝐹 𝐿 0.5 CFL=0.5 and t = 0.5 𝑡 0.5 t=0.5 s: L 1 superscript 𝐿 1 L^{1} errors and accuracy order for LW-CAT2p, p = 1 ,  2 ,  3 𝑝 1 2 3 p=1,\,2,\,3 .

	LW-CAT2		LW-CAT4		LW-CAT6
$\Delta x$	Error $\parallel \cdot {\parallel }_1$	Order $\parallel \cdot {\parallel }_1$	Error $\parallel \cdot {\parallel }_1$	Order $\parallel \cdot {\parallel }_1$	Error $\parallel \cdot {\parallel }_1$	Order $\parallel \cdot {\parallel }_1$
0.1053	3.34e-03		8.57e-04		5.49e-04
0.0526	8.82e-03	1.92	9.93e-05	3.11	3.53e-05	4.96
0.0263	2.28e-04	1.95	7.31e-06	3.76	1.01e-06	5.12
0.0132	5.69e-05	2.01	4.81e-07	3.93	1.94e-08	5.71
0.0066	1.35e-05	2.07	3.02e-08	3.99	3.21e-10	5.92
0.0033	2.71e-06	2.30	1.78e-09	4.08	4.99e-12	6.01

Equations

$$u_t+a{u}_x=0.$$

$$u_i^{n+1}=u_i^n+\sum _{k=1}^m\frac{(-1{)}^k{c}^k}{k!}\sum _{j=-p}^p{\delta }_{p,j}^k{u}_{i+j}^n,$$

$$f^{(k)}(x_i)\simeq D_{p,i}^k(f,\Delta x)=\frac{1}{\Delta x^k}\sum _{j=-p}^p{\delta }_{p,j}^{k}f(x_{i+j}),$$

$${\partial }_t^{k}u=(-1{)}^k{a}^k{\partial }_x^{k}u,k=1,2\dots$$

$$f^{(k)}(x_i+q\Delta x)\simeq A_{p,i}^{k,q}(f,\Delta x)=\frac{1}{\Delta x^k}\sum _{j=-p+1}^p{\gamma }_{p,j}^{k,q}f(x_{i+j}),$$

$$D_{p,i}^k(w_{\ast },\Delta x)$$

$$A_{p,i}^{k,q}(w_{\ast },\Delta x)$$

$${\partial }_x^{k}u(x_i+q\Delta x,t_n)$$

$${\partial }_t^{k}u(x_i,t_n+q\Delta t)$$

$$u_i^{n+1}=u_i^n+\sum _{k=1}^m\frac{(-1{)}^k{a}^k\Delta t^k}{k!}{D}_{p,i}^k(u_{\ast }^n,\Delta x).$$

$$f^{(k)}(0)\simeq D_{p,0}^k(f,1)=\sum _{j=-p}^p{\delta }_{p,j}^{k}f(j),$$

$$f^{(k)}(q)\simeq A_{p,0}^{k,q}(f,1)=\sum _{j=-p+1}^p{\gamma }_{p,j}^{k,q}f(j).$$

$$\sum _{j=-p}^p{j}^k{\delta }_{p,j}^k=k!,\sum _{j=-p}^p{j}^s{\delta }_{p,j}^k=0,s\ne k,0\le s,k\le 2p.$$

$$\sum _{j=-p+1}^p{j}^k{\gamma }_{p,j}^{k,0}=k!,\sum _{j=-p+1}^p{j}^s{\gamma }_{p,j}^{k,0}=0,s\ne k,0\le s,k\le 2p-1.$$

$$\sum_{j=-p+1}^{p}\gamma^{k,q}_{p,j}=\left\{\begin{array}[]{ll}1&\text{if $k=1$,}\\ 0&\text{otherwise.}\end{array}\right.$$

$$F_{p,j}(x)=\prod _{r=-p,r\ne j}^p\frac{(x-r)}{(j-r)},-p\le j\le p,$$

$${\delta }_{p,j}^k=F_{p,j}^{(k)}(0),$$

$$F_{p,j}(x)=\sum _{k=0}^{2p}\frac{{\delta }_{p,j}^k}{k!}{x}^k.$$

$${\delta }_{p,j}^k=(-1{)}^k{\delta }_{p,-j}^k;$$

$${\delta }_{p,0}^k=0\text{ if k is odd};$$

$$\sum _{j=-p}^p{\delta }_{p,j}^k{j}^{(2p+1)}=0\text{ if k is even};$$

$$\sum _{j=-p}^p{\delta }_{p,j}^k{j}^{(2p+2)}=0\text{ if k is odd}.$$

$$F_{p,-j}(x)=F_{p,j}(-x);$$

$${\delta }_{p,p}^k={\gamma }_{p,p}^{k-1,1/2};$$

$${\delta }_{p,j}^k={\gamma }_{p,j}^{k-1,1/2}-{\gamma }_{p,j+1}^{k-1,1/2},j=-p+1,\dots ,p-1;$$

$${\delta }_{p,-p}^k=-{\gamma }_{p,-p+1}^{k-1,1/2}.$$

$$f^{(k-1)}(1/2)\simeq A_{p,0}^{k-1,1/2}(f,1)=\sum _{j=-p+1}^p{\gamma }_{p,j}^{k-1}f(j),$$

$$f^{(k-1)}(-1/2)\simeq A_{p,-1}^{k-1,1/2}(f,1)=\sum _{j=-p+1}^p{\gamma }_{p,j}^{k-1}f(j-1),$$

$$f^{(k)}(0)\simeq A_{p,0}^{k-1,1/2}(f,1)-A_{p,-1}^{k-1,1/2}(f,1).$$

$$A_{p,0}^{k-1,1/2}(f,1)-A_{p,-1}^{k-1,1/2}(f,1)=f^{(k-1)}(1/2)-f^{(k-1)}(-1/2)=f^{(k)}(0),$$

$$g^{\prime }(0)\simeq g(1/2)-g(-1/2),$$

$$f^{(k)}(0)\simeq$$

$$+({\gamma }_{p,-p+1}^{k-1,1/2}-{\gamma }_{p,-p+2}^{k-1,1/2})f(-p+1)-{\gamma }_{p,-p+1}^{k-1,1/2}f(-p),$$

$$\sum _{j=-p+1}^p{\gamma }_{p,j}^{s,q}{\gamma }_{p,l}^{k-s,j}={\gamma }_{p,l}^{k,q},l=-p+1,\dots ,p.$$

Full text

1 Introduction

Lax-Wendroff methods for linear systems of conservation laws are based on Taylor expansions in time in which the time derivatives are transformed into spatial derivatives using the equations [5], [6],[12], [15]. The spatial derivatives are then discretized by means of centered high-order differentiation formulas. This procedure allows to derive numerical methods of order $2p$, where $p$ is an arbitrary integer, using a centered $(2p+1)$-point stencil that are $L^{2}$-stable (a review on these methods will be presented in Section 2).

This paper focuses on the extension of Lax-Wendroff methods to nonlinear systems of conservation laws. This problem is closely related to the design of numerical schemes based on the methods of lines in which the time discretization is performed by means of Taylor or Approximate Taylor methods. Many authors have focused on the design of this type of methods that can be an alternative to time discretizations based on SSP Runge-Kutta methods [9] that lead to some stability restrictions for orders bigger than three. The main difficulty to extend Lax-Wendroff methods to nonlinear problems come from the transformation of time derivatives into spatial derivatives using the equations. A first strategy to do this is given by the Cauchy-Kovalevskaya (CK) procedure. In [9] this procedure has been used together with WENO reconstructions for the spatial discretization. The main benefit compared to RK time discretizations is that only one WENO reconstruction is needed at each spatial cell per time step. On the other hand, the main drawback comes from the fact that the CK procedure leads to expressions whose number of terms grow exponentially what implies high computational costs and difficult implementations. In [14] an alternative has been presented based on an Approximate Taylor (AT) method in which the time derivatives are approximated using high-order centered differentiation formulas combined with Taylor approximations in time that are computed in a recursive way. The resulting method is easy to implement and shows a good performance. Nevertheless AT schemes are not proper generalizations of Lax-Wendroff methods: they have $(4p+1)$-point stencils and worse linear stability properties than the original Lax-Wendroff methods. Nevertheless, they can be stabilized by using one WENO reconstruction per spatial cell and time step, like in [9] the resulting methods give good results under a $CFL-0.5$ condition typically.

In order to design numerical methods that are proper generalization of Lax-Wendroff methods, a compact variant of the AT procedure is introduced here: first, the conservative expression of the high-order Lax-Wendoff methods is considered; then, the derivatives appearing in the expression of the numerical flux are computed using $2p$-point differentiation formulas. This strategy lead to Compact Approximate Taylor (CAT) methods that have $2p+1$-point stencil and order of accuracy $2p$. They reduce to the Lax-Wendroff method when applied to a linear systems and thus they are linearly $L^{2}$-stable under a $CFL-1$ condition. Nevertheless, the number of operations to perform a step is bigger than the original AT methods: this is due to the fact that the approximation of the time derivatives is local in the sense that they depend both on the point and on the stencil. However, this extra cost is compensated by better stability properties. As it happens with its linear counterpart, CAT methods lead to spurious oscillations near discontinuities. In order to cure them, two shock-capturing techniques are considered here: a flux-limiter technique and the use of a WENO reconstruction per cell and time step, as in [14] and [9].

This paper is organized as follows: in Section 2 a review of high order Lax-Wendroff methods for the linear transport equation is presented, including the study of the order and the $L^{2}$-stability as well as the computation and properties of the coefficients. Section 3 is devoted to their extension to nonlinear problems: first the AT technique is recalled and then CAT methods are presented. We show that they reduce to Lax-Wendroff methods when applied to a linear problem and we analyze the order of accuracy. In Section 4 the techniques considered to cure the spurious oscillations near the discontinuities are presented. In Section 5 CAT methods are compared in a number of test cases with WENO-RK methods and AT methods. The linear transport equation, Burgers equation, and the 1D compressible Euler system are considered. Future developments and conclusions are drawn in Section 6.

2 The high-order Lax-Wendroff method for linear problems

Let us first consider the linear scalar equation:

[TABLE]

We consider the numerical method:

[TABLE]

where $\{x_{i}\}$ are the nodes of a uniform mesh of step $\Delta x$; $u_{i}^{n}$ is an approximation of the point value of the solution at $x_{i}$ at the time $n\Delta t$, where $\Delta t$ is the time step; $p\geq 1$ is a natural number; $c={a\Delta t}/{\Delta x}$; and $\delta^{k}_{p,j}$ are the coefficients of the centered interpolatory formula of numerical differentiation based on a $(2p+1)$-point stencil:

[TABLE]

where $f^{(k)}$ represents the $k$-th derivative of a one-variable function $f$ and $f^{(0)}=f$. The expression of the numerical method is obtained by applying a Taylor expansion in time, and replacing time derivatives by space derivatives through the identities

[TABLE]

2.1 Formulas of numerical differentiation

Besides (3) the following family of interpolatory formulas based on a $2p$-point stencil will be used in this work:

[TABLE]

i.e. $A^{k,q}_{p,i}(f,\Delta x)$ is the numerical differentiation formula that approximates the $k$-th derivative at the point $x_{i}+q\Delta x$ using the values of the function at the $2p$ points $x_{i-p+1},\dots,x_{i+p}$. Observe that the coefficients, like in (3), do not depend on $i$.

Given a variable $w$, the following notation will be used:

[TABLE]

to indicate that the formulas are applied to the approximations of $w$, $w_{i}$, and not to its exact point values $w(x_{i})$. In cases where there are two or more indexes, the symbol ${*}$ will be used to indicate with respect to which the differentiation is applied. For instance:

[TABLE]

Using this notation, the algorithm (2) writes as follows:

[TABLE]

Let us discuss some properties of the coefficients of the numerical differentiation formulas (3) and (5) and some relations between them that will be used in that follows. Since the coefficients are independent of $\Delta x$ and $i$, we can consider, without loss of generality, the case $i=0$, $x_{0}=0$, $\Delta x=1$:

[TABLE]

[TABLE]

Since (7) is exact for polynomials of degree $\leq 2p$, by applying the formula to $x^{s}$, $s=0,\dots,2p$ at $x=0$, we get that the coefficients have to satisfy the equalities

[TABLE]

Analogously:

[TABLE]

[TABLE]

As it is well known, the coefficients $\delta^{k}_{p,j}$ are related to the Lagrange basis polynomials

[TABLE]

through the equalities:

[TABLE]

which allow us to write the Taylor expansion of $F_{p,j}$ centered at $x=0$ as follows:

[TABLE]

Proposition 1**.**

The coefficients $\delta^{k}_{p,j}$ of the formula (7), satisfy:

[TABLE]

• Proof.

(15) is deduced from the equality

[TABLE]

. Using (15) we get (16). (17) and (18) are deduced from (15) and (16).

$\Box$

Proposition 2**.**

For $k\geq 1$ the following relations hold:

[TABLE]

• Proof.

Let us consider the formulas

[TABLE]

[TABLE]

that are exact for polynomials of degree $\leq 2p-1$. Let us consider now the formula

[TABLE]

If $f$ is a polynomial of degree $2p$, then (22) and (23) are exact for $f$, furthermore

[TABLE]

where we have used that the formula

[TABLE]

is exact for polynomials of degree 1. Therefore, (24) coincide with (7). The proof is finished by writing (24) in the form

[TABLE]

and matching the weights.

$\Box$

Proposition 3**.**

Given $1\leq k\leq 2p-1$, $0\leq s\leq k$:

[TABLE]

• Proof.

The proof is similar to the one of the preceding in Proposition 2: consider the formula

[TABLE]

with

[TABLE]

check that it is exact for polynomials of degree $2p-1$; write it in the form:

[TABLE]

and match its weights with those of (8).

$\Box$

2.2 Conservative form of (2)

From the proof of Proposition 2 we deduce an alternative form for (3):

[TABLE]

Using this form in (6), the numerical method (2) can be written as:

[TABLE]

with

[TABLE]

Using (11) it is straightforward to verify that $F^{p}_{i+1/2}$ is a consistent numerical flux, what proves that (2) is a conservative method.

2.3 Computation of the coefficients: an iterative algorithm

Notice that (9) constitutes a $(2p+1)\times(2p+1)$ linear system with a Vandermonde matrix that can be used to compute $\delta^{k}_{p,i}$ . Nevertheless, as it is well-known, this system is ill-conditioned, so that it is recommendable to compute them by using an alternative algorithm: we adapt the recursive algorithm proposed in [3]. The following notation is adopted:

[TABLE]

Let us derive some recurrence formulas to compute the coefficients:

$\delta^{k}_{p,j}$ for $j=0,..,p-1$.

From (12), we obtain

[TABLE]

Using then the Taylor expansions (14) in (29) we get

[TABLE] 2. 2.

$\delta^{k}_{p,j}$ with $j=p$.

Substituting j=p in (12), we get

[TABLE]

and, using (14), we obtain:

[TABLE] 3. 3.

$\delta^{k}_{p,j}$ for ${j=-p,\dots,-1}$. $\eqref{p1}$ is used.

The algorithm is computed only once in increasing order of $p$. The coefficients $\gamma^{k,q}_{p,j}$ are computed using the algorithms described in [3],[2] and $\gamma^{k,1/2}_{p,j}$ is obtained from $\delta^{k+1}_{p,j}$.

2.4 Order of accuracy

Proposition 4**.**

The formula of numerical differentiation (3) has order of accuracy $\alpha_{k}-k$, with,

[TABLE]

• Proof.

Let $f$ be a function of class $C^{\alpha_{k}+1}$. Applying Taylor expansions and properties (9) and (17), we obtain:

[TABLE]

where

[TABLE]

$\Box$

Table 1 shows the order of (3) for different values of $p$ and $k$.

Proposition 5**.**

*The discretization error of the numerical method (2) is of order $O(\Delta t^{m+1}+\Delta x^{2p+1})$. *

• Proof.

Let $u$ be a smooth enough solution of (1). Using Proposition 4 we obtain

[TABLE]

where (4) has been used.

$\Box$

As a consequence, the order of accuracy of (2) is $\min(m,2p)$. Therefore, the optimal combination of these parameters is $m=2p$. From now on, we shall assume that this relation holds.

2.5 Modified equation and stability

Taking into account that $m=2p$ and (18), the local discretization error is as follows:

[TABLE]

Using (34) and (14) we get:

[TABLE]

where $q(c)$ is the polynomial of degree $\leq 2p$ that interpolates the points

[TABLE]

Since $q$ is clearly an odd function, we finally obtain:

[TABLE]

Reasoning in a similar way, we get:

[TABLE]

where $r$ is the polynomial of degree $\leq 2p$ that interpolates the points

[TABLE]

Using now (4), (35), and (36), the local discretization error can be written as follows:

[TABLE]

with

[TABLE]

Therefore, the numerical method solves with order $\mathcal{O}(\Delta^{2p+2})$ the following modified equation

[TABLE]

where

[TABLE]

In order to study the stability, let us look for an elementary solution $u(x,t)$ of (40) of the form

[TABLE]

where $\alpha$ is complex number. The following equality has to be satisfied:

[TABLE]

Therefore:

[TABLE]

The numerical method is thus expected to be stable if the real part is negative, i.e.

[TABLE]

or, equivalently

[TABLE]

$h_{2}$ is a pair polynomial of degree $2p+2$ such that

[TABLE]

Moreover, [math] is a double root of $h_{2}$ and $\pm 1,\dots,\pm p$ are single roots. Analyzing the change of signs of $h_{2}$, it can be easily checked that

[TABLE]

and thus (41) is satisfied if $c\in[0,1]$ (see Figure 1). This argument shows that the stability of the method under the standard CFL condition $c\leq 1$ can be formalized: extended details and a formal proof can be found in [8]. A study of the modified equation for the Lax Wendroff method can also be found in [13].

3 Extension to nonlinear problems

3.1 Approximate Taylor method

Following [14], instead of using the Cauchy-Kovaleskaya process to extend (27)-(28) to nonlinear problems

[TABLE]

we use the equalities

[TABLE]

To derive the expression of the numerical method, let us suppose that approximations

[TABLE]

are available. Then,

[TABLE]

being

[TABLE]

where, $\lceil\cdot\rceil$ denotes the ceiling function.

Using these approximations to approximate the Taylor expansion, we obtain the method

[TABLE]

Equivalently, using (26), the numerical method can be written in conservative form (27) with numerical flux

[TABLE]

being

[TABLE]

Now, to compute the approximations $\tilde{f}^{(k-1)}_{i}$, new Taylor expansions in time are used recursively as follows:

Compute $\tilde{f}^{(0)}_{i}=f(u^{n}_{i})$.

• (a)

For $k=1\dots 2p$:

• i.

Compute

[TABLE]

• ii.

Compute

[TABLE]

• iii.

Compute

[TABLE]

where

[TABLE]

Observe that Taylor expansions are used to approximate $f(u(x_{i},t_{n}+r\Delta t))$ and once these approximations have been computed, the centered formula of numerical differentiation (3) is used to approximate the temporal derivatives.

This method is order $2p$, but it is not a generalization of (2) in the sense that this latter method is not recovered if $f(u)=au$. To see this, consider $p=1$ and $f(u)=au$: it can be easily checked that (45) writes as follows

[TABLE]

which is different from the standard Lax-Wendroff method: (45) is a $(4p+1)$-point method whose stability properties are worse than those of the standard Lax-Wendroff method. (see [6]).

3.2 Compact Approximate Taylor method

In order to prevent the increase of the stencil observed for Approximate Taylor methods, we consider a modification based on the conservative form of the method. The numerical flux $F^{p}_{i+1/2}$ will be computed using only the approximations

[TABLE]

so that the values used to update $u_{i}^{n+1}$ are only those of the centered $(2p+1)$-point stencil, like in the linear case. In fact, we will show that this modification is a proper generalization of the Lax-Wendroff method for linear problems.

In order to be able to compute the numerical fluxes using only (49), for every $i$ we will compute local approximations of

[TABLE]

that will be represented by

[TABLE]

These approximations are local in the sense that $i_{1}+j_{1}=i_{2}+j_{2}$, does not imply that $\tilde{f}^{(k-1)}_{i_{1},j_{1}}=\tilde{f}^{(k-1)}_{i_{2},j_{2}}$. Once these approximations have been computed, the numerical flux is given by

[TABLE]

with

[TABLE]

Now, given $i$, to compute the approximations $\tilde{f}^{(k-1)}_{i,j}$, new Taylor expansions in time are used recursively as follows:

Compute $\tilde{f}^{(0)}_{i,j}=f(u^{n}_{i+j}),\quad j=-p+1,\dots,p.$

• (a)

For $k=1\dots 2p$:

• i.

Compute

[TABLE]

where

[TABLE]

• ii.

Compute

[TABLE]

• iii.

Compute

[TABLE]

with

[TABLE]

Notice that, unlike the Approximate Taylor methods (in which all the derivatives were approximated using the centered $(2p+1)$-point formula), in this algorithm the stencil $x_{i-p+1},\dots,x_{i+p}$ is used for the space derivatives and the stencil $t_{n-p+1},\dots,t_{n+p}$ for the time derivative.

Theorem 1**.**

*The compact approximate Taylor method reduces to (2) when $f(u)=au$. *

• Proof.

For $k>1$ we have:

[TABLE]

where (10) has been used. On the other hand:

[TABLE]

where (25) has been used. By recurrence:

[TABLE]

Next,

[TABLE]

where (25) has been used. Finally,

[TABLE]

what is the numerical flux (28) corresponding to (2), as we wanted to prove.

$\Box$

As a consequence, we obtain that the compact approximate Taylor method is linearly stable (in the $L^{2}$ sense) under the usual CFL condition

[TABLE]

Theorem 2**.**

*The compact approximate Taylor method is order $2p$. *

• Proof.

Let us perform a step of the method starting from the point values at time $t^{n}$, $u(x_{i},t_{n})$, of a smooth enough exact solution. We assume that $\Delta t/\Delta x$ remains constant.

First we have:

[TABLE]

Next:

[TABLE]

where

[TABLE]

is the first order Taylor polynomial in time of $u$ in $(x_{i+j},t_{n})$. Then

[TABLE]

where (10) has been used. This result can be extended by induction to every $k$ between 1 and $2p-1$ as follows:

[TABLE]

Using this equality we get:

[TABLE]

$\Box$

Remark: In the Approximate Taylor method proposed in [14] the derivatives $\tilde{u}^{(k+1)}_{i}$ are computed by applying the $2p_{k}+1$-point centered differentiation formula for first derivatives to $\tilde{f}^{(k)}_{i}$, where $p_{k}$ is given by (44): notice that $p_{k}$ decreases as $k$ increases. The same reduction of the stencil used to compute $\tilde{u}^{(k)}_{i,j}$ could be applied here, what would allow us to reduce the number of computations while preserving the overall order of accuracy. Nevertheless, the resulting method will not be an extension of the linear Lax-Wendroff method. On the other hand, the CPU reduction will be not significant.

3.3 Example: Fourth order compact approximate Taylor method

Since the method is conservative, we will only show in detail how to compute the numerical flux (50) and, to do this, it is enough to specify how to compute

[TABLE]

The procedure is as follows:

$\kappa^{1}_{i+1/2}$: First the assignment

[TABLE]

is done and then:

[TABLE]

• 2.

$\kappa^{2}_{i+1/2}$: The first order time derivatives of $u$ at the nodes $i-1,\dots,i+2$ are approximated by applying the corresponding differentiation numerical formula to $\tilde{f}^{(0)}_{i,j}$:

[TABLE]

Next first order Taylor expansions are used to approximate the values of the flux sixteen space-time local nodes: for $r=-1,\dots,2$

[TABLE]

Then, the first order time derivates of the flux at the nodes $i-1,\dots,i+2$ are approximated by applying the corresponding differentiation numerical formula to $\tilde{f}^{1,r}_{i,j}$:

[TABLE]

Finally;

[TABLE]

• 3.

$\kappa^{3}_{i+1/2}$: the second order time derivatives at the nodes are approximated by

[TABLE]

Second order Taylor expansions are used to compute the fluxes at the sixteen nodes in the space-time mesh: for $r=-1,\dots 2$

[TABLE]

Next, compute

[TABLE]

And finally;

[TABLE]

• 4.

$\kappa^{4}_{i+1/2}$: the third order time derivatives at the nodes are approximated by

[TABLE]

Compute the approximations of the fluxes: for $r=-1,\dots,2$

[TABLE]

Next, compute:

[TABLE]

Finally;

[TABLE]

If $f(u)=au$, then:

[TABLE]

which coincides with the numerical flux of the fourth order Lax-Wendroff in conservative form.

4 Shock-capturing techniques

Although the Compact Approximate Taylor methods are linearly stable in the $L^{2}$ sense under the usual $CFL$-1 condition, they may produce strong oscillations close to a discontinuity of the solution. The goal of this section is to modify the numerical method to avoid these oscillations. Two different techniques are considered here:

4.1 FLUX LIMITER-CAT methods

We consider the numerical method (27) with

[TABLE]

where $F^{L}_{i+1/2}$ is a first order robust numerical flux,  $F^{p}_{i+1/2}$ is given by (50), and $\varphi_{i+1/2}$ is the flux limiter function, see [4] [7], [12]. We consider here

[TABLE]

where $\varphi$ is the van Albada second version flux limiter:

[TABLE]

and

[TABLE]

where, $a_{i+1/2}$ is an estimate of the wave speed.

4.2 WENO-CAT methods

Following [9] WENO reconstructions of the flux are used to stabilize the method. The only differences with the algorithm described in Section 3.2 are the computation of $\tilde{u}^{(1)}_{i,j}$, that is now performed as follows:

[TABLE]

where $\hat{f}_{i+1/2}$ denotes the WENO flux splitting, reconstructions at $x_{i+1/2}$ of the flux function described in [10]. The expression of the numerical flux is then given by:

[TABLE]

4.3 Systems of conservation laws

Although the methods have been described in the scalar case for simplicity, they can be easily applied to systems using vector notation.

5 Numerical Experiments

In this section we apply to some scalar conservation laws and to the 1D Euler equation the following numerical methods:

LW-CAT$q$: Compact Approximate Taylor method of order $q$ (space and time).

• 2.

FL-CAT$q$: Compact Approximate Taylor method of order $q$ with flux limiter technique. The first order methods considered are Lax-Friedrich for scalar problems and HLL for systems.

• 3.

WENO$s$-CAT$q$: Compact Approximate Taylor method of order $q$ with WENO reconstructions of order $s$ to compute $\tilde{u}^{(1)}_{t,i}$.

• 4.

WENO$s$-RK$q$: WENO method of order $s$ for the space discretization and TVD-RK$q$ for the time discretization, see [10].

• 5.

WENO$s$-LWA$q$: Approximate Taylor method of order $q$ with WENO reconstructions of order $s$ to compute $\tilde{u}^{(1)}_{t,i}$, see [14].

5.1 Linear transport equation

We consider first (1) with $a=1$, in the space interval $[0,1]$, with initial condition

[TABLE]

and periodic boundary conditions. A uniform mesh with $N=80$ points is considered and the LW-CAT method (that, in this case, coincides with the Lax-Wendroff method) is applied for $p=1,\dots,5$.

Numerical simulations are shown in Figure 2: the $L^{2}$ stability of the scheme and the appearance of oscillations near the discontinuities can be observed. Next, we apply to the same problem LW-CAT4, FL-CAT4, WENO5-CAT4, WENO5-RK3, and WENO5-LWA5 methods. A general view is shown in Figure 3 together with a zoom of the area of interest. As it can be observed, the results given by WENO5-CAT4, WENO5-RK3 and WENO5-LWA5 are almost identical. Nevertheless, as it will be seen in the next test problem, WENO5-CAT4 still gives good results for CFL close to one, what is not the case for WENO5-RK3 or WENO5-LWA5.

Using the algorithm in Section 3.2 CAT methods are easily extended to any $2p$ order, nevertheless, this involves a significant increase of the CPU time simulation and flops (number of operations required), that should be considered, see Table 2.

Finally, we consider (1) in the space interval $[0,2]$ with initial condition,

[TABLE]

and periodic boundary conditions. Table 3 shows the error and the empirical order for LW-CAT2, LW-CAT4, LW-CAT6, and Table 4   for WENO5-RK3 and WENO5-LWA5 which coincides in all cases with the theoretical one. For smooth solutions WENO-CATq and FL-CATq reduce to the corresponding LW-CATq, so that the accuracy test is not necessary.

5.2 Burgers equation

We consider Burgers equation, i.e. (42) with

[TABLE]

When CAT methods are applied to approximate a discontinuous solution of this nonlinear problem, the oscillations appearing close to the shocks tend to grow and to spoil the numerical solution. Nevertheless, it is still possible to apply these methods by reducing the $CFL$ parameter (the reduction increases with $p$): for instance, Figure 4 shows the results obtained with CAT-LW2$p$, $p=1,2,3,4$ and $CFL={0.8,0.4,0.2,0.1}$, respectively, with initial conditions (62).

Next, the same test problem is solved using LW-CAT4, FL-CAT4, WENO5-RK3 and WENO5-LWA5 methods. Using $CFL=0.5$ we obtain numerical solutions without spurious oscillations for all the methods. Figure 5 shows a general view of solutions and the van Albada flux limiter function on every inter cell. In Figure 6 the results are compared with those obtained with $CFl=0.9$. From Figures 5, 6 we can conclude:

$CFL\leq 0.5$

• (a)

LW-CAT4 show oscillations near the discontinuities, but it is stable.

• (b)

FL-CAT4 is very diffusive near to the discontinuities, due to the selected first order accurate flux limiter function.

• (c)

WENO5-CAT4, WENO5-LWA4 and WENO5-RK3 show good results, stable and essentially the same values.

• 2.

$CFL>0.5$

• (a)

LW-CAT4: the amplitude of oscillations increases near the discontinuities. However, they remain stable.

• (b)

FL-CAT4: conversely to the previous $CFL$ condition, it shows acceptable solutions near the discontinuities.

• (c)

WENO5-CAT4 ,WENO5-LWA5 and WENO5-RK3 : slight oscillations appear near the discontinuities at the beginning of the simulations. Nevertheless, as the time increases, these oscillations tend to dismiss and the result remains acceptable and stable for WENO5-CAT4 while the solutions given by WENO5-LWA5 is very diffusive and the one given by WENO5-RK3 is overdamped.

Although FL-CAT4 shows better results for bigger $CFL$, it fails in smooth regions close to critical points and for systems (as it will be seen in Euler 1D equation). For $CFL\leq 1/2$, WENO5-LWA5 is faster than CAT methods (in the computational cost sense): a comparsion is shown in Table 5. However, we can obtain good solutions for CAT methods using $CFL>0.5$ and obtain similar CPU time computation, see Table 6.

In order to study the order of convergence, we consider again initial condition (63) and periodic boundary conditions. A reference solution at time $t=0.5$s (when the solution is still smooth) is obtained with WENO5-RK3 using a fine grid of $1400$ nodes. The errors and the empirical order are shown in Table 7: the numerical results verify the theoretical analysis.

5.3 1D Euler Equation

We solve the $1D$ Euler equation for gas dynamics

[TABLE]

with

[TABLE]

where $\rho$ is the density, $u$ the velocity, $E$ the total energy per unit volume, and $p$ the pressure. We assume an ideal gas with the equation of state,

[TABLE]

being $\gamma$ the ratio of specific heat capacities of the gas taken as 1.4 and $e$ is the internal energy per unit mass given by:

[TABLE]

We consider the space interval $[0,2]$ with the initial condition:

[TABLE]

and periodic boundary conditions. For this test we take $CFL=0.5$ and $t=0.5$s. We use a fine grid with $1400$-point mesh to compute LW-CAT8 as a reference solution. The results in Table 8 support the theoretically obtained accuracy.

Finally, two tests involving discontinuities are considered:

The Sod Shock tube problem. The initial condition is given by:

[TABLE]

Here, $x\in[-5,5]$, $CFL=0.5$, $t=1$s, and outflow boundary conditions are considered at both sides. For details of this problem see [11]. We compare FL-CAT4, WENO5-CAT4, WENO5-LWA5 and WENO5-RK3 using $450$ points. A reference solution is computed with WENO5-RK3 using a $1400$-point mesh. The flux limiter function is computed for every variable.

While all numerical solutions show stable and similar values over smooth regions (see Figure 7), the quality is different in the interest regions (a,b,c,..,f): an enlarged view of them can be seen in Figure 8. In the figure we can observe that the solution given by FL-CAT4 is the most diffusive one (due to the condition $CFL=0.5$) and that WENO5-LWA5 and WENO5-RK3 give essentially the same results. WENO5-CAT4 achieves some improvements, specially in a,c and d.

The Shu-Osher problem. The initial condition is given by:

[TABLE]

We consider the space interval $x\in[-5,5]$, $CFL=0.5$ and time $t=1$s. For details see [10] test 8. We compare FL-CAT4, WENO5-CAT4, WENO5-LWA5 and WENO5-RK3 using $450$- point mesh and a reference solution computed with WENO5-RK3 using a $1400$-point mesh.

Again, all solutions are closely similar and near to the reference solution with exception of FL-CAT4.

6 Conclusions

In this work, first a review of high order Lax-Wendroff methods for the linear transport equation has been presented, including the study of the order and the $L^{2}$-stability as well as the computation and properties of the coefficients. Next, an extension to nonlinear conservation laws has been introduced with arbitrary even order $2p$ of accuracy, the so-called Compact Approximate Taylor (CAT) methods. Unlike previous applications of Taylor methods to conservation laws, CAT methods have $2p+1$-point centered stencils, like Lax-Wendroff methods for linear problems. Moreover, since they inherit the stability properties of Lax-Wendroff methods, they are linearly $L^{2}$-stable under a CFL-1 condition. In order to prevent the spurious oscillations that appear close to discontinuities two shock-capturing techniques have been considered: a flux-limiter technique (FL-CAT methods) and WENO reconstruction for the first time derivative (WENO-CAT methods). We follow [14] in the second approach.

These new methods have been compared in a number of test cases with WENO-RK methods (Finite Differences WENO reconstructions in space, TVD-RK in time) and with the WENO-LW methods introduced in [14] (Finite Differences WENO reconstruction for the first time derivative, Approximate Taylor in time). The linear transport equation, Burgers equation, and the 1D compressible Euler system have been considered.For $CFL\leq 0.5$ all the numerical methods work correctly, and the results obtained with all the methods using WENO reconstructions are similar, while the FL-CAT method is more diffusive as expected. Nevertheless, CAT methods are more expensive in computational time and number of operations due to its local character (FL-CAT is less expensive than WENO-CAT as reconstructions are avoided). However, the extra computational cost of CAT methods is compensated by the fact that they still give good solutions with CFL values close to 1.

Future developments include:

Parallel implementation.

• 2.

Use of fast WENO reconstructions: see [1]

• 3.

Order adaptive CAT methods based on smooth indicators.

• 4.

Application to systems of balance laws.

• 5.

Extension to multidimensional problems.

7 Acknowledgements

“This project has received funding from the European Union’s Horizon 2020 research and innovation program, under the Marie Sklodowska-Curie grant agreement No 642768”.