An entropic fourier method for the Boltzmann equation

TL;DR

This paper introduces an entropic Fourier method for the Boltzmann equation that combines the advantages of discrete velocity and Fourier methods, ensuring positivity, entropy decay, and computational efficiency.

Contribution

It develops a novel entropic Fourier approach that preserves positivity and the H-theorem while enabling fast Fourier transform algorithms for the Boltzmann collision operator.

Findings

Preserves positivity of solutions

Satisfies a discrete H-theorem

Achieves second-order convergence rate

Abstract

We propose an entropic Fourier method for the numerical discretization of the Boltzmann collision operator. The method, which is obtained by modifying a Fourier Galerkin method to match the form of the discrete velocity method, can be viewed both as a discrete velocity method and as a Fourier method. As a discrete velocity method, it preserves the positivity of the solution and satisfies a discrete version of the H-theorem. As a Fourier method, it allows one to readily apply the FFT-based fast algorithms. A second-order convergence rate is validated by numerical experiments

Tables

Table 1. Table 1: The ℓ 1 subscript ℓ 1 \ell_{1} , ℓ 2 subscript ℓ 2 \ell_{2} and ℓ ∞ subscript ℓ \ell_{\infty} errors and convergence rates for the BKW solution at time t = 0.01 𝑡 0.01 t=0.01 with R = 6 𝑅 6 R=6 .

$N$	${\ell }_1$ error	rate	${\ell }_2$ error	rate	${\ell }_{\infty }$ error	rate
16	$4.68\times {10}^{-3}$		$3.23\times {10}^{-3}$		$3.12\times {10}^{-3}$
32	$1.72\times {10}^{-3}$	1.44	$1.36\times {10}^{-3}$	1.25	$1.40\times {10}^{-3}$	1.15
64	$5.54\times {10}^{-4}$	1.64	$4.56\times {10}^{-4}$	1.58	$5.57\times {10}^{-4}$	1.34
128	$1.55\times {10}^{-4}$	1.84	$1.29\times {10}^{-4}$	1.82	$1.73\times {10}^{-4}$	1.68
256	$4.05\times {10}^{-5}$	1.93	$3.42\times {10}^{-5}$	1.92	$4.73\times {10}^{-5}$	1.87
512	$1.03\times {10}^{-5}$	1.97	$8.76\times {10}^{-6}$	1.96	$1.22\times {10}^{-5}$	1.94

Table 2. Table 2: The ℓ 1 subscript ℓ 1 \ell_{1} error of EFM with fast algorithm in [ 19 ] for multiple choices of N 𝑁 N and M 𝑀 M .

$N$	$M=2$	$M=3$	$M=32$
16	$4.6852\times {10}^{-3}$	$4.6826\times {10}^{-3}$	$4.6830\times {10}^{-3}$
32	$1.7241\times {10}^{-3}$	$1.7244\times {10}^{-3}$	$1.7245\times {10}^{-3}$
64	$5.5368\times {10}^{-4}$	$5.5388\times {10}^{-4}$	$5.5394\times {10}^{-4}$
128	$1.5485\times {10}^{-4}$	$1.5488\times {10}^{-4}$	$1.5489\times {10}^{-4}$
256	$4.0513\times {10}^{-5}$	$4.0516\times {10}^{-5}$	$4.0517\times {10}^{-5}$

Table 3. Table 3: The ℓ 1 subscript ℓ 1 \ell_{1} , ℓ 2 subscript ℓ 2 \ell_{2} and ℓ ∞ subscript ℓ \ell_{\infty} errors and convergence rates for the BKW solution at time t = 0.01 𝑡 0.01 t=0.01 with R = 6 𝑅 6 R=6 .

$N$	${\ell }_1$ error	rate	${\ell }_2$ error	rate	${\ell }_{\infty }$ error	rate
16	$4.08\times {10}^{-3}$		$3.08\times {10}^{-3}$		$3.56\times {10}^{-3}$
32	$1.42\times {10}^{-3}$	1.52	$1.12\times {10}^{-3}$	1.47	$1.26\times {10}^{-3}$	1.50
64	$4.07\times {10}^{-4}$	1.80	$3.29\times {10}^{-4}$	1.76	$3.72\times {10}^{-4}$	1.76
128	$1.08\times {10}^{-4}$	1.91	$8.85\times {10}^{-5}$	1.90	$1.00\times {10}^{-4}$	1.89

Equations

$$\frac{\partial f}{\partial t}+v\cdot {\nabla }_{x}f=\mathcal{Q}[f,f],f(0,x,v)=f^0(x,v)$$

$$\mathcal{Q}[f,f](v)={\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{S}}^{d-1}}\mathcal{B}(v-v_{\ast },\omega )\left[f(v^{\prime })f(v_{\ast }^{\prime })-f(v)f(v_{\ast })\right]\mathrm{d}{v}_{\ast }\mathrm{d}\omega$$

$$v^{\prime }=\frac{v+v_{\ast }}{2}+\frac{|v-v_{\ast }|}{2}\omega ,v_{\ast }^{\prime }=\frac{v+v_{\ast }}{2}-\frac{|v-v_{\ast }|}{2}\omega$$

$$\mathcal{B}(v-v_{\ast },\omega )=b(|v-v_{\ast }|,\cos \theta ),\cos \theta =|(v-v_{\ast })\cdot \omega |/|v-v_{\ast }|.$$

$${\int }_{{\mathbb{R}}^d}\psi (v)\mathcal{Q}[f,f](v)\mathrm{d}v=$$

$$\frac{1}{4}{\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{S}}^{d-1}}\left(\psi (v)+\psi (v_{\ast })-\psi (v^{\prime })-\psi (v_{\ast }^{\prime })\right)\mathcal{B}[f(v^{\prime })f(v_{\ast }^{\prime })-f(v)f(v_{\ast })]\mathrm{d}{v}_{\ast }\mathrm{d}v\mathrm{d}\omega .$$

$${\int }_{{\mathbb{R}}^d}\mathcal{Q}[f,f](v)\mathrm{d}v=0,{\int }_{{\mathbb{R}}^d}\mathcal{Q}[f,f](v)v\mathrm{d}v=0,{\int }_{{\mathbb{R}}^d}\mathcal{Q}[f,f](v)|v{|}^2\mathrm{d}v=0,$$

$${\int }_{{\mathbb{R}}^d}\mathcal{Q}[f,f](v)\ln (f(v))\mathrm{d}v=$$

