Some preliminary results on a high order asymptotic preserving computationally explicit kinetic scheme

TL;DR

This paper presents a method to construct high-order, explicit kinetic schemes that are asymptotic preserving and stable for discontinuities, maintaining accuracy for smooth solutions with computational costs similar to explicit schemes.

Contribution

It introduces a novel high-order kinetic scheme on regular meshes that is asymptotic preserving, explicit, and includes a nonlinear stability method for discontinuities.

Findings

Scheme can be arbitrarily high order in space and time

Runs at least CFL one, with explicit computational costs

Stable for problems with discontinuities without losing accuracy

Abstract

In this short paper, we intend to describe one way to construct arbitrarily high order kinetic schemes on regular meshes. The method can be arbitrarily high order in space and time, run at least CFL one, is asymptotic preserving and computationally explicit, i.e., the computational costs are of the same order of a fully explicit scheme. We also introduce a non linear stability method that enables to simulate problems with discontinuities, and it does not kill the accuracy for smooth regular solutions.

Tables

Table 1. Table 1: List of Fourier symbols.

Operator	Symbol g
${\delta }_1$	$1-e^{-i\theta }$
${\delta }_2$	$\frac{1}{3}{e}^{i\theta }+\frac{1}{2}-e^{-i\theta }+\frac{1}{6}{e}^{-2i\theta }$
${\delta }_4^1$	$i\left(\frac{\sin \left(2\theta \right)}{6}+\frac{4}{3}\sin \theta \right)$
${\delta }_4^2$	$\frac{e^{i\theta }}{4}+\frac{5}{6}-\frac{3}{2}{e}^{-i\theta }+\frac{1}{2}{e}^{-2i\theta }-\frac{e^{-3i\theta }}{12}$

Table 2. Table 2: Stability conditions for the original scheme.

	First order	second order	Fourth order
${\delta }_1$	$✓$	$✓$	$✓$
${\delta }_2$	$✓$	$✓$	$✓$ if $\lambda \le 4.5$
${\delta }_4^1$	$✓$	$✓$ ($\left|G\right|=1$)	$✓$ ($\left|G\right|=1$)
${\delta }_4^2$	$✓$	$✓$	$✓$ if $\lambda \le \frac{9}{4}$
Analytical condition	$\lambda \Re \left(g\right)\ge 0$	$\lambda \Re \left(g\right)\ge 0$	$\lambda \Re \left(g-\frac{\lambda g^2}{6}\right)>0$

Table 3. Table 3: CFL number for stability of the DeC iterations. 0 0 means that the scheme is unconditionally unstable. If a real number x 𝑥 x is given, it means that the scheme is stable up to CFL x 𝑥 x , if ≥ x absent 𝑥 \geq x is written, this means that the scheme is stable for at least CFL x 𝑥 x (and slightly above indeed).

Scheme		$\#$ iterations
Order	$\delta$	1	2	3	4	5	6
2	${\delta }_1$	1	1	1	1	1	1
2	${\delta }_2$	0	$\ge 0.85$	$\ge 1.22$	$\ge 1.02$	$\ge 1.08$	$\ge 1.23$
2	${\delta }_4^1$	0	0	$\ge 1.45$	$\ge 1.45$	$\ge 0.002$	$\ge 0.01$
2	${\delta }_4^2$	0	$\ge 0.5$	$\ge 0.69$	$0.71$	$0.73$	$0.73$
3	${\delta }_1$	$6$	$\ge 1.5$	$\ge 1.87$	$\ge 2$	$\ge 2.23$	$\ge 2.48$
3	${\delta }_2$	0	0	$1$	$\ge 2.0447$	$\ge 2.17120$	$\ge 2.568$
3	${\delta }_4^1$	0	0	0	$\ge 1.6171$	$\ge 2.4727$	$\ge 2.9162$
3	${\delta }_4^2$	0	0	$\ge 0.1$	$\ge 1.3096$	$\ge 1.3955$	$\ge 1.8282$

Table 4. Table 4: Order of convergence for the convection problem and two wave model for order 1, 2 and 4, and 4th order with MOOD. The final time is T = 0.5 𝑇 0.5 T=0.5 . One can see that the two fourth order results are identical as expected.

First order
$h$	$L^1$	$r$	$L^2$	$r$	$L^{\infty }$	$r$
50	${\mathrm{2.75963992\hspace{0.33em}10}}^{-2}$	-	${\mathrm{3.89822088\hspace{0.33em}10}}^{-2}$	-	${\mathrm{2.49972343\hspace{0.33em}10}}^{-2}$	-
100	${\mathrm{1.31966826\hspace{0.33em}10}}^{-2}$	1.43	${\mathrm{1.86553914\hspace{0.33em}10}}^{-2}$	1.43	${\mathrm{1.18893785\hspace{0.33em}10}}^{-2}$	1.43
200	${\mathrm{6.94037229\hspace{0.33em}10}}^{-3}$	1.33	${\mathrm{9.81052034\hspace{0.33em}10}}^{-3}$	1.33	${\mathrm{6.26028096\hspace{0.33em}10}}^{-3}$	1.33
400	${\mathrm{3.47535103\hspace{0.33em}10}}^{-3}$	1.38	${\mathrm{4.91373939\hspace{0.33em}10}}^{-3}$	1.38	${\mathrm{3.13191721\hspace{0.33em}10}}^{-3}$	1.38
800	${\mathrm{1.73895701\hspace{0.33em}10}}^{-3}$	1.38	${\mathrm{2.45900149\hspace{0.33em}10}}^{-3}$	1.38	${\mathrm{1.56634545\hspace{0.33em}10}}^{-3}$	1.38
Second order
50	${\mathrm{4.83627617\hspace{0.33em}10}}^{-3}$	-	${\mathrm{6.70434069\hspace{0.33em}10}}^{-3}$	-	${\mathrm{4.40502120\hspace{0.33em}10}}^{-3}$	-
100	${\mathrm{1.21754361\hspace{0.33em}10}}^{-3}$	2.07	${\mathrm{1.70489808\hspace{0.33em}10}}^{-3}$	2.06	${\mathrm{1.10206485\hspace{0.33em}10}}^{-3}$	2.08
200	${\mathrm{3.05118738\hspace{0.33em}10}}^{-4}$	2.07	${\mathrm{4.29379230\hspace{0.33em}10}}^{-4}$	2.07	${\mathrm{2.75470491\hspace{0.33em}10}}^{-4}$	2.08
400	${\mathrm{7.60697367\hspace{0.33em}10}}^{-5}$	2.08	${\mathrm{1.07314205\hspace{0.33em}10}}^{-4}$	2.08	${\mathrm{6.85840860\hspace{0.33em}10}}^{-5}$	2.08
800	${\mathrm{1.89899602\hspace{0.33em}10}}^{-5}$	2.08	${\mathrm{2.68216700\hspace{0.33em}10}}^{-5}$	2.08	${\mathrm{1.71091069\hspace{0.33em}10}}^{-5}$	2.08
Fourth order
50	${\mathrm{0.201424\hspace{0.33em}10}}^{-4}$	-	${\mathrm{0.278979\hspace{0.33em}10}}^{-4}$	-	${\mathrm{0.183213\hspace{0.33em}10}}^{-4}$	-
100	${\mathrm{0.122376\hspace{0.33em}10}}^{-5}$	4.04	${\mathrm{0.171337\hspace{0.33em}10}}^{-5}$	4.02	${\mathrm{0.110818\hspace{0.33em}10}}^{-5}$	4.04
200	${\mathrm{0.758547\hspace{0.33em}10}}^{-7}$	4.01	${\mathrm{0.106742\hspace{0.33em}10}}^{-6}$	4.00	${\mathrm{0.684850\hspace{0.33em}10}}^{-7}$	4.01
400	${\mathrm{0.472475\hspace{0.33em}10}}^{-8}$	4.00	${\mathrm{0.666515\hspace{0.33em}10}}^{-8}$	4.00	${\mathrm{0.425979\hspace{0.33em}10}}^{-8}$	4.00
800	${\mathrm{0.294831\hspace{0.33em}10}}^{-9}$	4.00	${\mathrm{0.416433\hspace{0.33em}10}}^{-9}$	4.00	${\mathrm{0.265631\hspace{0.33em}10}}^{-9}$	4.00
Fourth order+MOOD
50	${\mathrm{0.201424\hspace{0.33em}10}}^{-4}$	-	${\mathrm{0.278979\hspace{0.33em}10}}^{-4}$	-	${\mathrm{0.183213\hspace{0.33em}10}}^{-4}$	-
100	${\mathrm{0.122376\hspace{0.33em}10}}^{-5}$	4.04	${\mathrm{0.171337\hspace{0.33em}10}}^{-5}$	4.02	${\mathrm{0.110818\hspace{0.33em}10}}^{-5}$	4.04
200	${\mathrm{0.758547\hspace{0.33em}10}}^{-7}$	4.01	${\mathrm{0.106742\hspace{0.33em}10}}^{-6}$	4.00	${\mathrm{0.684850\hspace{0.33em}10}}^{-7}$	4.01
400	${\mathrm{0.472475\hspace{0.33em}10}}^{-8}$	4.00	${\mathrm{0.666515\hspace{0.33em}10}}^{-8}$	4.00	${\mathrm{0.425979\hspace{0.33em}10}}^{-8}$	4.00
800	${\mathrm{0.294831\hspace{0.33em}10}}^{-9}$	4.00	${\mathrm{0.416433\hspace{0.33em}10}}^{-9}$	4.00	${\mathrm{0.265631\hspace{0.33em}10}}^{-9}$	4.00

Table 5. Table 5: Order of convergence for the convection problem and two wave model for order 1, 2 and 4 with MOOD. The final time is T = 10 𝑇 10 T=10 . The fourth order results with and without stabilisation are identical as expected.

First order
$h$	$L^1$	$r$	$L^2$	$r$	$L^{\infty }$	$r$
50	$0.374521524$	-	$0.529551387$	-	$0.337824076$	-
100	$0.222015068$	1.22	$0.313846916$	1.22	$0.200003594$	1.21
200	$0.121430904$	1.23	$0.171709701$	1.30	$0.109344706$	1.3
400	${\mathrm{6.35740533\hspace{0.33em}10}}^{-2}$	1.34	${\mathrm{8.99051651\hspace{0.33em}10}}^{-2}$	1.34	${\mathrm{5.72395548\hspace{0.33em}10}}^{-2}$	1.34
800	${\mathrm{3.25359367\hspace{0.33em}10}}^{-2}$	1.36	${\mathrm{4.60121371\hspace{0.33em}10}}^{-2}$	1.36	${\mathrm{2.92929020\hspace{0.33em}10}}^{-2}$	1.36
Second order
50	${\mathrm{9.76886451\hspace{0.33em}10}}^{-2}$	-	$0.135240585$	-	${\mathrm{8.88576061\hspace{0.33em}10}}^{-2}$	-
100	${\mathrm{2.43498404\hspace{0.33em}10}}^{-2}$	2.08	${\mathrm{3.40924263\hspace{0.33em}10}}^{-2}$	2.07	${\mathrm{2.20487341\hspace{0.33em}10}}^{-2}$	2.09
200	${\mathrm{6.06694631\hspace{0.33em}10}}^{-3}$	2.08	${\mathrm{8.53730459\hspace{0.33em}10}}^{-3}$	2.08	${\mathrm{5.47759095\hspace{0.33em}10}}^{-3}$	2.08
400	${\mathrm{1.51354610\hspace{0.33em}10}}^{-3}$	2.08	${\mathrm{2.13514664\hspace{0.33em}10}}^{-3}$	2.08	${\mathrm{1.36459176\hspace{0.33em}10}}^{-3}$	2.08
800	${\mathrm{3.77953198\hspace{0.33em}10}}^{-4}$	2.08	${\mathrm{5.33837010\hspace{0.33em}10}}^{-4}$	2.08	${\mathrm{3.40518804\hspace{0.33em}10}}^{-4}$	2.08
Fourth order
50	${\mathrm{0.399627\hspace{0.33em}10}}^{-3}$	-	${\mathrm{0.554329\hspace{0.33em}10}}^{-3}$	-	${\mathrm{0.363964\hspace{0.33em}10}}^{-3}$	-
100	${\mathrm{0.244527\hspace{0.33em}10}}^{-4}$	4.03	${\mathrm{0.342394\hspace{0.33em}10}}^{-4}$	4.01	${\mathrm{0.221427\hspace{0.33em}10}}^{-4}$	4.03
200	${\mathrm{0.151613\hspace{0.33em}10}}^{-5}$	4.01	${\mathrm{0.213344\hspace{0.33em}10}}^{-5}$	4.00	${\mathrm{0.136893\hspace{0.33em}10}}^{-5}$	4.01
400	${\mathrm{0.944521\hspace{0.33em}10}}^{-7}$	4.00	${\mathrm{0.133241\hspace{0.33em}10}}^{-6}$	4.00	${\mathrm{0.851587\hspace{0.33em}10}}^{-7}$	4.00
800	${\mathrm{0.589189\hspace{0.33em}10}}^{-9}$	4.00	${\mathrm{0.832214\hspace{0.33em}10}}^{-8}$	4.00	${\mathrm{0.530836\hspace{0.33em}10}}^{-8}$	4.00
Fourth order with Mood
50	${\mathrm{0.399627\hspace{0.33em}10}}^{-3}$	-	${\mathrm{0.554329\hspace{0.33em}10}}^{-3}$	-	${\mathrm{0.363964\hspace{0.33em}10}}^{-3}$	-
100	${\mathrm{0.244527\hspace{0.33em}10}}^{-4}$	4.03	${\mathrm{0.342394\hspace{0.33em}10}}^{-4}$	4.01	${\mathrm{0.221427\hspace{0.33em}10}}^{-4}$	4.03
200	${\mathrm{0.151613\hspace{0.33em}10}}^{-5}$	4.01	${\mathrm{0.213344\hspace{0.33em}10}}^{-5}$	4.00	${\mathrm{0.136893\hspace{0.33em}10}}^{-5}$	4.01
400	${\mathrm{0.944521\hspace{0.33em}10}}^{-7}$	4.00	${\mathrm{0.133241\hspace{0.33em}10}}^{-6}$	4.00	${\mathrm{0.851587\hspace{0.33em}10}}^{-7}$	4.00
800	${\mathrm{0.589189\hspace{0.33em}10}}^{-9}$	4.00	${\mathrm{0.832214\hspace{0.33em}10}}^{-8}$	4.00	${\mathrm{0.530836\hspace{0.33em}10}}^{-8}$	4.00

Table 6. Table 6: Convection problem: error for the second order scheme with different ε 𝜀 \varepsilon

$\epsilon$	0		${10}^{-6}$		${10}^{-4}$		${10}^{-3}$		${10}^{-2}$
$\log \Delta x$	$\log L^2$	slope	$\log L^2$	slope	$\log L^2$	slope	$\log L^2$	slope	$\log L^2$	slope
-2.995	-1.332	-	-1.332	-	-1.332	-	-1.400	-	-2.049	-
-3.688	-2.655	1.908	-2.655	1.908	-2.655	1.907	-2.710	1.889	-3.314	1.825
-4.382	-4.030	1.984	-4.031	1.984	-4.031	1.984	-4.063	1.951	-4.647	1.922
-5.075	-5.415	1.997	-5.415	1.997	-5.415	1.996	-5.410	1.943	-5.999	1.951
-5.768	-6.801	1.999	-6.801	1.999	-6.800	1.998	-6.749	1.930	-7.365	1.971

Table 7. Table 7: Convection problem: error for the fourth order scheme with different ε 𝜀 \varepsilon

$\epsilon$	0		${10}^{-6}$		${10}^{-4}$		${10}^{-3}$		${10}^{-2}$
$\log \Delta x$	$\log L^2$	slope	$\log L^2$	slope	$\log L^2$	slope	$\log L^2$	slope	$\log L^2$	slope
-2.995	-3.615	-	-3.615	-	-3.701	-	-3.701	-	-4.445	-
-3.688	-6.385	3.995	-6.385	3.996	-6.475	3.996	-6.475	4.001	-7.212	3.992
-4.382	-9.158	4.000	-9.159	4.001	-9.253	4.000	-9.253	4.007	-9.991	4.008
-5.075	-11.93	4.000	-11.93	3.999	-12.03	4.000	-12.03	4.007	-12.76	4.004
-5.768	-14.70	4.000	-14.70	4.000	-14.80	4.000	-14.81	4.006	-15.53	4.001

Table 8. Table 8: CFL number for stability of the DeC iterations given by Gauss-Seidel and Jacobi methods. 0 0 means that the scheme is unconditionally unstable. x 𝑥 x means that the scheme is stable up to CFL x 𝑥 x , ≥ x absent 𝑥 \geq x means that the scheme is stable for at least CFL x (and slightly above).

Iterations	1	2	3	4	5
Symbol	Gauss Seidel
$g_1$	1.5	$1.276906714$	$1.167201858$	$1.197067146$,	$1.152628955$
$g_2$	0	$\ge 1.65$	$\ge 1.47$	$\ge 1.435$	$\ge 1.55$
$g_4^1$222always 1 for $x=\pi$	0	$\ge 0.926$	$\ge 1.775$	0	0
$g_4^2$	0	$\ge 0.917$	$0.8754013933$	$\ge 0.89$	$\ge 0.86$
Symbol	Jacobi
$g_1$	$1$	$1$	$1.256372663$	$1.392646782$	$1.774161172$
$g_2$	$0$	$\ge 0.87$	$\ge 1.625$	$\ge 1.744$	$\ge 2.06$
$g_4^1$ 333always 1 for $x=\pi$	$0$	$0$	$\ge 1.25$	$\ge 2.06$	$\ge 2.52$
$g_4^2$	$0$	$0$	$\ge 0.905$	$\ge 1.044$	$\ge 1.321$

Equations

$$\frac{\partial \mathbf{u}}{\partial t}+\frac{\partial \mathbf{f}(\mathbf{u})}{\partial x}=0$$

$$\mathbf{u}(x,0)={\mathbf{u}}_0(x),$$

$$\frac{\partial \mathbf{F}}{\partial t}+\Lambda \frac{\partial \mathbf{F}}{\partial x}=\frac{\mathbb{M}(\mathbb{P}\mathbf{F})-\mathbf{F}}{\epsilon }$$

$$\begin{array}[]{l}\dfrac{\partial u}{\partial t}+\dfrac{\partial v}{\partial x}=0\\ \dfrac{\partial v}{\partial t}+a^{2}\dfrac{\partial u}{\partial x}=\dfrac{f(u)-v}{\varepsilon}$$

\end{array}$$

$$\begin{array}{rl}\frac{\partial f_1}{\partial t}+a\frac{\partial f_1}{\partial x}& =\frac{{\mathbb{M}}_1-f_1}{\epsilon },\\ \frac{\partial f_2}{\partial t}-a\frac{\partial f_2}{\partial x}& =\frac{{\mathbb{M}}_2-f_2}{\epsilon },\end{array}$$