$$\frac{1}{4}{\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{S}}^{d-1}}\ln \left(\frac{f(v)f(v_{\ast })}{f(v^{\prime })f(v_{\ast }^{\prime })}\right)\mathcal{B}[f(v^{\prime })f(v_{\ast }^{\prime })-f(v)f(v_{\ast })]\mathrm{d}{v}_{\ast }\mathrm{d}v\mathrm{d}\omega \le 0$$

$$\mathcal{Q}[f,f](v)={\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{R}}^d}\tilde{\mathcal{B}}(y,z)\delta (y\cdot z)[f(v+y)f(v+z)-f(v)f(v+y+z)]\mathrm{d}y\mathrm{d}z.$$

$$\tilde{\mathcal{B}}(y,z)=2^{d-1}\mathcal{B}\left(y+z,\frac{y-z}{|y-z|}\right)|y+z{|}^{2-d}.$$

$$\frac{\partial f}{\partial t}=\mathcal{Q}[f,f]$$

$$X=\{h\cdot k|k=(k_1,\dots ,k_d),-n\le k_1,\dots ,k_d\le n\}.$$

$$\frac{\mathrm{d}{F}_r(t)}{\mathrm{d}t}=Q_r(t):=\sum _{p,q,s\in X}{A}_{pq}^{rs}\left(F_p(t)F_q(t)-F_r(t)F_s(t)\right),r\in X.$$

$$A_{pq}^{rs}\ne 0\Rightarrow p+q=r+s\text{and}|p{|}^2+|q{|}^2=|r{|}^2+|s{|}^2$$

$$A_{pq}^{rs}=A_{qp}^{rs}=A_{rs}^{pq}.$$

$$\sum _{r\in X}{Q}_r{\psi }_r=\frac{1}{4}\sum _{p,q,r,s\in X}{A}_{pq}^{rs}\left({\psi }_r+{\psi }_s-{\psi }_p-{\psi }_q\right)\left(F_p{F}_q-F_r{F}_s\right).$$

$$\sum _{r\in X}{Q}_r=0,\sum _{r\in X}{Q}_{r}r=0,\sum _{r\in X}{Q}_r|r{|}^2=0.$$

$$\sum _{r\in X}{Q}_r\ln (F_r)=\frac{1}{4}\sum _{p,q,r,s\in X}{A}_{pq}^{rs}\ln \left(\frac{F_r{F}_s}{F_p{F}_q}\right)\left(F_p{F}_q-F_r{F}_s\right)\le 0$$

$${\mathcal{Q}}^R[f,f](v):={\int }_{B_R}{\int }_{B_R}\tilde{\mathcal{B}}(y,z)\delta (y\cdot z)[f(v+y)f(v+z)-f(v)f(v+y+z)]\mathrm{d}y\mathrm{d}z,$$

$$f(v)=\sum _{k\in {\mathbb{Z}}^d}{\hat{f}}_k{E}_k(v),{\hat{f}}_k=\frac{1}{(2T{)}^d}{\int }_{{\mathcal{D}}_T}f(v)E_{-k}(v)\mathrm{d}v,$$

$${\mathcal{Q}}^R[f,f](v)=\sum _{l,m\in {\mathbb{Z}}^d}\left(\hat{B}(l,m)-\hat{B}(m,m)\right){\hat{f}}_l{\hat{f}}_m{E}_{l+m}(v),v\in {\mathcal{D}}_T,$$

$$\hat{B}(l,m):={\int }_{B_R}{\int }_{B_R}\tilde{\mathcal{B}}(y,z)\delta (y\cdot z)E_l(y)E_m(z)\mathrm{d}y\mathrm{d}z,l,m\in {\mathbb{Z}}^d.$$

$$\hat{B}(l,m)=\hat{B}(m,l)=\hat{B}(l,-m).$$

$${\hat{\mathcal{Q}}}^R[f,f{]}_k=\sum _{l,m\in {\mathbb{Z}}^d}1(l+m-k)\left(\hat{B}(l,m)-\hat{B}(m,m)\right){\hat{f}}_l{\hat{f}}_m,k\in {\mathbb{Z}}^d.$$

$$K:=\{k|k=(k_1,\dots ,k_d),-n\le k_1,\dots ,k_d\le n\},$$

$${\mathbb{P}}_N=\mathrm{span}\left\{E_k(v)|k\in K\right\}\subset L_{\mathrm{per}}^2({\mathcal{D}}_T).$$

$${\hat{Q}}_k^{\mathsf{G}}:=\sum _{l,m\in K}1(l+m-k)\left(\hat{B}(l,m)-\hat{B}(m,m)\right){\hat{F}}_l{\hat{F}}_m,k\in K,$$