$$\mathbf{F}=\left(\begin{array}{c}{f}_1\\ f_2\end{array}\right),\Lambda =\left(\begin{array}{cc}a& 0\\ 0& -a\end{array}\right),\mathbb{P}\mathbf{F}=f_1+f_2\text{ and }\mathbb{M}=\left(\begin{array}{c}{\mathbb{M}}_1\\ {\mathbb{M}}_2\end{array}\right)$$

$${\mathbb{M}}_1+{\mathbb{M}}_2=f_1+f_2=u,a({\mathbb{M}}_1-{\mathbb{M}}_2)=f(u),$$

$${\mathbb{M}}_1(f,a)=\frac{1}{2}\left(f_1+f_2+\frac{f(u)}{a}\right),{\mathbb{M}}_2(f,a)=\frac{1}{2}\left(f_1+f_2-\frac{f(u)}{a}\right).$$

$$\frac{\partial \mathbf{F}}{\partial t}+\frac{1}{\Delta x}\Lambda \delta \mathbf{F}=\frac{\mathbb{M}(\mathbb{P}\mathbf{F})-\mathbf{F}}{\epsilon }+O(\Delta x^q).$$

$$\frac{dU}{dt}+\mathcal{F}(U)=\frac{\mathcal{G}(U)}{\epsilon }$$

$$\begin{array}[]{c|c}c&A\\ \hline\cr 0&b^{T}\end{array}\quad\text{and }\quad\begin{array}[]{c|c}\tilde{c}&\tilde{A}\\ \hline\cr 0&\tilde{b}^{T}\end{array}$$

$$\begin{array}{rl}{U}_0& =U^n\\ ⋮& \\ U_k& =U_0+\Delta t\sum _{j=1}^{k-1}{a}_{kj}\mathcal{F}(U_j)+\frac{\Delta t}{\epsilon }\sum _{j=1}^s{\tilde{a}}_{kj}\mathcal{G}(U_j)\\ ⋮& \\ U^{n+1}& =U^n+\Delta t\sum _{j=1}^s{b}_j\mathcal{F}(U_j)+\frac{\Delta t}{\epsilon }\sum _{j=1}^s{\tilde{b}}_j\mathcal{G}(U_j)\end{array}$$

$${\int }_{t_n}^{t_n+c_i\Delta t}\phi (s)ds\approx \Delta t\sum _{j=0}^q{a}_{ij}\phi (t_n+c_j\Delta t).$$

$$a_{ij}={\int }_0^{c_i}{\ell }_j(s)ds,$$

$$\sum _{j=0}^q{a}_{ij}=c_i.$$

$$\mathbf{F}_{k}^{n,j}-\mathbf{F}_{k}^{n,0}+\dfrac{\Delta t}{\Delta x}\bigg{(}\sum_{l=0}^{q}a_{il}\Lambda\delta_{k}\mathbf{F}^{n,l}\bigg{)}-\mu\sum_{l=0}^{q}a_{il}\big{(}\mathbb{M}\mathbb{P}\mathbf{F}_{k}^{n,l}-\mathbf{F}_{k}^{n,l}\big{)}=0,\quad j=1,\ldots,q$$

$$\mathbf{\mathcal{F}}_{k}=\big{(}\mathbf{F}_{k}^{n,1},\ldots,\mathbf{F}_{k}^{n,q}\big{)}^{T},\quad\mathbf{\mathcal{F}}_{k}^{(0)}=\big{(}\mathbf{F}_{k}^{n,0},\ldots,\mathbf{F}_{k}^{n,0}\big{)}^{T}=\big{(}\mathbf{F}_{k}^{n},\ldots\mathbf{F}_{k}^{n}\big{)},$$

$$\begin{split}\mathbf{\mathcal{F}}_{k}\!\!-\mathbf{\mathcal{F}}_{k}^{(0)}+\dfrac{\Delta t}{\Delta x}\Lambda A\delta_{k}\mathbf{\mathcal{F}}\!&-\mu A\big{(}\mathbb{M}(\mathbb{P}\mathbf{\mathcal{F}}_{k})\!-\mathbf{\mathcal{F}}_{k}\big{)}\!+\!\dfrac{\Delta t}{\Delta x}\Lambda\mathbf{a}_{0}\otimes\delta_{k}\mathbf{F}^{n,0}\!\!\\ &-\mu\mathbf{a}_{0}\otimes(\mathbb{M}(\mathbb{P}\mathbf{F}_{k}^{n,0})\!-\mathbf{F}_{k}^{n,0})\!=0,\end{split}$$

$$\mathbb{M}(\mathbb{P}\mathbf{\mathcal{F}})=\big{(}\mathbb{M}(\mathbb{P}\mathbf{F}^{1}),\ldots,\mathbb{M}(\mathbb{P}\mathbf{F}^{q})\big{)}^{T}.$$

$$A=\left(\begin{array}{ccc}{a}_{11}& \dots & a_{1q}\\ ⋮& ⋮& ⋮\\ a_{q1}& \dots & a_{qq}\end{array}\right)$$

$${\mathbf{a}}_0=\left(\begin{array}{c}{a}_{0q}\\ ⋮\\ a_{01}\end{array}\right).$$

$$\mathbf{F}_{k}^{n,j}-\mathbf{F}_{k}^{n}+c_{j}\dfrac{\Delta t}{\Delta x}\Lambda\delta_{k}\mathbf{F}^{n,0}-\mu\sum_{l=0}^{q}a_{jl}\big{(}\mathbb{M}(\mathbb{P}\mathbf{F}_{k}^{n,l})-\mathbf{F}_{k}^{n,l}\big{)}=0,\quad j=1,\ldots,q.$$

$$\mathbf{\mathcal{F}}_{k}-\mathbf{\mathcal{F}}_{k}^{(0)}+\dfrac{\Delta t}{\Delta x}C\Lambda\delta_{k}\mathbf{\mathcal{F}}^{(0)}-\mu A\big{(}\mathbb{M}(\mathbb{P}\mathbf{\mathcal{F}}_{k})-\mathbf{\mathcal{F}}_{k}\big{)}-\mu\mathbf{a}_{0}\otimes(\mathbb{M}\mathbb{P}\mathbf{F}_{k}^{n,0}-\mathbf{F}_{k}^{n,0})=0$$

$$\big{[}\mathcal{L}^{1}(\mathbf{\mathcal{F}})\big{]}_{k}:=\mathbf{\mathcal{F}}_{k}-\mathbf{\mathcal{F}}_{k}^{(0)}+\dfrac{\Delta t}{\Delta x}C\Lambda\delta_{k}\mathbf{\mathcal{F}}^{(0)}-\mu A\big{(}\mathbb{M}(\mathbb{P}\mathbf{\mathcal{F}}_{k})-\mathbf{\mathcal{F}}_{k}\big{)}-\mu\mathbf{a}_{0}\otimes(\mathbb{M}(\mathbb{P}\mathbf{F}^{n,0}_{k})-\mathbf{F}^{n,0}_{k}),$$

$$\begin{split}\big{[}\mathcal{L}^{2}(\mathbf{\mathcal{F}})\big{]}_{k}:=&\mathbf{\mathcal{F}}_{k}-\mathbf{\mathcal{F}}_{k}^{(0)}+\dfrac{\Delta t}{\Delta x}\Lambda A\delta_{k}\mathbf{\mathcal{F}}-\mu A\big{(}\mathbb{M}(\mathbb{P}\mathbf{\mathcal{F}}_{k})-\mathbf{\mathcal{F}}_{k}\big{)}+\dfrac{\Delta t}{\Delta x}\Lambda\mathbf{a}_{0}\otimes\delta_{k}\mathbf{F}^{n,0}\\ &-\mu\mathbf{a}_{0}\otimes(\mathbb{M}(\mathbb{P}\mathbf{F}^{n,0}_{k})-\mathbf{F}^{n,0}_{k}).\end{split}$$

$${\delta }_k\mathcal{F}={\hat{\mathcal{F}}}_{k+1/2}-{\hat{\mathcal{F}}}_{k-1/2}$$

$$\mathcal{L}^{2}(\mathbf{\mathcal{F}})_{k}-\mathcal{L}^{1}(\mathbf{\mathcal{F}})_{k}=\dfrac{\Delta t}{\Delta x}\Lambda A\big{(}\delta_{k}\mathbf{\mathcal{F}}-\delta_{k}\mathbf{\mathcal{F}}^{(0)}\big{)},$$

$${\mathcal{L}}^1({\mathcal{F}}^{(p+1)})={\mathcal{L}}^1({\mathcal{F}}^{(p)})-{\mathcal{L}}^2({\mathcal{F}}^{(p)}),$$

$$\Vert {\mathcal{F}}^{(q)}-{\mathcal{F}}^{\star }\Vert \le C\Delta t^q,$$

$$\mathbb{P}{\mathcal{F}}_k=\mathbb{P}{\mathcal{G}}_k+\mathbb{P}{\mathcal{F}}_k^{(0)}-\frac{\Delta t}{\Delta x}\mathbb{P}C\Lambda {\delta }_k{\mathcal{F}}^{(0)}=\mathbb{P}{\mathcal{G}}_k+\mathcal{K},$$

$$\mathcal{K}=\mathbb{P}{\mathcal{F}}_k^{(0)}-\frac{\Delta t}{\Delta x}\mathbb{P}C\Lambda {\delta }_k{\mathcal{F}}^{(0)}.$$

Full text

{AMS}

65M12; 65L04; 65M60

1 Introduction

Let us specify first the context. We are given the PDE

[TABLE]

with $\mathbf{u}\in\mathbb{R}^{p}$ and ${\mathbf{f}}:\mathbb{R}^{p}\rightarrow\mathbb{R}^{p}$. a Lipschitz continuous flux It is known, at least since the work of Jin [1] and then Natalini [2] and co-workers, that this system can formally be seen as the limit for $\varepsilon\rightarrow 0$ of a relaxation system:

[TABLE]

with $\mathbf{F}\in\mathbb{R}^{k\times p}$, $\mathbb{M}$ is a Maxwellian and $\mathbb{P}$ is a linear operator such that $\mathbb{P}\mathbb{M}(\mathbb{P}\mathbf{F})=\mathbb{P}\mathbf{F}$. The constant matrix $\Lambda$ and the flux ${\mathbf{f}}$ are linked by $\mathbb{P}\Lambda\mathbb{M}(\mathbb{P}\mathbf{F})={\mathbf{f}}(\mathbb{P}\mathbf{F})$. The simplest example, due to Jin and Xin [1], is

[TABLE]

that can be rewritten in the form (1.2) with:

[TABLE]

i.e., where

[TABLE]

where the Maxwellian is defined from the relations

[TABLE]

i.e.

[TABLE]

We know that $a$ must be larger than the max of $|f^{\prime}(u)|$ because of the Whitham sub-characteristic condition, obtained via a formal Chapman Enskog expansion. Another argument is, as shown by [3], that under this condition the two Maxwellian $\mathbb{M}_{1}$ and $\mathbb{M}_{2}$ satisfy a monotonicity condition, i.e. the BGK model becomes compatible with entropy inequalities.

The questions we address in this paper are the following: given a system (1.1) and a regular grid of spatial step $\Delta x>0$, can we construct a computationally explicit scheme that solves (1.2) with uniform accuracy of order $r>0$ for all $\varepsilon>0$ and with a CFL condition, based on the matrix $\Lambda$, that is larger than $1$111 Initially, the first author was motivated by understanding in a better way the LBM method, even though the answer is not about the LBM method at all. The only remaining property between what we look for and the LBM method is the CFL condition.. The answer is yes, and this paper proposes a simple construction in one dimension. With computationally explicit we mean that the solution of a certain scheme does not require any nonlinear solver, nor the inversion of a mass matrix.

High order accurate methods for kinetic problems à la Shi-Jin has received a lot of attention in the recent years. For a long time the state of the art was that of second second order in time and space finite volume with TVD like stabilisation, see e.g [4]. For higher than second order, one may mention [5] where a splitting approach is adopted with a regular CFL stability condition for the overall finite volume scheme, [6] where relaxed upwind schemes are proposed running up to CFL $=1$ and up to third order in time/space, again in a finite volume context. In [7], a WENO approach is proposed. In [8] a discontinuous Galerkin approximation of the system (1.2) is developed (with a temporal scheme allowing very large CFL number). In the kinetic literature where the fluid system is represented with the BGK approximation, so that dense and less dense flows can be simulated, there has also been a large effort towards high order schemes with asymptotic preserving properties. One may mention [9] for hyperbolic systems with diffusion, [10] where a high order conservative semi-Lagrangian technique is developed.

We want to go beyond that, with very simple and cheap numerical schemes that are potentially arbitrary high order and run at CFL $=1$, with an accuracy that is independent of the relaxation parameter $\varepsilon$. The format of the paper is as follows. We first introduce the general method which amounts to describe the discretisation of $\Lambda\dfrac{\partial\mathbf{F}}{\partial x}$ and a time discretisation. We take into account the source term. The scheme resulting from this discretisation is fully implicit. The next step is to show that, thanks to the operator $\mathbb{P}$, and using a particular time discretisation, we can make it computationally explicit, and high order accurate, independently of the parameter $\varepsilon$. Several choices of $\Lambda$ and Maxwellians $\mathbb{M}$ are described. We also address the question of the non linear stabilisation of the method when discontinuities appear. Several numerical examples, covering scalar and system cases, are then proposed to show the relevance of the method. The accuracy is checked for the scalar case.

2 General discretisation principle

Starting from (1.2), the idea is to discretise first in space $\Lambda\dfrac{\partial\mathbf{F}}{\partial x}$. This introduces an error which we assume to be $O(\Delta x^{q})$,

[TABLE]

The second step is to discretise in time, so that we expect that the resulting scheme will be of order $p$ in space and time, at least for moderate values of $\varepsilon$. The problem is then two-fold: (i) how to define the discretisation operator $\delta$ for which a minimum requirement is the semi discrete linear stability when there is no source term, (ii) how to discretise in time so that the accuracy is uniform in time and $\varepsilon$. We first discuss the issue of time discretisation, then space discretisation.

2.1 Time discretisation

One may use IMEX Runge-Kutta schemes, and more precisely SSP IMEX Runge-Kutta schemes, to have a better control of the stability properties of the method. Rewriting (1.2a) as the sum of a non stiff term and a stiff one

[TABLE]

an IMEX method is defined by two Butcher’s tableaux

[TABLE]

where the first one is for non stiff part, while the second one is for the stiff part:

[TABLE]

with various compatibility conditions so that a given order is reached, see [11, Chapter IV]. Anticipating a bit, if there exists a linear operator $\mathbb{P}$ such that $\mathbb{P}\mathcal{G}=0$ as here, we see that, applying $\mathbb{P}$ to (2.6), a necessary condition is that the explicit RK scheme defined by the explicit part is itself SSP. Since we want to have a running CFL number of at least one, this needs that the SSP RK scheme must have a CFL number of at least $1+\epsilon$, $\epsilon>0$. To our knowledge there are some explicit SSP RK schemes satisfying this condition, inter alia [12], but they are not generalizable to arbitrarily high order of accuracy, and no IMEX versions are available.

For this reason, we use an IMEX deferred correction (DeC) method. It is a general way of building arbitrarily high order Runge Kutta schemes. It also allows more freedom in the spatial discretization, for instance, allowing the use of lumped mass matrix [13]. Its implicit and IMEX versions allow to use a combination of more traditional low order IMEX schemes and arbitrarily high order implicit RK schemes, obtaining arbitrarily high order IMEX schemes. We leave the study of SSP version of these schemes for future research. The final IMEX DeC scheme we obtain is computationally explicit and it is also matrix-free.

2.1.1 Deferred Correction

The DeC is an iterative procedure that was proposed and developed in its explicit version in [14] and in an implicit version in [15]. It was applied to hyperbolic PDE, for instance, in [13], with a new formalism which makes the proof of its properties more straightforward. An IMEX version of this algorithm applied to hyperbolic PDE is available in [16], and the algorithm we discuss in the following is a modification of this one. With the notation of [13], the DeC uses two operators: one high order accurate $\mathcal{L}^{2}$, which defines a fully implicit method, and a low order easy to solve $\mathcal{L}^{1}$ operator. The process allows to approximate with arbitrary accuracy the solution of the high order operator $\mathcal{L}^{2}$, with the simplicity of the operator $\mathcal{L}^{1}$. We start with the description of the high order operator $\mathcal{L}^{2}$.

Let us consider $q+1$ points in $[0,1]$, $c_{0}=0<c_{1}\ldots<c_{i}<\ldots<c_{q}=1$ and the quadrature formula

[TABLE]

More precisely, if $\{\ell_{j}\}$ are the Lagrange polynomials associated to the partition $\{c_{j}\}_{j=0}^{q}$, if we take

[TABLE]

the quadrature formula is of order $q+1$. We will always require that the quadrature formula are consistent, i.e.

[TABLE]

Considering $x_{k}$, a grid point, and setting $\mathbf{F}_{k}^{n,j}\approx\mathbf{F}(x_{k},t_{n}+c_{j}\Delta t)$ and $\mathbf{F}_{k}^{n,0}=\mathbf{F}_{k}^{n}$, an approximation of (1.2) is:

[TABLE]

where $\mu=\tfrac{\Delta t}{\varepsilon}$ and $\tfrac{\delta\mathbf{F}}{\Delta x}$ is a consistent approximation of $\dfrac{\partial\mathbf{F}}{\partial x}$. We will set $\mathbf{F}_{k}^{n+1}=\mathbf{F}_{k}^{n,q}$. The relations (2.8) can be rewritten in matrix form, setting

[TABLE]

and neglecting the index of the timestep $n$, as

[TABLE]

where by abuse of language we have written

[TABLE]

The matrix $A$ is

[TABLE]

and we have

[TABLE]

As a result, (2.9) is implicit, and in general non linear, because of the Maxwellian. In order to simplify the resolution, we consider a simpler scheme, where the source term discretisation remains the same and the forward Euler method is used on each sub-time step:

[TABLE]

We rewrite this as:

[TABLE]

where C=\text{ diag }\big{(}c_{1},\ldots,c_{q}\big{)} and $\mathbf{\mathcal{F}}^{(0)}=(\mathbf{F}^{n,0},\ldots,\mathbf{F}^{n,0})^{T}$.

This leads to the introduction of two operators $\mathcal{L}^{1}$ and $\mathcal{L}^{2}$ acting on $\mathbf{\mathcal{F}}=(\ldots,\mathbf{\mathcal{F}}_{k},\mathbf{\mathcal{F}}_{k+1},\ldots)$ and defined as:

[TABLE]

and

[TABLE]

So that (2.11) is $\mathcal{L}^{1}(\mathbf{\mathcal{F}}^{n,j})_{k}=0$ while (2.9) is $\mathcal{L}^{2}(\mathbf{\mathcal{F}}^{n,j})_{k}=0$. In order to have more structure, we will require that $\delta_{k}\mathbf{\mathcal{F}}$ has the following difference form:

[TABLE]

where $\widehat{\mathbf{\mathcal{F}}}_{k+1/2}$ depends on $P$ arguments, is consistent with $\mathbf{\mathcal{F}}$ and uniformly Lipschitz continuous with respect to its arguments. Examples will be given in section 2.2.

Thanks to (2.7), we see that

[TABLE]