$$\{\begin{array}{rl}& \frac{\mathrm{d}{\hat{F}}_k}{\mathrm{d}t}={\hat{Q}}_k^{\mathsf{G}}=\sum _{l,m\in K}1(l+m-k)\left(\hat{B}(l,m)-\hat{B}(m,m)\right){\hat{F}}_l{\hat{F}}_m,\\ & {\hat{F}}_k(t=0)={\hat{F}}_k^0,\end{array}$$

$$\begin{array}{rl}\hat{B}(l,m)& ={\int }_{B_R}{\int }_{{\mathbb{S}}^{d-1}}\mathcal{B}(g,\omega )E_l\left(\frac{1}{2}(g+|g|\omega )\right)E_m\left(\frac{1}{2}(g-|g|\omega )\right)\mathrm{d}g\mathrm{d}\omega ,\end{array}$$

$$F_p=\sum _{k\in K}{\hat{F}}_k{E}_k(p),{\hat{F}}_k=\frac{1}{N^d}\sum _{p\in X}{F}_p{E}_{-k}(p),$$

Full text

{AMS}

65M70, 65R20, 76P05

1 Introduction

Gas kinetic theory describes the statistical behavior of a large number of gas molecules in the joint spatial and velocity space. It has been widely used to model gases outside the hydrodynamic regime, for example in the field of rarefied gas dynamics. Let $f(t,{x},{v})$ be the mass density distribution of the particles, depending on the time $t\in\mathbb{R}^{+}$, position ${x}\in\mathbb{R}^{d}$ ($d\geq 2$) and microscopic velocity ${v}\in\mathbb{R}^{d}$. Based on the molecular chaos assumption, the Boltzmann equation

[TABLE]

for the evolution of $f(t,{x},{v})$ was derived in [6] and has served as the fundamental equation in the gas kinetic theory. When modeling the binary interaction between the particles, the Boltzmann collision operator $\mathcal{Q}[f,f]$ takes the form

[TABLE]

for monatomic gases, where

[TABLE]

are the post-collisional velocities of two particles with pre-collisional velocities ${v}$ and ${v}_{*}$, and $\omega$ is the angular parameter of the collision. Here the variables $t$ and ${x}$ are omitted for simplicity and we shall continue doing so when focusing only on the collision term. The collision kernel $\mathcal{B}$ is non-negative and usually takes the form

[TABLE]

The Boltzmann equation (1) guarantees that $f(t,{x},{v})$ remains non-negative if the initial value $f(t=0,{x},{v})$ is also non-negative [28]. The symmetry of the collision term (2) and the fact $\,\mathrm{d}{v}\,\mathrm{d}{v}_{*}=\,\mathrm{d}{v}^{\prime}\,\mathrm{d}{v}^{\prime}_{*}$ imply that for any function $\psi(\cdot)$:

[TABLE]

Setting $\psi({v})=1,{v},|{v}|^{2}$ gives rise to the conservation of the mass, momentum and energy

[TABLE]

respectively. The famous H-theorem that states the monotonicity of the entropy

[TABLE]

can also be obtained by setting $\psi({v})=\ln(f({v}))$.

By introducing new variables ${y}={v}^{\prime}-{v}$ and ${z}={v}^{\prime}_{*}-{v}$ and carrying out algebra calculations, the Boltzmann collision operator can be rewritten as (see [7, 22, 19] for details)

[TABLE]

This is the well-known Carleman representation [7] of the Boltzmann collision operator, where $\tilde{\mathcal{B}}(y,z)$ is related to $\mathcal{B}({v}-{v}_{*},\omega)$ in (3) by

[TABLE]

Though the Boltzmann equation serves as the fundamental equation in gas dynamics, its high dimensional nature and the complexity of the collision operator pose difficulties for its numerical solution. A classical method is the direct Monte Carlo simulation [2], which uses simulated particles to mimic gas molecules and handles the collisions in a stochastic way. Though it treats the high dimensionality effectively, the convergence order is low and the numerical solution appears rather oscillatory.

With the rapid growth of computing power, it has become more practical to solve the Boltzmann equation with deterministic methods. For all deterministic approaches, the complexity of the collision integral poses the most serious difficulty for numerical computation. Therefore, this paper focuses on the spatially homogeneous case

[TABLE]

for simplicity.

In the past decades, several deterministic schemes have been developed for the Boltzmann collision term. Two methods that have attracted the most attention are the discrete velocity method (DVM) [13, 27, 4, 22] and the Fourier Galerkin method (FGM) [5, 23, 24], which will be reviewed in what follows.

1.1 Discrete velocity method

The discrete velocity method (DVM) assumes that the particle velocity takes only values from a finite set. Consider a domain $\mathcal{D}_{T}=[-T,T]^{d}$ for the velocity variable $v$, that is discretized uniformly with step size $h=2T/N$ for a positive integer $N$ (which is assumed to be odd for simplicity). By adopting the $d$-dimensional multi-index notation $k=(k_{1},\ldots,k_{d})$, one can denote the set of discrete velocity samples by

[TABLE]

where $N=2n+1$. In the rest of this paper, we use the lower-case letters $p,q,r$ and $s$ to denote the discrete velocity samples in ${X}$.

Using $F_{r}(t)$ for $r\in{X}$ as the numerical approximations of the distribution function $f(t,{v})$ at the points in ${X}$, the governing equations of DVM for $F_{r}(t)$ are

[TABLE]

Here $Q_{r}(t)$ serves as an approximation to $\mathcal{Q}[f,f](t,r)$. The coefficients $A_{pq}^{rs}$ are non-negative constants and satisfy the conservation relations

[TABLE]

and the symmetry property

[TABLE]

The property (12) shows that the collisions in DVM also satisfy the momentum and energy conservation, and the property (13) implies that the collisions are also reversible as in the Boltzmann equation. These two facts guarantee that DVM maintains a number of fundamental physical properties of the continuous Boltzmann equation, such as (a) the positivity of the distribution function, (b) the exact conservation of mass, momentum and energy, and (c) a discrete H-theorem.

More precisely, the values $F_{r}(t)$ for $r\in{X}$ are always non-negative if the initial values $F_{r}(t=0)$ are non-negative [28, 26]. The symmetry relation (13) implies that

[TABLE]

Combining (14) and the relations (12) gives rise to the conservation of mass, momentum and energy in the discrete sense

[TABLE]

After letting $\psi_{r}=\ln(F_{r})$ in (14), one obtains a discrete version of the H-theorem for DVM

[TABLE]

using the non-negativity of the coefficients $A_{pq}^{rs}$ and the monotonicity of the $\ln$ function. Notice that, in this argument, the symmetry relations (13), the non-negativity of $A_{pq}^{rs}$, and the non-negativity of the initial values are all essential for derivation of the H-theorem.

As a direct discretization for the high-dimensional integral of the collision term, DVM has a rather high computational cost $O(N^{2D+\delta})$ for some $0<\delta\leq 1$ [8]. It is also difficult to achieve decent convergence rate due to the insufficiency collision pairs in the Cartesian grid used for velocity discretization. More precisely, in the 2D case, the rate of convergence of DVM introduced in [13] is only $O((1/\log h)^{p})$ with $p<1/2-1/\pi$ [9]. For the 3D case, the best rate of convergence of DVM is also slower than the first order [21, 22].

The method in [18, 29] tries to improve the accuracy of DVM by interpolation. While the mass, momentum and energy are conserved in this scheme, positivity and the H-theorem fail to hold. The fast algorithm in [20] reduces the computational cost of DVM to $O(\bar{N}^{d}N^{d}\log(N))$ with some parameter $\bar{N}\ll N$ for the hard sphere molecules, but it abandons the conservation of momentum and energy.

1.2 Fourier Galerkin method

The Fourier-based methods assume that the distribution function $f(t,v)$ is supported (in the $v$ variable) in a ball $B_{R/2}$ centered at the origin with radius $R/2$. Under this assumption, it makes sense to focus on the functions $f(t,v)$ with $\mathrm{Supp}(f)\subset B_{R/2}$. For those functions, $\mathrm{Supp}(\mathcal{Q}[f,f])\subset B_{\sqrt{2}R/2}$ and the collision term $\mathcal{Q}[f,f]$ reduces to a truncated version $\mathcal{Q}^{R}[f,f]$ defined as

[TABLE]

where the superscript in $\mathcal{Q}^{R}[f,f]$ denotes the truncation radius. In order to obtain a spectral approximation to the collision term, one restricts the domain of the distribution function $f(v)$ to the cube $\mathcal{D}_{T}=[-T,T]^{d}$ with $T\geq\frac{3\sqrt{2}+1}{4}R$ in order to reduce aliasing. One then extends it periodically to the whole space. (See [24, 19] for details of the derivation). After periodization, $f({v})$ can be written as a Fourier series

[TABLE]

where $E_{k}({v})=\exp\left(\frac{\boldsymbol{\mathrm{i}}\pi}{T}k\cdot{v}\right)$. Substituting (18) into (17) gives rise to the following representation of the truncated collision operator:

[TABLE]

where

[TABLE]

It is easy to check that the coefficients $\hat{B}(l,m)$ are real and satisfy the symmetry relations

[TABLE]

In terms of the Fourier expansion,

[TABLE]

Here $\boldsymbol{\mathrm{1}}(\cdot)$ is the indicator function, equal to $1$ at the origin and [math] otherwise.

The Fourier Galerkin method (FGM) in the literature (e.g. [24]) starts with a finite square grid of Fourier modes

[TABLE]

and the subspace

[TABLE]

We shall denote the grid points in $K$ with lower-case letters $j,k,l,m$. FGM approximates the collision term (22) by projecting it to the subspace $\mathbb{P}_{N}$

[TABLE]

where $\hat{F}_{k}$ for $k\in K$ serve as the approximation of the Fourier modes $\hat{f}_{k}$ of the exact solution.

Putting together the above discussion, one arrives at the equations of the discrete Fourier coefficients $\hat{F}_{k}(t)$ for $k\in K$ for FGM [24]

[TABLE]

where $\hat{F}_{k}^{0}$ are the Fourier coefficients of the initial condition $f^{0}({v})$ restricted on $\mathcal{D}_{T}$.

Remark 1.1**.**

The above description of the Fourier Galerkin method is based on the Carleman representation of the Boltzmann collision operator (7). Starting from the classical form (2), one can also derive a relation similar to (25) (see [24] for details), while the definition of $\hat{B}(l,m)$ is slightly different:

[TABLE]

*where $T\geq\frac{3+\sqrt{2}}{4}R$ in the definition of $E_{k}(\cdot)$. It is straightforward to check that these coefficients also satisfy the symmetry relation (21). *

FGM achieves spectral accuracy, although the computational cost is still as high as $O(N^{2d})$ [24]. Two fast algorithms [19, 11] reduced the cost to $O(MN^{d}\log(N))$ for the hard sphere molecules (the Maxwell molecules for 2D case) [19] and to $O(MN^{d+1}\log(N))$ for general collision kernels [11], where $M$ is the number of points used for discretizing the unit sphere $\mathbb{S}^{d-1}$.

Compared to DVM, the solution of FGM loses most of the aforementioned physical properties, including positivity, the conservation of momentum and energy, and the H-theorem. In [25], Pareschi and Russo proposed a positivity preserving regularization of FGM by using Fejér filter at the expense of spectral accuracy. Despite this, the solution fails to satisfy the H-theorem. The loss of the conservation can be fixed by a spectral-Lagrangian strategy [12].

1.3 Motivation

DVM preserves a number of physical properties (such as positivity of the solution, the H-theorem, and exact conservation of mass, momentum and energy) but suffers from high computational costs and low accuracies. FGM enjoys spectral accuracies and lower computational costs but sacrifices almost all physical properties except the mass conservation. In this paper, we aim for a trade-off between the physical properties and the spectral accuracy.

As a fundamental property of the solution to the Boltzmann equation, the positivity of the distribution function helps establish the H-theorem, which is one of the crucial properties in order to guarantee the well-posedness of the discrete system. Therefore, it makes sense to maintain the positivity and the H-theorem, as long as it does not sacrifice significantly other properties such as numerical accuracy and efficiency. This paper is an initial study in this direction.

To achieve this goal, we first study carefully the reason behind the loss of the H-theorem in FGM by comparing it with DVM. With a few novel modifications to FGM, we propose an entropic Fourier method (EFM) that preserves the positivity, the mass conservation, and the H-theorem. In addition, the computational cost of this new method is the same as the one of FGM.

The rest of the paper is organized as follows. In Section 2, we first outline the key steps to develop EFM and state the main results of the paper. The details of the derivation and some deeper understandings of the model are provided in Section 3. Section 4 presents the implementation of EFM and the numerical results. The paper ends with a discussion in Section 5.

2 Main result

This section outlines the overall procedure of our derivation and lists some key results. Detailed derivation and investigation will be given in Section 3.

Aiming at developing an entropic Fourier method for the homogeneous Boltzmann equation, one works mainly with the evolution of the Fourier coefficients $\hat{F}_{k}(t)$. Recall the discrete Fourier transform

[TABLE]

where ${X}$ defined in (10) and $K$ defined in (23) are the sets of uniform samples in the velocity space and the Fourier domain, respectively.

Using the discrete Fourier transform, one can instead treat the point values $F_{p}$ as the degrees of freedom and write the numerical scheme in the DVM form (11). According to the derivation in Section 1.1, the following condition is required in order to guarantee the H-theorem for DVM:

Condition 1**.**

The DVM defined in (11) satisfies

the coefficients $A_{pq}^{rs}$ satisfy the symmetry relation $A_{pq}^{rs}=A_{qp}^{rs}=A_{rs}^{pq}$; 2. 2.

the coefficients $A_{pq}^{rs}$ are non-negative, i.e $A_{pq}^{rs}\geq 0$; 3. 3.

the initial values are non-negative, i.e $F_{p}(t=0)\geq 0$ for any $p\in{X}$.

The general idea of our approach is to revise the existing Fourier Galerkin method (FGM) so that Condition 1 is fulfilled. Below we list the steps that lead to a numerical scheme that satisfies the H-theorem.

Apply the Fourier collocation method to (19). This leads to an approximation to (22) in the form

[TABLE]

where $\boldsymbol{\mathrm{1}}_{N}(l):=\boldsymbol{\mathrm{1}}(l\bmod N)$ and $\hat{B}_{N}(l,m):=\hat{B}(l\bmod N,m\bmod N)$. Here $\bmod$ is the symmetric modulo function, i.e. each component of $l\bmod N$ ranging from $-n$ to $n$ (recall $N=2n+1$). Using the relation between the Fourier coefficients and the values on collocation points (28), we can rewrite (29) as

[TABLE]

where $A_{pq}^{rs}$ (given in (45)) is determined by the Fourier modes of the collision kernel $\hat{B}_{N}(\cdot,\cdot)$ and satisfies the symmetry relation (Condition 1.1). 2. 2.

A careful study shows that $A_{pq}^{rs}$ fails to be non-negative. This can be fixed by applying a positivity preserving filter to $\hat{B}_{N}(l,m)$, i.e.,

[TABLE]

where ${\sigma}_{N}(l)$ is the tensor-product of $d$ one-dimensional modified Jackson filter [17, 30]. The modified collision term takes the following form in the Fourier domain

[TABLE]

Using ${(A^{\sigma})}^{rs}_{pq}$ to denote the coefficients determined by the new kernel modes $\hat{B}_{N}^{\sigma}(l,m)$ and writing

[TABLE]

one can verify that both the symmetry relation (Condition 1.1) and the non-negativity (Condition 1.2) are satisfied. 3. 3.

To guarantee the positivity of the initial values (Condition 1.3), we adopt interpolation rather than orthogonal projection while discretizing the initial distribution function.

Main result.

Summarizing the outline given above, we arrive at a new entropic Fourier method (EFM) that takes the following simple form

[TABLE]

This method preserves several key physical properties, as guaranteed by the following theorem.

Theorem 2.1**.**

If $f(t=0,{v})\geq 0$ for ${v}\in\mathbb{R}^{d}$, then the solution $F_{r}(t)=\sum_{k\in{K}}\hat{F}_{k}(t)E_{k}(r)$ for $r\in{X}$ of (32) satisfies for $\forall t>0$,

[TABLE]

The proof is presented in Section 3.5. Due to the positivity-preserving filter (58), the numerical accuracy of EFM in approximating the collision operator is second order (see Section 3.4 for details).

Another important result for EFM is the existence of fast algorithms. For FGM, the fast algorithms proposed in [19, 11] are based on the approximation of the kernel $\hat{B}_{N}(\cdot,\cdot)$:

[TABLE]

with the number of terms $M\ll N^{d}$. Since the filtered kernel $\hat{B}_{N}^{\sigma}(l,m)$ turns out to have a similar approximation

[TABLE]

these fast algorithms still apply. Moreover, when the above approximation is applied, the H-theorem still holds. Detailed discussion will be given in Section 4.1.

3 Entropic Fourier method

As shown in Section 1.1, a discrete H-theorem can be obtained from the classical DVM, where the associated entropy function can be considered as a numerical quadrature for the integral of $f\ln f$. This requires the positivity of the distribution function, which can be guaranteed by the positivity of the discrete collision kernel $A_{pq}^{rs}$. In general, to preserve the Boltzmann entropy in the numerical scheme, the positivity of the numerical solution needs to be enforced in a certain sense due to the presence of $\ln f$ in the entropy function. However, in FGM, there is no guarantee of any form of positivity in the numerical solution, and hence the H-theorem does not hold.

In this paper, rather than enforcing the non-negativity of the whole distribution function, we take a collocation approach and focus on the non-negativity only at the collocation points. Based on this idea, we start from a collocation method for the homogeneous Boltzmann equation and write it as a DVM of the function values defined at the collocation points. One then tries to alter the coefficients to match the requirements in Condition 1 so that the H-theorem can be subsequently derived.

The three steps listed in Section 2 are elaborated in the first three subsections below. After that, Section 3.4 compares the entropic Fourier method (EFM) with other Fourier methods.

3.1 Fourier collocation method in a DVM form

The mechanisms of DVM and FGM are quite different: DVM is concerned with the values of the distribution on discrete points, whereas FGM (25), as a Galerkin method, works on the Fourier modes of the distribution function. It is not straightforward how to link these two methods. Alternatively, we will consider another type of Fourier methods — the collocation method (also known as the pseudospectral method).

3.1.1 Fourier collocation method

In the Fourier collocation method (FCM), the collision term on the set $X$ is evaluated directly using (19):

[TABLE]

where $\hat{B}_{N}(l,m):=\hat{B}(l\bmod N,m\bmod N)$, and $\bmod$ is the symmetric modulo function, i.e. each component of $l\bmod N$ ranging from $-n$ to $n$ (recall $N=2n+1$). Since the above equation only uses the value of $\hat{B}_{N}$ in $K$, one can use $\hat{B}$ and $\hat{B}_{N}$ interchangeably here. Here we note that $\hat{B}_{N}(l,m)$ satisfy the symmetry relation (21).

The corresponding Fourier modes can be obtained by an inverse discrete Fourier transform:

[TABLE]

where $\boldsymbol{\mathrm{1}}_{N}(l):=\boldsymbol{\mathrm{1}}(l\bmod N)$.

If the initial value is smooth enough, due to the smoothing effect of the Boltzmann collision operator [1], both the Fourier Galerkin and the Fourier collocation methods have spectral accuracy [8]. Moreover, in some cases, FCM (39) is numerically more efficient, especially for the fast summation algorithms in [19, 11]. For example, in [19], the following approximation of $\hat{B}_{N}(l,m)$ is considered:

[TABLE]

where $M\in\mathbb{N}^{+}$ is the total number of quadrature points on the sphere. Then the collision term in this Galerkin method can be approximated by

[TABLE]

To evaluate (41) efficiently, one needs to utilize FFT-based convolutions. To obtain these coefficients, one needs the zero-padding technique to avoid aliasing. If one uses the same method to evaluate $\hat{Q}_{k}$ in (39), then no zero padding is needed. Therefore the collocation method shortens the length of vectors used in the Fourier transform, which makes the algorithm faster.

3.1.2 DVM form

To link FCM with DVM, we split the collision term (38) into the gain part ($+$) and the loss part ($-$):

[TABLE]

Noticing $E_{l+m}(r)=\sum_{k\in{K}}\boldsymbol{\mathrm{1}}_{N}(l+m-k)E_{k}(r)$ and plugging (28) into the gain part yields

[TABLE]

Since $\frac{1}{N^{d}}\sum_{s\in{X}}E_{j}(s)=\boldsymbol{\mathrm{1}}_{N}(j)$, one can sum over $j$ to get $\frac{1}{N^{d}}\sum_{j\in K,s\in{X}}E_{j}(s)=1$. With this equation, one introduces two new indices to (43) by multiplying its right hand side with $\frac{1}{N^{d}}\sum_{j\in K,s\in X}E_{j}(s)$:

[TABLE]

If one introduces

[TABLE]

then the gain term is

[TABLE]

Apparently, such term does take the form of the gain term of DVM (11).

For the loss term, the identity

[TABLE]

leads to the following derivation

[TABLE]

In summary, one can write FCM in the following DVM form

[TABLE]

with $A_{pq}^{rs}$ given in (45). Finally, the symmetry relation (Condition 1.1)

[TABLE]

holds, as this can be easily seen by the symmetry relation of $\hat{B}_{N}(l,m)$ and switching the indices in (45).

3.2 Positivity preserving

As remarked earlier, in order to obtain an H-theorem for FCM, one needs to ensure that all the coefficients $A_{pq}^{rs}$ are non-negative (Condition 1.2). Below, we first show that $A_{pq}^{rs}$ as defined in (49) fail to be non-negative, and then apply a filter to recover non-negativity.

3.2.1 Failure of positivity preservation in FCM

We start by simplifying the coefficients $A_{pq}^{rs}$ based on (45)

[TABLE]

where one uses $j=l+m-k\bmod N$. By performing a change of variables $i=(m-k)\bmod N$ and $j=(l-k)\bmod N$, we arrive at

[TABLE]

By introducing

[TABLE]

which is by definition a periodic function with period $\mathcal{D}_{T}$, one can write compactly

[TABLE]

In order to check whether $A_{pq}^{rs}$ is non-negative, one just needs to check whether $G(\cdot,\cdot)$ is non-negative on the collocation points in ${X}$. To get a better understanding of the function $G(\cdot,\cdot)$ as defined in (53), one applies the definition of $\hat{B}_{N}(\cdot,\cdot)$ to obtain

[TABLE]

Here $\chi_{N}$ is the Dirichlet kernel over $\mathcal{D}_{T}$ defined by

[TABLE]

and its discrete Fourier transform $\hat{\chi}_{N}(k)$ is equal to 1 for $k\in{K}$ and 0 on $\mathbb{Z}^{d}\setminus{K}$. By introducing a periodic function in $\mathcal{D}_{T}$

[TABLE]

one can write

[TABLE]

where the convolution is defined periodically in $\mathcal{D}_{T}\times\mathcal{D}_{T}$. Equivalently, $G({y},{z})$ is also the truncated Fourier expansion of $H({y},{z})$ by keeping only the frequencies in ${K}$.

Although $H({y},{z})$ is non-negative in the weak sense, its truncated Fourier approximation $G({y},{z})$ fails to be so. For example, the values of $G$ for the kernel $\tilde{\mathcal{B}}({y},{z})\equiv\frac{1}{\pi}$, $R=6$ in 2D are plotted in Figure 1. This clearly shows that negative values appear as expected. Therefore in general, the H-theorem does not hold for FCM.

3.2.2 Filtering

In the previous subsection, one can see that if $\chi_{N}(\cdot)$ were a non-negative function, then $G({y},{z})$ would be non-negative for any ${y}$ and ${z}$. Thus, in order to get non-negative coefficients, a possible way is to replace the function $\chi_{N}$ by a non-negative one. Note that $\chi_{N}({v})$ is Dirichlet kernel, which is an approximation of the Dirac delta. As $N\rightarrow+\infty$, the function $\chi_{N}(\cdot)$ tends to Dirac delta weakly in an oscillatory way. As pointed out in [10], the oscillation breaks the non-negativity of the solution.

As mentioned earlier in Section 2, we adopt the one-dimensional modified Jackson filter [17, 30] given by

[TABLE]

where $N=2n+1$ and $-n\leq\beta\leq n$. By a slight abuse of notation, the $d$-dimensional modified Jackson filter for a multiindex $k=(k_{1},\ldots,k_{d})\in{K}$ is defined through tensor-product

[TABLE]

The modified kernel $\chi^{\sigma}_{N}(v)$ can then be defined as

[TABLE]

Once $\chi_{N}$ is replaced with $\chi_{N}^{\sigma}$ in (55), the function $G({y},{z})$ is substituted with $G^{\sigma}({y},{z}):=(G*(\chi^{\sigma}_{N}\otimes\chi^{\sigma}_{N}))({y},{z})$. A direct calculation shows that

[TABLE]

where $\hat{B}_{N}^{\sigma}(l,m):=\hat{B}_{N}(l,m)\sigma_{N}(l)\sigma_{N}(m)$ as defined in (31). With $G^{\sigma}({y},{z})\geq 0$ guaranteed, one can mimic (54) and define

[TABLE]

which are apparently non-negative. Since replacing $\hat{B}_{N}(l,m)$ with $\hat{B}_{N}^{\sigma}(l,m)$ does not affect the symmetry relation (21), the new coefficients $(A^{\sigma})^{rs}_{pq}$ also satisfies $(A^{\sigma})^{rs}_{pq}=(A^{\sigma})^{rs}_{qp}=(A^{\sigma})^{pq}_{rs}$.

From the above discussion, we now define the entropic collision term to be

[TABLE]

in the Fourier domain. In the velocity domain, it is equal to

[TABLE]

3.3 Initial condition

In order to obtain the H-theorem, one needs to make sure that the initial data are non-negative. Since the collision term in (62) depends only on the values at the collocation points in $X$, it is sufficient to make the initial data non-negative at these collocation points. Consequently, it is natural to use sampling rather than orthogonal projection while preparing the discrete initial data $F^{0}_{r}$ for $r\in{X}$. More precisely,

[TABLE]

where $\varphi^{\epsilon}\geq 0$ is a mollifier, such that $\|f^{0}-\varphi^{\epsilon}*f^{0}\|_{L^{2}}<\epsilon$ for $\epsilon$ sufficiently small. Once $\{F^{0}_{r}\}$ are ready, the corresponding Fourier coefficients $\{\hat{F}^{0}_{k}\}$ are computed via a fast Fourier transform.

At this point, all ingredients of the entropic Fourier method (EFM) are ready. The Cauchy problem of EFM takes the following form as a DVM

[TABLE]

Equivalently in the Fourier domain, EFM takes the form

[TABLE]

Remark 3.1**.**

The same technique can be applied to the Fourier method derived from the classical form of the Boltzmann collision operator (2) to obtain an entropic Fourier method. In fact, from (27) and following the definition of $G({y},{z})$ in (53), one can directly obtain

[TABLE]

*Again, by (54), the positivity of $A_{pq}^{rs}$ depends only on the positivity of $G({y},{z})$ at the collocation points. Therefore, replacing $\chi_{N}$ with $\chi_{N}^{\sigma}$ does the job. *

3.4 Comparison

The derivation of EFM switched frequently between the language of Fourier methods and DVM for different purposes. As we have shown in (64) and (65), EFM can be regarded either as a Fourier method or as a special DVM.

In what follows, we provide a comparison between the entropic Fourier method (EFM), the FGM (Fourier Galerkin method), and the Fourier collocation method (FCM). To set up a uniform notation, let $\mathcal{Q}[\cdot;\cdot,\cdot]$ be the general collision operator

[TABLE]

with a collision kernel $C(\cdot,\cdot)$. Thus the truncated collision term (17) can be written as $\mathcal{Q}(H;f,f)$ using the definition of $H$ in (57).

Notice that a special feature of a function in $\mathbb{P}_{N}$ is that it is uniquely defined via its function values at points in $X$ defined in (10). Therefore, $\{F_{p}|p\in{X}\}$ can be regarded both as discrete set of values or the samples from the smooth periodic $F(v)\in\mathbb{P}_{N}\subset L_{\mathrm{per}}^{2}(\mathcal{D}_{T})$. By introducing two operators

[TABLE]

for the space $L_{\mathrm{per}}^{2}(\mathcal{D}_{T})$, the three methods are different approximations of $\mathcal{Q}(H;f,f)$ with different initial values:

[TABLE]

where $\mathcal{I}_{N}$ is the interpolation operator defined in (63).

The list (67) to (69) clearly shows how we change from FGM to EFM in our derivation. The last line (69) also shows that EFM provides an approximation of the original binary collision operator in the language of spectral methods. Below we will briefly review the basic properties of all the three methods.

The method (67) stands for FGM as described in Section 1.2. In the derivation, the kernel $K$ is not explicitly projected. However, (25) shows that the discrete collision operator depends only on $(\mathcal{P}_{N}\times\mathcal{P}_{N})H$. By replacing the projection operator applied to $\mathcal{Q}$ with interpolation as in (39), we arrive at the Fourier collocation method (68) introduced in Section 3.1. Since a direct projection of $H$ does not preserve the positivity of the kernel, the negative part of the discrete kernel may cause a violation of the H-theorem. Nevertheless, both these two methods have spectral accuracy in the velocity space.

To ensure the positivity of the discrete kernel, the filter $\mathcal{S}_{N}^{\sigma}\otimes\mathcal{S}_{N}^{\sigma}$ is applied in (69), and thus positive coefficients (60) are obtained. The method (69) also ensures the positivity of the approximation of $F$ at collocation points, and thus the discrete H-theorem follows.

However, the filter $S_{N}^{\sigma}$ has a smearing effect, which reduces the order of convergence. For any smooth periodic function $f\in L_{\mathrm{per}}^{2}(\mathcal{D}_{T})$, the $L^{2}$-error $\|f-\mathcal{S}_{N}^{\sigma}f\|_{2}$ is $O(N^{-2})$ [16, Chapter 4], and therefore EFM is at most second-order. On the other hand, if one splits the collision term of EFM into the gain part and the loss part and let $F=\mathcal{P}_{N}f$, then

[TABLE]

for any $r\in{X}$. Following the boundedness of the truncated collision operator proven in [24], one concludes that

[TABLE]

Hence, EFM has second-order accuracy in approximating the truncated collision operator, if the distribution is smooth enough (due to the smoothing effect of the Boltzmann collision operator [1], we only the initial value is smooth enough). This order of convergence will be also numerically verified in the next section.

3.5 Proof of Theorem 2.1

This subsection provides the proof of Theorem 2.1.

Proof 3.2** (Proof of Theorem 2.1).**

The argument in Section 3.2.2 indicates that if $f(t=0,{v})\geq 0$, then the coefficients $A_{pq}^{rs}$ and the initial values $F^{0}_{r}$ satisfy all the three conditions in Condition 1.

The symmetry relation (21) of $\hat{B}(l,m)$ and the definition of $\hat{B}_{N}^{\sigma}(l,m)$ indicates $\boldsymbol{\mathrm{1}}_{N}(l+m)\left(\hat{B}_{N}^{\sigma}(l,m)-\hat{B}_{N}^{\sigma}(m,m)\right)=0$, i.e. $Q^{\sigma}_{0}=0$. Noticing that the zero frequency $\hat{F}_{0}(t)=\frac{1}{N^{d}}\sum_{r\in X}F_{r}(t)$, one can directly obtain the conservation of mass (33).

Since $f^{0}(v)\geq 0$, $F^{0}_{r}\geq 0$ by construction. If there exists $t^{\prime}>0$ and $r\in{X}$ such that $F_{r}(t^{\prime})=0$ and $F_{p}(t^{\prime})\geq 0$ for any other $p\in{X}$, then

[TABLE]

which indicates $F_{r}(t)\geq 0$ for all $t>0$ and $r\in{X}$.

The symmetry relation and non-negativity of $(A^{\sigma})_{pq}^{rs}$ indicate the discrete H-theorem (35)

[TABLE]

*This completes the proof. *

4 Numerical tests

This section describes several numerical tests to demonstrate the properties of EFM and to compare with the Fourier Galerkin method (FGM) in [24] and the positivity preserving spectral method (PPSM) in [25].

4.1 Implementation

It is pointed out in Section 2 that the fast algorithms in [19, 11] can be applied to EFM without affecting the H-theorem. To show this, one needs to check that the fast algorithms do not violate the first two conditions in Condition 1.

These fast algorithms are based on an approximation of $\hat{B}_{N}(l,m)$ (37) of the following form.

[TABLE]

After the filtering is applied, one obtains a similar approximation for $\hat{B}_{N}^{\sigma}(l,m)$

[TABLE]

It can be verified that both kernels satisfy the symmetry relation

[TABLE]

which indicates Condition 1.1 is valid for the fast algorithms.

To see that the fast algorithms do not affect the non-negativity of $G^{\sigma}({y},{z})$, we use the fast algorithm in [19] with $d=2$ and $\tilde{B}=1$ as an example. The first step of this algorithm writes ${y}$ and ${z}$ in (20) in the polar coordinates ${y}=\rho{e}$ and ${z}=\rho_{*}{e}_{*}$:

[TABLE]

Let $\psi_{R}(l,{e})=\int_{-R}^{R}E_{l}(\rho{e})\,\mathrm{d}\rho$, then

[TABLE]

Integrating it with respect to ${e}_{*}$ yields

[TABLE]

Substituting (76) into (59) gives rise to

[TABLE]

Let $\phi^{\sigma}_{R}({y},{e})=\int_{-R}^{R}\chi_{N}^{\sigma}(\rho{e}-{y})\,\mathrm{d}\rho$. Apparently, $\phi^{\sigma}_{R}({y},{e})\geq 0$ due to $\chi_{N}^{\sigma}({y})\geq 0$ for any ${y}\in\mathbb{R}^{2}$. Then integrating (78) with respect to ${e}_{*}$ yields

[TABLE]

The idea of the fast algorithm is to replace the integration in (77) with a quadrature formula. More precisely (77) is approximated by

[TABLE]

Similarly to (80), one obtains

[TABLE]

Since $\phi^{\sigma}_{R}({y},{e})\geq 0$ for any ${y}\in\mathbb{R}^{2}$, ${e}\in\mathbb{S}^{1}$, $G^{\sigma}_{\mathsf{fast}}({y},{z})\geq 0$ for any ${y},{z}\in\mathbb{R}^{2}$. Hence, the fast algorithm does not destroy the non-negativity of $G^{\sigma}_{\mathsf{fast}}({y},{z})$.

As we pointed out in Section 3.1, an aliased convolution can be directly used to calculate (32). Since the accuracy of EFM is only second order, the smoothing filter is the main source of the error. In the fast algorithms, the number $M$ in (72) perhaps can be smaller than that in [19, 11].

In the above discussion, we only study the case when $N$ is odd. The case of even $N$ values can be reduced to the odd $(N-1)$ case by setting the coefficient of a mode $k$ to be zero [math] if any component of $k=(k_{1},\ldots,k_{d})$ is equal to $-N/2$.

For the time discretization, the third-order strong stability-preserving Runge-Kutta method proposed in [14] is employed in the discretization of time. In all tests, the time step is chosen as $\Delta t=0.01$.

4.2 Numerical results

The test problems used here are solutions of the space-homogeneous Boltzmann equation for Maxwell molecules ($\mathcal{B}({g},\omega)=\frac{1}{2\pi}$ in 2D and $\mathcal{B}({g},\omega)=\frac{1}{4\pi}$ in 3D).

Example 1 ($2D$ BKW solution).

The first example is the well-known $2D$ BKW solution, obtained independently in [3] and [15]. The exact solution takes the form

[TABLE]

where $S=1-\exp(-t/8)/2$. The BKW solution allows one to check the accuracy, positivity of the solution, and the entropy of the proposed method. Here we set the truncation radius $R=6$ in the tests.

Both FGM (26) or FCM (39) result in good approximations of the exact solution at the collocation points. However, the solutions are not non-negative. Even if we use interpolation rather than orthogonal projection to prepare the initial values, the solutions of these methods still fail to be non-negative. The positivity of the solution is measured by the positivity error defined via

[TABLE]

Figure 2 shows that both FGM and FCM fail to preserve the positivity of the solution at the collocation points no matter whether the initial value is given by orthogonal projection or interpolation. On the other hand, the positivity error of EFM is strictly zero, thanks to the modification in (59).

Table 1 summarizes the $\ell_{1}$, $\ell_{2}$ and $\ell_{\infty}$ errors of (64) at time $t=0.01$. Here the $\ell_{\mathrm{p}}$ relative errors for $\mathrm{p}=1,2,\infty$ are defined by

[TABLE]

where $f(q)$ is the exact solution at $q\in{X}$. The numerical results also show that the convergence rate of EFM is of second order.

As discussed in Section 4.1, the fast algorithm in [19] can be applied to EFM to accelerate the computation. In (77), the integration on $[0,\pi)$ can be reduced to $[0,\pi/2)$ and $M$ is equal to the number of samples within $[0,\pi/2)$. Table 2 presents the $\ell_{1}$ error for multiple values of $M$, noticing that $M=2$ is good enough in practice while in [19] the authors suggest $M\geq 4$.

In order to demonstrate that the proposed method satisfies the H-theorem numerically, we define a time-dependent discrete entropy function

[TABLE]

The evolution of the entropy, plotted in Figure 3, shows that as the number of the discrete points $N$ increases the discrete entropy converges to the one of the exact solution.

As a comparison with PPSM, Figure 4 presents the numerical solutions in the $v_{1}$ direction of PPSM and EFM at $t=1$ with $N=32$. The smoothing filter used for EFM results in much less dissipation, thus leading to a much better agreement with the exact solution. Finally, Figure 5 shows that as $N$ increases the solution of EFM converges rapidly to the exact solution.

Example 2 (Bi-Gaussian initial value).

Another frequent example is a problem with the bi-Gaussian initial value

[TABLE]

where ${u}_{1}=(-2,0)^{\mathsf{T}}$ and ${u}_{2}=(2,0)^{\mathsf{T}}$. This is solved for the Maxwell molecules (2D in velocity) with radius $R=6$. Figure 6 shows the numerical results of PPSM and EFM. The reference solution is calculated by the Fourier spectral method with $N=400$ and $R=8$. It is clear that EFM solution is much closer to the reference solution. Figure 7 demonstrates that as $N$ increases EFM converges rapidly to the reference solution.

Example 3 (Discontinuous initial value).

The initial condition given by

[TABLE]

in this example is discontinuous. Here $\rho_{1}=\frac{6}{5}$ and the values of $\rho_{2}$, $T_{1}$ and $T_{2}$ are uniquely determined by the three conditions

[TABLE]

The profile of the reference solution is presented in Figure 8, which is computed by EFM with $N=2048$. Due to the discontinuity in the initial value, the spectral accuracy of FGM is lost. In addition, the Gibbs phenomenon leads to oscillations in the initial value of FGM. In Figure 9, the plots around the discontinuity demonstrate that EFM has a much better agreement as compared to FGM. The oscillations in FGM solutions exhibit large errors and the amplitude of the oscillation decreases slowly as $N$ increases. On the contrary, there is no oscillation for EFM around the discontinuity and the solution is always kept non-negative.

Example 4 ($3D$ BKW solution).

The solution of this example is the exact $3D$ BKW solution, given by

[TABLE]

where $S=1-2\exp(-t/6)/5$. Similarly to the $2D$ case, we first check the accuracy of EFM. At time $t=0.01$, the $\ell_{1}$, $\ell_{2}$ and $\ell_{\infty}$ errors and the convergence rates are listed in the Table 3. Similar to the $2D$ case, the convergence rate is of the second order and the errors are rather small.

As a comparison with PPSM, Figure 10 presents the numerical solutions on the $v_{1}$ direction of PPSM and EFM at $t=1$ with $N=32$. The plots clearly show that the smoothing filter used in EFM results in much less dissipation, thus leading to a better agreement with the exact solution.

5 Discussion

EFM proposed in this paper is a trade-off between accuracy and preservation of physical properties. The resulting scheme can be viewed both as a discrete velocity method and a Fourier method. In terms of the convergence rate, it is better than DVM but slower than FGM. In terms of physical properties, it guarantees positivity, mass conservation and a discrete H-theorem, while the momentum and energy conservation is lost. Regarding the computational cost, fast algorithms in [19, 11] remain valid for EFM. As to the future work, we plan to study how to mitigate momentum and energy loss, where higher order accuracy is needed for long time simulation. The numerical implementation of the spatially inhomogeneous setting is also in progress.

Acknowledgments

The authors thank Jingwei Hu for discussion on the filter and fast algorithm on Boltzmann collision term.