the important fact is that $\varepsilon$ plays no role here.

We will solve the problem (2.9) with the following defect correction (DeC) method:

• •

Set, for any $k$, $\mathbf{\mathcal{F}}_{k}^{(0)}=(\mathbf{F}_{k}^{n},\ldots,\mathbf{F}_{k}^{n})^{T}$,

• •

Solve for $p=0,\ldots,M-1$ the problem

[TABLE]

• •

Set $\mathbf{\mathcal{F}}^{n+1}=\mathbf{\mathcal{F}}^{(M)}.$

We remark that the operator $\mathcal{L}^{2}$ will never be solved directly as it will be applied to the previously computed iteration $\mathbf{\mathcal{F}}^{(p)}$. The DeC procedure will converge to the solution of the $\mathcal{L}^{2}(\mathbf{\mathcal{F}}^{*})=0$ problem by solving iteratively (2.15). We show that if the problem (2.9) has a unique solution $\mathbf{\mathcal{F}}^{\star}$ and taking $M=q$, we have a formal error of $\Delta t^{q}$, i.e., for a norm to be defined,

[TABLE]

so that the formal accuracy is the same as solving exactly (2.9). Before doing that, we have first to explain how we solve for $\mathcal{L}^{1}$ and, hence, (2.15), then we show the error estimate (and define the proper norm).

2.1.2 Solution of $\mathcal{L}^{1}(\mathbf{\mathcal{F}})=\mathbf{\mathcal{G}}$ and (2.15).

Let us first start with $\mathcal{L}^{1}(\mathbf{\mathcal{F}})=\mathbf{\mathcal{G}}$ for any $\mathcal{G}\in$. Applying $\mathbb{P}$ to this equation, we get, for any $k\in\mathbb{Z}$,

[TABLE]

with

[TABLE]

The found equation is explicit for $\mathbb{P}\mathbf{\mathcal{F}}$ and we can in practice compute this term and use it to obtain the solution of the whole operator. Substituting $\mathbb{P}\mathbf{\mathcal{F}}$ into the Maxwellian, we obtain

[TABLE]

where all the unknown terms $\mathbf{\mathcal{F}}$ depends only linearly on some coefficients that we can collect on the left hand side

[TABLE]

Now, if $\text{Id}_{q\times q}+\mu A$ is invertible, we can compute only once and store its inverse to obtain an easy solution of the whole problem, i.e.,

[TABLE]

In this way, we are able to find a solution of the system $\mathcal{L}^{1}(\mathbf{\mathcal{F}})=\mathcal{G}$ in a computationally explicit way: the source term is split into the linearly implicit part and the Maxwellian evaluated in the previously computed $\mathbb{P}\mathbf{\mathcal{F}}$. The remaining terms are the explicit right hand side $\mathcal{G}$, the explicit convection term $\delta_{k}\mathbf{\mathcal{F}}^{(0)}$ and the explicit part of the high order time integration of the source. There have been many other works using similar techniques in order to explicitly solve the implicit discretization of the kinetic system (1.1b), inter alia [17, 4, 18, 16]. This DeC approach combining the two operators, has the advantage of being arbitrarily high order without building complicated structures. Indeed, we can use this computationally explicit solution to solve (2.15).

Setting

[TABLE]

we can apply (2.16) directly. However, since the source term discretisation is the same, we have simplifications. Indeed, after little algebra, (2.15) is rewritten as:

[TABLE]

and we can apply the same technique. We first apply $\mathbb{P}$,

[TABLE]

so we know explicitly $\mathbb{P}\mathbf{\mathcal{F}}^{(p+1)}$, and then we can solve

[TABLE]

and then

[TABLE]

Again, we see that the method is computationally explicit.

Next, we show the error estimate, and then we comment more on (2.19), in particular when $\varepsilon\rightarrow 0$.

2.1.3 Error estimate

If $\varphi:\mathbb{R}\rightarrow\mathbb{R}^{p}$ is $C^{1}(\mathbb{R})^{p}$ and has a compact support, we can consider the discrete version of its $L^{2}$ and $H^{1}$ norms:

[TABLE]

where $D_{i}\varphi=\tfrac{\varphi_{i+1}-\varphi_{i}}{\Delta x}.$

We will establish error estimates that are valid in a given but arbitrary compact $I=[a,b]$ with discrete equivalent of $L^{2}_{loc}$ and $H^{-1}_{loc}$ estimates:

[TABLE]

and we note that for $\varphi\in C^{1}_{0}([a,b])^{p}$, we have a Poincaré like inequality

[TABLE]

We first show that

Lemma 2.1**.**

If $\widehat{\mathbf{\mathcal{F}}}_{k+1/2}=\sum\limits_{l=-p}^{q}\alpha_{l}\mathbf{\mathcal{F}}_{k+l}$ and letting $C=\max\limits_{-p\leq l\leq q}|\alpha_{l}|\times\max\limits_{i}|\lambda_{i}|$, we have

[TABLE]

Proof 2.2**.**

We have, from (2.14) and since $\delta_{k}\mathbf{\mathcal{F}}=\widehat{\mathbf{\mathcal{F}}}_{k+1/2}-\widehat{\mathbf{\mathcal{F}}}_{k-1/2}$, we have, using that $\varphi$ has a compact support,

[TABLE]

We remark that the norm of the coefficients $A$ is smaller or equal to 1.

We also have the following lemma on $\mathcal{L}_{1}$:

Lemma 2.3**.**

We assume that the Maxwellian is Lipschitz continuous and that there exists $C,C^{\prime}>0$ such that for all $\varepsilon>0$,

[TABLE]

Let us consider $\mathbf{\mathcal{F}},\mathbf{\mathcal{F}}^{\prime}$ and $\mathbf{\mathcal{G}},\mathbf{\mathcal{G}}^{\prime}$ such that

[TABLE]

Then, there exists $\alpha>0$, independent of $\mathbf{\mathcal{F}},\mathbf{\mathcal{F}}^{\prime}$, $\varepsilon$ and $I$ such that

[TABLE]

and

[TABLE]

Proof 2.4**.**

We have the explicit solution $\mathbf{\mathcal{F}}$ and $\mathbf{\mathcal{F}}^{\prime}$ from (2.16), and we see that

[TABLE]

so that if $\|~{}.~{}\|$ is any of the two norms, we have

[TABLE]

The constant $\alpha$ depends on $C$, $C^{\prime}$, $\Lambda$ and the Lipschitz constant of the Maxwellian. Remember also that in (2.16), $\mathcal{K}$ depends only on $\mathbf{\mathcal{F}}_{0}$, $\frac{\Delta t}{\Delta x}$ and $\Lambda$. It is independent of $\varepsilon$.

Then, wrapping all together, we have the following proposition:

Proposition 2.5**.**

Under the assumptions of lemma 2.1 and 2.3, if $\mathbf{\mathcal{F}}^{\star}$ is the unique solution of $\mathcal{L}^{2}(\mathbf{\mathcal{F}})=0$, there exists $\theta$ independent of $\varepsilon$ such that we have, for all $p\in\mathbb{N}$

[TABLE]

Proof 2.6**.**

*We first have, since $\mathcal{L}^{2}(\mathbf{\mathcal{F}}^{\star})=0$ *

[TABLE]

so that combining the inequalities of lemmas 2.1, 2.3 and the Poincaré like inequality, we obtain the result.

Remark 2.7** (Comments about inequality (2.21)).**

This result shows that after $p+1$ iteration, the error is $O(\Delta x^{p+1})=O(\Delta t^{p+1})$ if a CFL-like condition is available. Of course it is better that $\theta\Delta x<1$ for the inequality to be effective, so we may experience a reduction of the CFL number. This reduction needs to be studied case by case, however this also show that the overall cost of the method is of the order of an explicit one. This is why we name this computationally explicit.

2.1.4 Asymptotic preservation

We can show that the presented method is asymptotic preserving (AP), starting from the Chapman–Enskog expansion of the model (1.2a). Let us define $\mathbf{u}^{\varepsilon}:=\mathbb{P}\mathbf{F}$, we obtain that

[TABLE]

Proposition 2.8**.**

The discretisation given by (2.19) is consistent with the limit model (2.22) up to an $\mathcal{O}(\varepsilon)$.

Proof 2.9**.**

Now, using first (2.19) and then (2.18), defining $\mathbf{u}^{(p),\varepsilon}=\mathbb{P}\mathbf{\mathcal{F}}^{(p)}$ and recalling that $\mu=\frac{\Delta t}{\varepsilon}$, $\mathbb{M}(\mathbb{P}\mathbf{u})=\mathbf{u}$, $\mathbb{P}\Lambda\mathbb{M}(\mathbf{u})={\mathbf{f}}(\mathbf{u})$, by induction on the subtimesteps $p$

[TABLE]

we can extend the formal expansion also in the discrete case, i.e.,

[TABLE]

The final result is a discretisation in space and time of the asymptotic model given by (2.22), if the the spatial discretisation is consistent with the space derivative.

Remark 2.10**.**

One can proceed further and prove that, both in the discrete and the continuous case, the next term of the Chapman Enskog expansion is a diffusive term under Whitham’s subcharacteristic condition of $\Lambda^{2}-\partial_{u}{\mathbf{f}}(u)$ being positive definite. We can also prove that the discretisation is consistent also with that term up to an $\mathcal{O}(\varepsilon^{2})+\mathcal{O}(\Delta t^{2})$ if the spatial discretisation is at least consistent. We refer to [16] for the details of such computations for the sake of brevity.

2.1.5 Examples of $\mathcal{L}^{2}$ time discretisation.

Here, we will consider second and fourth order approximation in time in the $\mathcal{L}^{2}$ operator, namely the Crank-Nicholson method and the fourth order one that uses the points $t_{n}$, $t_{n}+\tfrac{\Delta t}{2}$ and $t_{n+1}$. They are described by their matrices $A$,

• •

Second order

[TABLE]

with $\mathbf{F}^{n,0}=\mathbf{F}^{n}$ and $\mathbf{F}^{n,1}\approx\mathbf{F}(t^{n+1})$. Writing the $\mathcal{L}^{2}$ discretisation of the time derivative applied to $\mathbf{F}$, we would have

[TABLE]

i.e.,

[TABLE]

This is Crank-Nicholson. We see that

[TABLE]

are uniformly bounded in $[0,2]$ for any $(\Delta t,\varepsilon)$.

• •

The fourth order scheme is obtained by

[TABLE]

where $\mathbf{F}^{n,0}=\mathbf{F}^{n}$, $\mathbf{F}^{n,1}\approx\mathbf{F}(t_{n}+\frac{\Delta t}{2})$ and $\mathbf{F}^{n,2}\approx\mathbf{F}(t^{n+1}).$ We see that

[TABLE]

so the matrix is invertible. It is also easy to see that the matrices

[TABLE]

are uniformly bounded.

In fact, in this case, the operator $\mathcal{L}^{2}$ corresponds to the scheme Lobatto III, which is 4th order accurate [11]. For that reason, we will use this temporal scheme in conjunction with a fourth order spatial approximation.

2.2 Space discretisation: Definition of the $\delta$ operator

The only question left is how to define a stable scheme. As we have seen in section 2.1.4, the scheme is asymptotic preserving. Under Whitham’s subcharacteristic conditions the relaxation terms introduces diffusion which further stabilize the scheme. Hence, we focus only on the stability of the convection scheme, which will also guarantee the stability of the full scheme. The answer for the fully non linear convection problem is out of reach, at least for this paper, so we will rely on a classical linear stability analysis. The stability of the convection schemes splits into two sub-questions: is the convection scheme defined by $L^{2}=0$ conditionally or unconditionally stable, and then, is the convection scheme defined by the DeC iteration (2.16) stable, and under which conditions. In the next section, we will provide 3 examples with increasing accuracy, and sketch a general method.

The matrix $\Lambda$ is diagonal. In [19], the author considers the transport equation

[TABLE]

and shows that if $a<0$ and

[TABLE]

then the order is at most $2\min(r+1,s)$ and in addition the only stable methods are those defined for $r=s$ or $s=r+1$ or $s=r+2$. If $a>0$, we set

[TABLE]

while in that case $r=s$ or $r=s+1$ or $r=s+2$. We will only consider these approximations. Following [19], we have

[TABLE]

and

[TABLE]

Remark 2.11** (Conservation).**

We note that we can always write

[TABLE]

with

[TABLE]

Proof 2.12**.**

Assuming that $\hat{\mathbf{f}}_{i+1/2}=\sum\limits_{j=-r}^{s-1}\beta_{j}u_{i+j}$ for any $i$, we write

[TABLE]

so that $\beta_{j}=-\sum\limits_{l\geq j+1}\alpha_{l}$, using that $\sum\limits_{l=-s}^{r}\alpha_{l}=0$.

This means that the approximations (2.16) and (2.19), in the limit $\varepsilon\rightarrow 0$ is always conservative since $\Lambda$ is diagonal, and thanks to (2.13).

We list some possible choices for $\delta$:

• •

First order approximation: this is the upwind scheme. If $a>0$, we take $\delta_{1}u_{j}=u_{j}-u_{j-1}$, while if $a<0$, $\delta_{1}u_{j}=u_{j+1}-u_{j}$. If $a=0$, of course $\delta_{1}u_{j}=0$. The flux is

[TABLE]

• •

Second order: for $a<0$,

[TABLE]

so that

[TABLE]

with

[TABLE]

This corresponds to the $[r,r+2]$ approximation with $r=-1$. In term of flux, we have (for $a<0$):

[TABLE]

For $a>0$, we have

[TABLE]

so all in all

[TABLE]

• •

Fourth order: if $r=s=2$, and for any $a$

[TABLE]

hence

[TABLE]

and if $r=1$, $s=3$ and $a<0$,

[TABLE]

In term of flux, we have:

• –

for $\delta_{4}^{1}$,

[TABLE]

• –

for $\delta_{4}^{2}$,

[TABLE]

3 Stability analysis

We study the stability of the discretisation of the homogeneous problem. Since $\Lambda$ is diagonal, it is enough to look at the scalar conservation problem. We first look at the implicit method defined by $\mathcal{L}^{2}=0$, and then at the DeC iteration that is constructed on top of it. This is done by Fourier analysis, we can assume that $a>0$ and the Fourier symbol of $\delta$ is $g$. The table 1 displays the symbols of the operators.

The next step is to evaluate the amplification factors of the method, first without DeC iteration, then with DeC iteration.

3.1 First order in time

For a first order scheme the $\mathcal{L}^{2}$ operator can be written as an implicit Euler method, though being computationally explicit, while the DeC iteration, which consists of one step, resemble the explicit Euler method with CFL constrained $0\leq\lambda\leq 1$, where $\lambda=a\Delta t/\Delta x$. For the $\mathcal{L}^{2}=0$ operator, by Fourier transform, we have $\hat{u}^{n+1}-\hat{u}^{n}+\lambda g\hat{u}^{n+1}=0$, so that the amplification factor is $G=\dfrac{1}{1+\lambda g}$ which is of modulus $\leq 1$ if

[TABLE]

If $\lambda\rightarrow 0^{+}$, we see that $\Re(g)\geq 0$ is a necessary condition, while if $\lambda\rightarrow 0^{-}$, $\Re(g)\leq 0$. In all cases, $\lambda\Re(g)\geq 0$ is a necessary condition. Writing $g=a+ib$, and assuming that $\lambda\neq 0$, we see that this condition reads:

[TABLE]

We also see that

[TABLE]

so that $\lambda\Re(g)\geq 0$ is a necessary and sufficient condition for stability. The Table 2 provides the stability condition for the first, second and fourth order schemes. For the rest of the discussion we consider $a=1$ and in case it is different, one has to rescale $\lambda$ accordingly, as classically done for CFL conditions.

3.2 Second order in time

In that case the $\mathcal{L}^{2}=0$ scheme reads:

[TABLE]

for which the amplification factor is simply

[TABLE]

We have $|G|\leq 1$ if and only if

[TABLE]

Again, the Table 2 provides the stability condition for the first, second and fourth order schemes.

The DeC iteration is

[TABLE]

so that

[TABLE]

and we see that

[TABLE]

3.3 Fourth order in time

Here, the $\mathcal{L}^{2}=0$ scheme reads:

[TABLE]

so that the Fourier transform gives

[TABLE]

with

[TABLE]

and we have to look at $\max\{|G_{1}|,|G_{2}|\}\leq 1$ for the calculation of $\hat{u}^{n+1/2}$ and $\hat{u}^{n+1}$ to be stable. We have

[TABLE]

Then with obvious notations, the DeC iteration is

[TABLE]

The Fourier analysis gives:

[TABLE]

The amplification vector, $G_{p}$ after the $p$-th iteration is defined by

[TABLE]

We note that, setting $\mathbf{\theta}=\begin{pmatrix}\theta_{1}^{1}&\theta_{2}^{1}\\ \theta_{1}^{2}&\theta_{2}^{2}\end{pmatrix}$,

[TABLE]

and $\rho(\mathbf{\theta})=\dfrac{1}{2\sqrt{3}}$. So using the spectra decomposition of $\theta$ which has two complex and distinct eigenvalues, we have that

[TABLE]

We get finally $\mu_{1}=0.4115783562$, $\mu_{2}=0.1188124373$, $\mu_{3}=0.03429819635$, $\mu_{4}=0.009901036444$, hence the convergence is very quick.

3.4 Summary of the stability analysis

Combining these expressions with the actual form of the Fourier symbol of $\delta$, we get the results of table 2.

Now, we turn our attention on the DeC iteration. For the second order in time approximation, we first have

[TABLE]

so we get

[TABLE]

hence if $|G|\leq 1$, a sufficient condition is that

[TABLE]

For the fourth order scheme, we have similarly

[TABLE]

but it is more complicated to get an analytical condition. So we rely on Maple.

The stability conditions are summarised in Table 3.

4 Wave model

We have to specify the diagonal matrix $\Lambda$ and the Maxwellians $\mathbb{M}$. We will use two kinds of wave models:

• •

A two waves model. In that case,

[TABLE]

with $a\geq\max_{i}\rho({\mathbf{f}}^{\prime}(u_{i}))$. Setting $\mathbf{u}^{\varepsilon}=\mathbb{P}\mathbf{F}$, we have ${\mathbf{f}}(\mathbf{u}^{\varepsilon})=\mathbb{P}\Lambda\mathbf{F}$ and we know explicitly $\mathbb{M}=(\mathbb{M}_{1},\mathbb{M}_{2})$:

[TABLE]

• •

A three waves model, where

[TABLE]

In that case, the Maxwellian is $\mathbb{M}=(\mathbb{M}_{1},\mathbb{M}_{2},\mathbb{M}_{3})$ and we have

[TABLE]

so we need to specify $\mathbb{M}_{2}$.

For the scalar problems, we will use the two wave model that reveals itself sufficient. For the fluid problems, we will show that the two wave models is not perfect, and hence the three wave model needs to be considered.

In the case of 3 waves, let us specify $\mathbb{M}_{2}$. Following [3], we know that the sub-characteristic condition is equivalent to the monotonicity of the Maxwellians: they need to be differentiable and have only positive eigenvalues. In [3, 2], it is proposed to use

[TABLE]

where ${\mathbf{f}}(\mathbf{u}^{\varepsilon})={\mathbf{f}}_{+}(\mathbf{u}^{\varepsilon})+{\mathbf{f}}_{-}(\mathbf{u}^{\varepsilon})$, ${\mathbf{f}}_{\pm}$ are differentiable, $\nabla_{\mathbf{u}}{\mathbf{f}}_{+}(\mathbf{u})$ has only positive eigenvalues, while $\nabla_{\mathbf{u}}{\mathbf{f}}_{-}(\mathbf{u}^{\varepsilon})$ has only negative eigenvalues. A possible choice, inspired by the Enquist-Osher-Solomon flux, is

[TABLE]

but the integral (or the path integral for system) must be evaluated. In the case of the Euler equations, we give a second one that does not necessitate to evaluate an integral. In the case of the Euler equations,

[TABLE]

we propose to use a Maxwellian that relies on the van Leer flux splitting [20]. It is purely algebraic and defined by:

if $M=\dfrac{u}{c}\leq-1$, with $c^{2}=\gamma\frac{p}{\rho}$, then ${\mathbf{f}}_{-}(\mathbf{u}^{\varepsilon})={\mathbf{f}}(\mathbf{u}^{\varepsilon})$, ${\mathbf{f}}_{+}(\mathbf{u}^{\varepsilon})=0$, 2. 2.

If $M\geq 1$, then ${\mathbf{f}}_{+}(\mathbf{u}^{\varepsilon})={\mathbf{f}}(\mathbf{u}^{\varepsilon})$, ${\mathbf{f}}_{-}(\mathbf{u}^{\varepsilon})=0$, 3. 3.

if $-1\leq M\leq 1$, then

[TABLE]

and ${\mathbf{f}}_{+}={\mathbf{f}}-{\mathbf{f}}_{-}$

The eigenvalues of ${\mathbf{f}}_{\pm}$ are bounded by

[TABLE]

Note that $\tfrac{\gamma+3}{2\gamma+|M|(3-\gamma)}\leq\tfrac{\gamma+3}{2\gamma}$ for $|M|\leq 1$. For $\gamma=1.4$, $\tfrac{\gamma+3}{2\gamma}=\tfrac{11}{7}\approx 1.57$.

5 Non linear stabilisation

If the solution is expected to be non smooth, then one can expect the occurrence of spurious oscillations. Sometimes, oscillations are acceptable, provided they do not lead to the crash of the simulation. In order to get rid of them, or to control them, we have adopted the MOOD technique initially designed in [21] with some improvements described in [22]. We have adapted it our way in order to get results that are formally of order $p+1$ in space and time, here $p=1,2,3$.

MOOD is an *a posteriori *corrector of high order numerical methods [21, 22]. MOOD requires a sequence of schemes ordered from the most accurate/less stable one to the low order/more reliable one. It also requires a series of criteria that the solution should fulfill, e.g. physical admissibility, discrete minimum principle or numerical errors. After having performed a step of the most accurate scheme, it checks the criteria on each cell/degree of freedom, and detects the areas where the criteria are not met. There, we switch to the next scheme in a cascade style, which is supposed to be more stable and reliable. We proceed iteratively until either the criteria are met or the most reliable/less accurate parachute scheme is used. The parachute scheme should analytically guarantee all the criteria.

In the following we describe how the criteria must be verified on the described spatial discretisation, while the list of the schemes that we use consists always of 2 schemes (the considered one and the upwind discretisation $\delta_{1}$ as parachute scheme) and we specify directly in the numerical simulations which criteria will be considered, as they are problem dependent.

We proceed as follows: at the time step $t_{n}$, we have the values $(\mathbf{F}_{k}^{n})_{k}$. For now on, we drop the superscript $n$, since there is no ambiguity. In the DeC iteration (2.8), with the spatial scheme defined by $\delta_{p}$, writes (with the convention that $\mathbf{F}^{(l),0}=\mathbf{F}^{0}$ for $l=0,\ldots,q-1$)

[TABLE]

from which we get

[TABLE]

The increment $\delta_{k}\mathbf{F}^{l}$ is the difference of two terms, and we write

[TABLE]

with

[TABLE]

One equivalent way to rephrase the conservation is

[TABLE]

and the right hand side of this relation is independent of the order $p$. It is equivalent because we see that

[TABLE]

Using this we rewrite (5.30) as:

[TABLE]

In practice, we compute for each interval $[k,k+1]$

[TABLE]

and then apply (5.32a).

In the simplified version of the MOOD algorithm we use, we consider only two spatial approximations, namely the first order one defined by $\delta_{1}$, and the high order one define by $\delta_{p}$, $p=2$ or $3$ in this paper. The idea is to use as often as possible the highest order scheme, and to use the low order one to correct potential problems. Knowing the $\{\mathbf{F}_{k}^{(p),j}\}_{k}$, we first compute for each interval $[k,k+1]$ the quantities defined by (5.33) with the high order residuals. Then we test the values of the results using a set of criteria, applied on $(\widetilde{\mathbb{P}\mathbf{F}_{k}^{(p+1),j}})_{k+1/2}$ and $(\widetilde{\mathbb{P}\mathbf{F}_{k+1}^{(p+1),j}})_{k+1/2}$. This set of criteria is explained in the next paragraph. If both $(\widetilde{\mathbb{P}\mathbf{F}_{k}^{(p+1),j}})_{k+1/2}$ and $(\widetilde{\mathbb{P}\mathbf{F}_{k+1}^{(p+1),j}})_{k+1/2}$ pass the tests, this element is declared sane, else un-sane. This enable to identify a set $\mathcal{I}$ of un-sane elements $[k,k+1]$ where the criteria are not met, and we store the residual $\{\Phi_{k}^{[k,k+1],(p),j},\Phi_{k+1}^{[k,k+1],(p),j}\}$ for the sane elements. We then repeat the procedure for the un-sane elements with the lowest order scheme. At the end of the procedure, we have evaluated residuals, that we still denote by $\{\Phi_{k}^{[k,k+1],(p),j},\Phi_{k+1}^{[k,k+1],(p),j}\}$, even though they are potentially evaluated by different schemes. We then compute $\mathbb{P}\mathbf{F}_{k}^{(p+1),j}$ by (5.30). There is no problem of conservation since (5.31) holds true.

Now we describe the criteria we apply to $(\widetilde{\mathbb{P}\mathbf{F}_{k}^{(p+1),j}})_{k+1/2}$ and $(\widetilde{\mathbb{P}\mathbf{F}_{k+1}^{(p+1),j}})_{k+1/2}$, following the ideas of [21, 22] with some small adaptation to the context. When specific tests are done on a variable, we denote this variable by $\xi$. For a scalar problem, $\xi$ is simply the conserved variable. In the case of the Euler equations, we test this on some primitive variables: the density and the energy, and for some severe problems, the velocity. We can add as many criteria as needed.

We first check if $(\widetilde{\mathbb{P}\mathbf{F}_{k}^{(p+1),j}})_{k+1/2}$ and $(\widetilde{\mathbb{P}\mathbf{F}_{k+1}^{(p+1),j}})_{k+1/2}$ lies in the invariance domain if relevant: in the case of the Euler equation, we check if the density and the internal energy are both positive. If not, we set the criteria to .FALSE. on this element. In that case we jump to the next element, else we look for the next criterion. 2. 2.

We check if the solution is not locally constant. Taking $\nu=\Delta x^{3}$ and $\mathcal{S}$ the stencil defined by the operator $\delta$, we check if

[TABLE]

If this is true, the criteria is kept to .TRUE., else it is set to .FALSE. and we jump to the next element 3. 3.

We check if a new extrema is created or not, by comparing with the solution at the previous time step, in a neighbourhood extended to the right and the left by one cell: we are running at CFL 1.

1. (a)

We first test if $\xi_{k}^{n+1},\xi_{k+1}^{n+1}\in[\min\limits_{l\in\mathcal{S}}\xi_{k+l}+\epsilon,\max\limits_{l\in\mathcal{S}}\xi_{k+l}-\epsilon]\cap[\min\limits_{l\in\mathcal{S}}\xi_{k+1+l}+\epsilon,\max\limits_{l\in\mathcal{S}}\xi_{k+1+l}-\epsilon]$. If this is true, we jump to the next element, 2. (b)

else, denoting by $P_{j}$ the Lagrange interpolation polynomial that interpolates $\{\xi_{j+l}\}_{l\in\mathcal{S}}$

• •

we compute $\xi^{\prime}=P_{k}^{\prime}(x_{k})$, $\xi^{\prime}_{L}=P^{\prime}_{k}(x_{k}-\frac{\Delta x}{2})$, \xi^{\prime,k-1/2}_{min/max}=\min/\max\big{(}P^{\prime}_{k}(x_{k}-\frac{\Delta x}{2}),P^{\prime}_{k-1}(x_{k}-\frac{\Delta x}{2})\big{)} then

• –

If $\xi^{\prime}_{L}<\xi^{\prime}$, $\alpha_{L}=\min(1,\dfrac{\xi^{\prime,k-1/2}_{max}-\xi^{\prime}}{\xi^{\prime}_{L}-\xi^{\prime}})$

• –

if $\xi^{\prime}_{L}=\xi^{\prime}$, $\alpha_{L}=1$

• –

if $\xi^{\prime}_{L}<\xi^{\prime}$, $\alpha_{L}=\min(1,\dfrac{\xi^{\prime,k-1/2}_{min}-\xi^{\prime}}{\xi^{\prime}_{L}-\xi^{\prime}})$

• •

We compute $\xi^{\prime}=P_{k}^{\prime}(x_{k})$, $\xi^{\prime}_{R}=P^{\prime}_{k}(x_{k}+\frac{\Delta x}{2})$, \xi^{\prime,k+1/2}_{min/max}=\min/\max\big{(}P_{k+1}^{\prime}(x_{k}+\frac{\Delta x}{2}),P^{\prime}_{k}(x_{k}+\frac{\Delta x}{2})\big{)} then

• –

If $\xi^{\prime}_{R}<\xi^{\prime}$, $\alpha_{R}=\min(1,\dfrac{\xi^{\prime,k+1/2}_{max}-\xi^{\prime}}{\xi^{\prime}_{R}-\xi^{\prime}})$

• –

if $\xi^{\prime}_{R}=\xi^{\prime}$, $\alpha_{R}=1$

• –

if $\xi^{\prime}_{R}<\xi^{\prime}$, $\alpha_{R}=\min(1,\dfrac{\xi^{\prime,k+1/2}_{min}-\xi^{\prime}}{\xi^{\prime}_{R}-\xi^{\prime}})$

• •

We set $\alpha=\min(\alpha_{L},\alpha_{R})$

• •

If $\alpha=1$, then we have a true extrema, keep the criteria to .TRUE. and jump to the next element. Else, we set the criteria to .FALSE. and jump to the next element.

The idea behind the step 3 is described in [22] and is also related to [23]: we try to check if the gradient of the interpolation $\xi$ lies in the interval \left[\min\big{(}\xi^{\prime,k-1/2}_{min},\xi^{\prime,k+1/2}_{min}\big{)},\max\big{(}\xi^{\prime,k-1/2}_{max},\xi^{\prime,k+1/2}_{max}\big{)}\right].

Remark 5.1** (Stability).**

The von Neumann stability study of Section 3 does not hold directly on the MOOD algorithm but it is clear that we are dealing with a combination of the stabilities of the schemes used in the MOOD cascade. Hence, if we choose CFL conditions that guarantee the stability of both the high order schemes and of the parachute scheme (upwind CFL=1 in our case), then we know that the global MOOD scheme will be von Neumann stable.

6 Numerical examples

6.1 Scalar problems

The first problem is the transport equation with periodic boundary conditions

[TABLE]

where the initial condition is

[TABLE]

A two waves model is used with $a=1.01$ (so a little larger that the actual maximum speed. We always proceed as such for scalar and system cases. We make a convergence test for short and long final times, namely $T=0.5$ and $T=10$. The CFL number, with respect to the wave model maximum speed, is always set to $1$. In both cases, we see that the expected order of accuracy is obtained. Note that for the second and fourth order schemes, the non linear stabilisation does not detect any troubled point, giving exactly the same error as in the non stabilised case. We show the results for the fourth order schemes. All the calculations are done with the two waves model. Note that the first order scheme, with the initial condition given by the Maxwellian, is nothing more than the Lax-Friedrichs scheme, for second order in time approximation. For fourth order in time, since the equilibrium relaxation is more complex, we get a different scheme. Note that the non linear stabilisation procedure of section 5 does not flag any cell.

The convergence tables are obtained against the solution of the asymptotic model, proving moreover that the model is asymptotic preserving as expected.

The figure 1 shows some results for the Burgers equation

[TABLE]

with the initial condition (6.34). This generates an unsteady shock wave, so a priori more challenging than a steady one. The non linear stabilisation performs correctly.

The last scalar example is the Buckley-Leverett equation

[TABLE]

again with the same initial condition (6.34). The flux is non convex, so the problem is a bit more challenging.

We have run the simulations with 100 spatial points, until time $T=1$, with the 2 waves model we have considered. The first order (O1), second order (O2), fourth order (O4), second order with non linear stabilisation (O2M), and fourth order with non linear stabilisation (O4M) are displayed in figure 2, together with a reference solution computed with 1000 points and the first order scheme: remember that this corresponds to the Lax Friedrichs scheme, and it satisfies all entropy inequalities. This guaranties that the scheme converges. The non linearly stabilized solution have a correct behavior.

6.2 Uniform order of accuracy with respect to $\varepsilon$

In order to test the convergence of the scheme also in non asymptotic regimes, we use again the scalar advection equation

[TABLE]

with initial condition

[TABLE]

The CFL number is set to $1$, the time and space order are set to 2 and 4, the final time is $T=1$, the boundary conditions are periodic. The values of $\varepsilon$ are $\{0,10^{-6},10^{-5},10^{-4},10^{-3},10^{-2}\}$. We evaluate the order of convergence with the following standard procedure: if $u_{\Delta x}$, $u_{\Delta x/2}$ are the numerical solution evaluated for consecutive meshes, the order $\alpha$ is, for the norm $\|~{}.~{}\|$,

[TABLE]

This test does not have obvious result, as order reduction phenomena are common for IMEX schemes when the space discretization and the relaxation variable are of the same order. Nevertheless, we see on tables 6 and 7 that the convergence order does not depend on $\varepsilon$. It is optimal. Also qualitatively, we see that for small enough values of $\varepsilon$, the solution obtained for $\varepsilon=0$ is almost indistinguishable from $0<\varepsilon\ll 1$.

6.3 Euler equations

In this section we test our scheme on Euler equations (4.29). We set $\gamma=1.4$ and we run some standard cases: the Sod case and the Shu-Osher case.

6.3.1 Sod test case

The Sod problem consists of a Riemann problem defined by the following initial conditions:

[TABLE]

The final time is $T=0.16$. We have used the 3-waves model described above. The mesh resolution is of 100 elements, and the CFL is again 1 in all cases. From figure 3, we see that the results are of good quality, at least compared with more standard methods.

For the sake of completeness, we have made the same simulation with the two wave model in figure 4.

We see a stair case solution of the first order in space which is typical for the Lax-Friedrichs scheme. Comparing the solutions, the 2 waves model provides results of lower quality with respect to the 3 waves one. For that reason, we will not consider it anymore for the Euler equations.

6.3.2 Shu-Osher test case

The conditions of the Shu-Osher test are

[TABLE]

on the domain $[-5,5]$ and the final time of the problem is $T=1.8$. The reference solution is obtained with 10.000 points and 4th order limited scheme. It is difficult to see any modification in the solution if we use more grid points, this is why we consider this solution as the reference solution. We display only the solution with the MOOD stabilisation technique, however, we have tried two different strategies. The figures labeled as OXMood, where X=2 or 4, use the full strategy of section 5. The physical variables are the density and the pressure, nothing is tested on the velocity. In the figures labeled as OXMoodNaN, with X=2 or 4, we only check if the solution lies in the invariance domain, i.e., density and pressure stay positive and we do not encounter NaN values. In figure 5 we plot the results for 800 points. In figure 6 we compare results for 200, $400$ and $800$ mesh points only for the MoodNaN strategy.

From figure 6, we see that with $800$ points, there is hardly no difference with between the O4MoodNaN solution and the reference one.

7 Conclusion

In this paper, simplifying a method described in [16], we show how to construct a class of kinetic numerical methods that can run at least at CFL 1. They can handle in a simple manner hyperbolic problems, and in particular compressible fluid mechanics one. These methods are always locally conservative and thus can handle correctly discontinuities. We have described a rather simple stabilisation mechanism which can be further improved or changed: it is not really the core of the proposed method. Our methodology can be arbitrarily high order and can use CFL number larger or equal to unity on regular Cartesian meshes. Extension to the multidimensional case will be the topic of future works. In particular, our implementation of these methods indicates that they can be potentially very fast. The parallelisation should also be straightforward. This, however, has to be confirmed in several spatial dimensions.

Acknowledgment

R.A. would like to thanks Prof. Li-Shi Luo (Old Dominion, USA) and Professor Sagaut (Aix-Marseilles University, France) for their encouragement. Finally, he would like to thank Professors Lallemand and D’Humières for the discussion we had during his PhD studies. This paper shows that this has never been forgotten. I take this opportunity to also thanks Prof. Chi-Wang Shu for very helpfull discussions about time stepping. D.T. has been partially funded by ITN ModCompShock project funded by the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No 642768.

Appendix A Another iteration technique for the fourth order in time case

The iteration (2.15) solved by (2.19) reads (when we have no source term or after the application of the projector operator):

[TABLE]

and after the application of the Fourier transform, we have

[TABLE]

so that the amplification factor satisfies

[TABLE]

This can be seen as the Jacobi iteration for solving the system

[TABLE]

By analogy, we can define a Gauss-Seidel iteration by

[TABLE]

whose Fourier transform is

[TABLE]

and hence

[TABLE]

In both cases, $G^{(0)}=e$.

We can study the stability of the Gauss–Seidel iteration, and we recall the results of Jacobi’s for comparison. Denoting by $g_{1}$ (resp. $g_{2}$, $g_{4,1}$, $g_{4,2}$) the Fourier symbol of the operators $a\delta_{1}$ (resp. $a\delta_{2}$, $a\delta_{4}^{1}$, $a\delta_{4}^{2}$), we get the results of table 8.

Remark A.1** (A few remarks about table 8.).**

**

• •

For $g_{4}^{1}$, the amplification factor is always equal to $1$ when $x=\pi$, and strictly below $1$ under the condition stated above.

• •

For $g_{1}$ and 3 iterations, the CFL condition can be computed exactly. It is $\frac{1}{2}\sqrt[3]{4+\sqrt{17}}-\frac{1}{2}\,{\frac{1}{\sqrt[3]{4+\sqrt{17}}}}+\frac{1}{2}\approx 1.256372663$.

From this results, we see that there is no fundamental reason to prefer Gauss-Seidel iteration to the Jacobi one; the coding of the Gauss-Seidel is also slightly more involved. However, this conclusion holds true only for the schemes we have considered here, and might not be true for others.