A bi-fidelity method for the multiscale Boltzmann equation with random parameters

TL;DR

This paper introduces a bi-fidelity stochastic collocation method for efficiently solving the multiscale Boltzmann equation with random parameters, combining low- and high-fidelity models to accurately approximate macroscopic quantities.

Contribution

The paper develops a bi-fidelity approach that leverages a low-fidelity Euler model and a few high-fidelity Boltzmann simulations for efficient uncertainty quantification.

Findings

The bi-fidelity method accurately captures macroscopic quantities.

The approach reduces computational cost significantly.

Numerical experiments confirm high efficiency and accuracy.

Abstract

In this paper, we study the multiscale Boltzmann equation with multi-dimensional random parameters by a bi-fidelity stochastic collocation (SC) method developed in [A. Narayan, C. Gittelson and D. Xiu, SIAM J. Sci. Comput., 36 (2014); X. Zhu, A. Narayan and D. Xiu, SIAM J. Uncertain. Quantif., 2 (2014)]. By choosing the compressible Euler system as the low-fidelity model, we adapt the bi-fidelity SC method to combine computational efficiency of the low-fidelity model with high accuracy of the high-fidelity (Boltzmann) model. With only a small number of high-fidelity asymptotic-preserving solver runs for the Boltzmann equation, the bi-fidelity approximation can capture well the macroscopic quantities of the solution to the Boltzmann equation in the random space. A priori estimate on the accuracy between the high- and bi-fidelity solutions together with a convergence analysis is established. Finally, we present extensive numerical experiments to verify the efficiency and accuracy of our proposed method.

Figures

Equations

$${\partial }_{t}f+v\cdot {\nabla }_{x}f=\frac{1}{\epsilon }\mathcal{Q}(f,f),(x,v)\in \Omega \times {\mathbb{R}}^{d_v},$$

$$\mathcal{Q}(f,f)(v)={\int }_{{\mathbb{R}}^{d_v}}{\int }_{{\mathbb{S}}^{d_v-1}}B(|v-v_{\ast }|,\cos \theta )\left(f(v^{\prime })f(v_{\ast }^{\prime })-f(v)f(v_{\ast })\right)d\sigma d{v}_{\ast }.$$

$$\{\begin{array}{ll}& v^{\prime }=\frac{v+v_{\ast }}{2}+\frac{|v-v_{\ast }|}{2}\sigma ,\\ & v_{\ast }^{\prime }=\frac{v+v_{\ast }}{2}-\frac{|v-v_{\ast }|}{2}\sigma ,\end{array}$$

$$B(v,v_{\ast },\sigma )=b|v-v_{\ast }{|}^{\gamma },-d_v<\gamma \le 1,$$

$$m(v)={\left(1,v,\frac{|v{|}^2}{2}\right)}^{\top },$$

$${\int }_{{\mathbb{R}}^{d_v}}\mathcal{Q}(f,f)m(v)dv=0,$$

$${\int }_{{\mathbb{R}}^d}\mathcal{Q}(f,f)\ln fdv\le 0.$$

$$M(v{)}_{\rho ,u,T}=\frac{\rho }{(2\pi T{)}^{\frac{d_v}{2}}}\exp \left(-\frac{|v-u{|}^2}{2T}\right),$$

$$\rho ={\int }_{{\mathbb{R}}^{d_v}}f(v)dv,u=\frac{1}{\rho }{\int }_{{\mathbb{R}}^{d_v}}f(v)vdv,T=\frac{1}{d_v\rho }{\int }_{{\mathbb{R}}^{d_v}}f(v)|v-u{|}^{2}dv.$$

$$\displaystyle\left\{\begin{array}[]{l}\displaystyle\partial_{t}f+v\cdot\nabla_{x}f=\frac{1}{\varepsilon}\mathcal{Q}_{z}(f,f),\\[4.0pt] \displaystyle f(0,x,v,z)=f_{I}(x,v,z),\qquad(x,v)\in\Omega\times\mathbb{R}^{d_{v}},z\in I_{z}.\end{array}\right.$$

$${\partial }_t\left(\begin{array}{c}\rho \\ \rho u\\ E\end{array}\right)+{\nabla }_x\cdot \left(\begin{array}{c}\rho u\\ \rho u\otimes u+\mathbb{P}\\ Eu+\mathbb{P}u+\mathbb{Q}\end{array}\right)=0,$$

$$E=⟨\frac{1}{2}|v{|}^{2}f⟩=\frac{1}{2}\rho |u{|}^2+\frac{d_v}{2}\rho T,$$

$$⟨g⟩={\int }_{{\mathbb{R}}^{d_v}}g(v)dv.$$

$$\mathbb{P}=p\mathbb{I},\mathbb{Q}=0,$$

$${\partial }_t\left(\begin{array}{c}\rho \\ \rho u\\ E\end{array}\right)+{\nabla }_x\cdot \left(\begin{array}{c}\rho u\\ \rho u\otimes u+p\text{I}\\ (E+p)u\end{array}\right)=0,$$

$$U^L(\gamma )=\text{span}\{u^L(\gamma )\}=\text{span}\{u^L(z_1),\cdots ,u^L(z_k)\}.$$

$$u^H(\gamma )=[u^H(z_1),\cdots ,u^H(z_k)],U^H(\gamma )=\text{span}\{u^H(\gamma )\}.$$

$$u^L(z)\approx {\mathcal{P}}_{U^L({\gamma }_N)}{u}^L=\sum _{k=1}^N{c}_k(z)u^L(z_k).$$

$$u^B(z)=\sum _{k=1}^N{c}_k(z)u^H(z_k).$$

$$z_{i_k}=\arg \underset{z\in \Gamma }{\max }\text{dist}(u^L(z),U^L({\gamma }_{k-1})),{\gamma }_k={\gamma }_{k-1}\cup z_{i_k},$$

$$u^L(z)\approx {\mathcal{P}}_{U^L({\gamma }_N)}[u^L(z)]=\sum _{k=1}^N{c}_k(z)u^L(z_{i_k}),$$

$${\mathbf{G}}^L\mathbf{c}=\mathbf{f},\mathbf{f}=(f_k{)}_{1\le k\le N},f_k=⟨u^L(z),u^L(z_{i_k}){⟩}^L,$$

$$({\mathbf{G}}^L{)}_{ij}=⟨u^L(z_{i_k}),u^L(z_{j_k}){⟩}^L,1\le k,i,j\le N,$$

$$u^B(z)=\sum _{k=1}^N{c}_k(z)u^H(z_{i_k}).$$

$$\frac{f^{n+1}-f^n}{\Delta t}+v\cdot {\nabla }_x{f}^n=\frac{\mathcal{Q}(f^n,f^n)-{\beta }^n(M^n-f^n)}{\epsilon }+\frac{{\beta }^n(M^{n+1}-f^{n+1})}{\epsilon }.$$

$$f^{n+1}=\frac{1}{1+\frac{\Delta t}{\epsilon }{\beta }^n}\left[f^n-\Delta tv\cdot {\nabla }_x{f}^n+\frac{\Delta t}{\epsilon }\left(\mathcal{Q}(f^n,f^n)-{\beta }^n(M^n-f^n)+{\beta }^n{M}^{n+1}\right)\right],$$

$${\beta }^n=\sup \left|\frac{\mathcal{Q}(f^n,f^n)}{f^n-M^n}\right|,$$

$$\frac{W^{n+1}-W^n}{\Delta t}+{\nabla }_x\cdot ⟨vm{f}^n⟩=0.$$

$$\frac{W^{n+1}-W^n}{\Delta t}+{\nabla }_x\cdot ⟨vm{M}^n⟩=0,$$

$$f=\mathcal{M}+\epsilon \sqrt{\mathcal{M}}h,$$

Full text

Keywords. Boltzmann equation, uncertainty, bi-fidelity models, multiple scales, stochastic collocation

1 Introduction

Kinetic equations are widely used in classical fields such as rarefied gas, plasma physics, astrophysics, also in emerging areas such as semiconductor device modeling, social and biological sciences. They model the non-equilibrium dynamics of a large number of particles from a statistical point of view [13]. The Boltzmann equation, as one of the most fundamental kinetic equations [11], is an integro-differential equation describing the time evolution of probability density distribution of particles in rarefied gas.

There have been extensive studies on the Boltzmann and related kinetic models, both in theory and numerical computation [67, 21]. One of the main computational challenges is that the kinetic problems often encounter multiple temporal and spatial scales, characterized by the Knudsen number, the dimensionless mean free path, that may vary in orders of magnitude in the computational domain, covering the regimes from fluid, transition to rarefied. Asymptotic-preserving (AP) schemes, which preserve the asymptotic transition from one scale to another at the discrete level, have been shown to be an effective computational paradigm in the past two decades [38, 39]. They allow efficient numerical approximations in all regimes–coarse mesh and large time steps can be used even in the fluid dynamic regime, without numerically resolving the small Knudsen number. For the space inhomogeneous Boltzmann equation, AP schemes were first designed using BGK-operator penalty based method [26]. Other approaches include the exponential integrator based methods [20] or micro-macro decomposition [44, 28]. Despite the successful development of AP solvers in recent years, the complexity and memory cost for computing the collision term still remains challenging [58, 52, 29, 27].

Another challenge, which has been ignored in the community until recently, is the issue of uncertainties in kinetic models. Kinetic equations, derived from $N$-body Newton’s equations via the mean-field limit [8], typically contain an integral operator modeling interactions between particles. Calculating the collision kernel from first principles is extremely complicated and not possible for complex particle systems, thus only empirical formulas are used for general particles [12], which means that the collision kernel contains many uncertainties. Other sources of uncertainties may also come from inaccurate measurements of the initial or boundary data, forcing or source terms. See for example [41, 23, 34, 50, 17, 48, 49] for recent efforts on uncertainty quantification for kinetic equations, in particular [64, 35], where numerical schemes for the Boltzmann equation with multi-dimensional random inputs have been studied.

One widely used method in uncertainty quantification is stochastic collocation (SC) method, especially in conjunction with the gPC expansion and high-performance grids such as sparse grids. There have been many works developed, for example [1, 36, 54, 56, 31, 68, 6]. One of the challenges central to collocation approaches is the simulation cost. For many complex systems, in particular, the multiscale Boltzmann equation with multi-dimensional random inputs we are studying, an accurate high-fidelity deterministic simulation can be so time-consuming and memory demanding that only a few high-fidelity simulations can be afforded. As many stochastic algorithms such as SC require repetitive implementations of the deterministic solver, the overall accurate stochastic simulation can be difficult and even computationally infeasible.

Fortunately, there usually exist some approximate, less complex low-fidelity models for practical problems. Compared to the high-fidelity models, these low-fidelity models usually contain simplified physics and/or are simulated on a coarser physical mesh, and consequently, own a cheaper computational cost. Although their accuracy may not be high, the low-fidelity models are designed in such a way that they can resolve or capture certain important features of the underlying problem and produce reliable and qualitative predictions. Despite numerous multi-fidelity algorithms have been developed in different communities from different perspectives [30, 14, 71, 42, 66, 61, 60, 59, 55, 25, 53, 72, 73], there have not been many attempts for kinetic equations with uncertainty in the multi-fidelity setting, except for the recent work [22]. In [22], the authors take the steady state or approximated time-dependent solution with the asymptotic behavior close to the fluid limit as control variate models to accelerate the convergence of standard Monte Carlo methods.

Despite numerous multi-fidelity algorithms have been developed for different applications [63, 62, 70, 69, 57, 45], the goal of our work is to adapt the bi-fidelity method developed in [53, 72, 73] to efficiently approximate the high-fidelity solutions of the Boltzmann equation with multi-dimensional random parameters and multiple scales. In rarefied gas dynamics, fluid models are derived when the mean free path of a particle is very small compared to the typical macroscopic length. One can perform a Hilbert or Chapman-Enskog expansion ([11]) of the solution to the Boltzmann equation in powers of the Knudsen number $\varepsilon$. At the leading order in $\varepsilon$, the distribution function approaches a local equilibrium – a Maxwellian whose parameters are the fluid variables (density, mean velocity and temperature) governed by the compressible Euler equations. When $\varepsilon$ is small, the Euler equations provide us a good accuracy in the physical space and numerical efficiency. Motivated by the above observation, we take advantage of this multiscale nature of the kinetic problem and choose the low-fidelity model based on the Euler equations in this work. More specifically, we connect an $\varepsilon$-dependent microscopic model (Boltzmann equation) and its macroscopic model (Euler equations) through the corresponding macroscopic quantities of the Boltzmann equation. Our numerical experiments demonstrate that the bi-fidelity solutions can approximate the high-fidelity solutions well at a much reduced computational cost.

This paper is outlined as follows: In Section 2, we give an introduction of the Boltzmann equation and its macroscopic equations. Section 3 reviews the framework of SC method with multifidelity models and discusses how it is adapted to our problem under study. Section 4 establishes the convergence of the bi-fidelity approximation to the high-fidelity solution under suitable assumptions. In section 5, we provide extensive numerical experiments to illustrate the effectiveness and efficiency of our proposed method, where kinetic, fluid and mixed regimes are carefully examined. We conclude the paper in Section 6.

2 Introduction of the Boltzmann equation

2.1 The Boltzmann equation with uncertainty

We first give an introduction to the classical (deterministic) Boltzmann equation, known as one of the most celebrated kinetic equations for rarefied gas. A dimensionless form reads

[TABLE]

where $f(t,x,v)$ is the probability density distribution function, modeling the probability of finding a particle at time $t>0$, at position $x$ in a bounded domain $\Omega\subset\mathbb{R}^{d_{x}}$, with velocity $v\in\mathbb{R}^{d_{v}}$, where $d_{x}$ and $d_{v}$ are the dimensions of the $x$ and $v$ variables. Periodic boundary condition is considered. The parameter $\varepsilon$ is the Knudsen number defined as the ratio of the mean free path over a typical length scale such as the size of the spatial domain. The collision operator $\mathcal{Q}$ is a quadratic integral operator modeling the binary elastic collision between particles, and is given by

[TABLE]

$(v,v_{\ast})$ and $(v^{\prime},v_{\ast}^{\prime})$ are the velocity pairs before and after the collision, during which the momentum and energy are conserved; thus $(v^{\prime},v_{\ast}^{\prime})$ can be expressed in terms of $(v,v_{\ast})$ as follows:

[TABLE]

with the vector $\sigma$ the scattering direction varying on the unit sphere $\mathbb{S}^{d_{v}-1}$. The collision kernel $B$ is a non-negative function of the form $B(v,v_{\ast},\sigma)=B(|v-v_{\ast}|,\cos\theta)$, where $\cos\theta=\frac{\sigma\cdot(v-v_{\ast})}{|v-v_{\ast}|}$ is the deviation angle. We consider the variable hard sphere (VHS) model [7], with a commonly used form for the collision kernel:

[TABLE]

where $b$ is a positive constant, $\gamma>0$ corresponds to the hard potential, and $\gamma<0$ is the soft potential. Denote the vector

[TABLE]

then

[TABLE]

which correspond to the conservation of mass, momentum and kinetic energy of the collision operator.

The celebrated Boltzmann’s H-theorem gives the dissipation of entropy ([11]):

[TABLE]

Furthermore, the equality holds if and only if $f$ reaches the equilibrium state

[TABLE]

which is known as the Maxwellian. Here $\rho$, $u$ and $T$ are the density, bulk velocity and temperature, respectively:

[TABLE]

There are many sources of uncertainties in the Boltzmann equation, such as the initial data, boundary data, and collision kernel. We introduce the Boltzmann equation with uncertainty

[TABLE]

Here $z\in I_{z}$ is a $d$-dimensional random parameter with probability distribution $\pi(z)$ known in priori characterizing the uncertainty in the system. Without loss of generality, we only consider periodic boundary condition in space throughout this paper.

2.2 The macroscopic fluid equations

When the Knudsen number $\varepsilon>0$ becomes very small, the macroscopic fluid dynamics describing the evolution of averaged quantities such as the density $\rho$, momentum $\rho u$ and temperature $T$ of the gas, namely, the compressible Euler or Navier-Stokes equations, become adequate [2, 8].

Multiplying (2.11) by $m(v)$ and integrating with respect to $v$, by using the conservation property of $\mathcal{Q}$ given in (2.6), one gets a non-closed system of conservation laws

[TABLE]

where $E$ is the total energy defined by

[TABLE]

with $\langle\,\cdot\,\rangle$ denoted as a velocity average of the argument,

[TABLE]

Here $\mathbb{P}=\langle(v-u)\otimes(v-u)f\rangle$ is the pressure tensor, and $\mathbb{Q}=\frac{1}{2}\langle(v-u)|v-u|^{2}f\rangle$ is the heat flux vector. Note that the variables $\rho$, $u$ and $E$ in (2.12) depend on the random parameter $z$.

When $\varepsilon\to 0$, $f\to M(v)_{\rho,u,T}$. We can approximate $f$ by $M(v)_{\rho,u,T}$ and use the expression (2.7), $\mathbb{P}$ and $\mathbb{Q}$ become

[TABLE]

where $p=\rho\,T$ is the pressure, $\mathbb{I}$ is the identity matrix. Then (2.12) reduces to the compressible Euler equations of gas dynamics for a mono-atomic gas:

[TABLE]

which is known as a first order approximation with respect to $\varepsilon$ to the Boltzmann equation (2.11). By the Chapman-Enskog expansion, the compressible Navier-Stokes equations give a second order approximation in $\varepsilon$ to the distribution function of the Boltzmann equation [11].

3 A stochastic collocation method with bi-fidelity models

In this section, we first briefly review an efficient bi-fidelity approximation to the high-fidelity solution studied in [53, 72], then we shall discuss the motivation of our choices of the low-fidelity model in our current study.

3.1 A bi-fidelity Algorithm

Assume we have access to the high-fidelity solutions $u^{H}(z)$ and low-fidelity solutions $u^{L}(z)$. Let $M$ be the number of affordable low-fidelity simulation runs, which can be very large. $N$ denotes the number of high-fidelity simulation runs that can be afforded and is often very small, i.e., $M\gg N$. Let $\gamma=\{z_{1},\cdots,z_{k}\}$, $k\geq 1$ be a set of sample points in $I_{z}$. Denote the low-fidelity snapshot matrix $u^{L}(\gamma)=[u^{L}(z_{1}),\cdots,u^{L}(z_{k})]$ and the corresponding low-fidelity approximation space

[TABLE]

Similarly, the high-fidelity snapshot matrix and the correponding high-fidelity approximation as follows:

[TABLE]

The bi-fidelity algorithm for approximating the high-fidelity solution consists of offline and online stages. In the offline stage, we employ the cheap low-fidelity model to explore the parameter space to find the most important parameter points. Within the online stage, we learn the approximation rule from the low-fidelity model for any given $z$, and apply it to construct the bi-fidelity approximation. The detailed algorithm is summarized in Algorithm.1.

Most of the steps in this algorithm are straightforward. It would be instructional to provide details for Step 3 (point selection) and Step 6 (bi-fidelity reconstruction).

Point selection. To select the subset $\gamma_{N}$, we shall search the parameter space by the greedy algorithm proposed in [53, 72]. Start with a trivial subspace $\gamma_{0}=\varnothing$, and assume that the first $k-1$ important points $\gamma_{k-1}=\{z_{i_{1}},\cdots,z_{i_{k-1}}\}\subset\Gamma$ have been selected. We shall choose the next point $z_{i_{k}}\in\Gamma$ as the point that maximizes the distance between its corresponding low-fidelity solution and the approximation space $U^{L}(\gamma_{k-1})$, spanned by the low-fidelity solutions on the existing point set $\gamma_{k-1}$, i.e.,

[TABLE]

where $\text{dist}(v,W)$ is the distance function between $v\in u^{L}(\Gamma)$ and subspace $W\subset U^{L}(\Gamma)$. The greedy procedure essentially serves the purpose of searching the linear independent basis set in the parameterized low-fidelity solution space. We remark that the whole algorithm allows an efficient implementation by standard linear algebra operations. See [53, 72] for more technical details.

Bi-fidelity approximation. In the offline stage, we have constructed the low- and high-fidelity approximation space, $U^{L}(\gamma_{N})$ (step 3) and $U^{H}(\gamma_{N})$ (step 4), respectively. During the online stage, for any given sample point $z\in I_{z}$, we shall project the corresponding low-fidelity solution $u^{L}(z)$ onto the low-fidelity approximation space $U^{L}(\gamma_{N})$:

[TABLE]

where $\mathcal{P}_{V}$ is the projection operator onto a Hilbert space $V$ and the corresponding projection coefficients $\{c_{k}\}$ are computed by the following projection:

[TABLE]

where ${\bf G}^{L}$ is the Gramian matrix of $u^{L}(\gamma_{N})$, defined by

[TABLE]

with $\langle\cdot,\cdot\rangle^{L}$ the inner product associated with the approximation space $U^{L}(\gamma_{N})$.

These low-fidelity coefficients $\{c_{k}\}$ serve as the surrogate of the corresponding high-fidelity coefficients of $u^{H}(z)$. Therefore, the sought bi-fidelity approximation of $u^{H}(z)$ can be constructed as follows:

[TABLE]

We emphasize that if the low-fidelity model can mimic the variations of the high-fidelity model in the parameter space, the low-fidelity coefficients can be a good approximation of the corresponding high-fidelity coefficients for a given sample $z$. We refer interested readers to [53, 72, 33] for details of the error analysis and justifications.

It is worth noting that since the number of low-fidelity basis is typically small ($\mathcal{O}(10)$ in our numerical tests), the cost of computing the low-fidelity projection coefficients by solving the linear system (3.2) is negligible. The dominant cost of the online step is one low-fidelity simulation run. If the low-fidelity solver is much cheaper than the high-fidelity solver, the speedup during the online stage can be significant.

3.2 The high- and low-fidelity models in our problem

Our purpose is to efficiently approximate high-fidelity solutions for the uncertain Boltzmann equation (2.11) for a fixed $z$, which is solved by a deterministic AP solver discussed in section 3.2.1. It is well known that existing solvers for deterministic kinetic equations are time-consuming and memory demanding due to its high-dimensional nature in the physical space. With the random parameter, it is more challenging to fully sweep the multi-dimensional parameter space by solving the Boltzmann equation repeatedly, especially given the complicated nonlinear collision operator in our model.

To mitigate this computational cost, we consider to choose the compressible Euler equations (2.14) as our low-fidelity model. It is a first-order approximation to the Boltzmann equation, which can mimic the variations of macroscopic quantities of the Boltzmann equation in the fluid regime up to a certain accuracy. Besides, it is worth noting that the macroscopic quantities do not depend on the velocity $v$ in (2.12). Therefore, solving the deterministic Euler equation is much easier and more efficient in terms of memory and computational time compared to solving the deterministic Boltzmann equation (2.1). These facts motivate us to choose the Euler equation as the low-fidelity model in our numerical experiments. A comparison of the computation cost (CPU time) for the two models are given in Section 5.

3.2.1 A high-fidelity solver

To solve the high-fidelity model Boltzmann equation, we shall resort to a high-fidelity asymptotic-preserving (AP) solver. There have been many works in developing robust numerical schemes for kinetic equations in the framework of the asymptotic-preserving scheme, see for example [5, 44, 43, 40, 16, 28]. As pointed in [26], AP scheme for the kinetic equation has two major merits:

1. as the Knudsen number $\varepsilon$ go to zero, it automatically becomes a consistent and stable scheme for the limiting fluid equation, with the stability condition independent of $\varepsilon$ (i.e., $\Delta t$ independent of $\varepsilon$); 2) the implicit collision terms can be implemented explicitly, free of Newton-type nonlinear algebraic solvers. Compared with multi-physics domain decomposition methods [18], AP schemes avoid the coupling of physical equations of different scales where coupling conditions and interface locations are difficult to determine. In contrast with many existing multiscale solvers, the AP schemes only require solving one equation – the kinetic equation and it becomes a robust macroscopic solver automatically when $\varepsilon\to 0$.

For our problem, we shall employ an AP scheme developed in [26] for the deterministic rescaled Boltzmann equation (2.1) as our high-fidelity solver. The main idea of [26] is to penalize the collision term $\mathcal{Q}(f,f)$ by the BGK operator $P(f)=M-f$ which can be inverted easily, thus the scheme can be solved explicitly. Let the initial distribution function be $f_{\text{in}}$ and consider periodic boundary conditions. The basic scheme consists of the following two major steps:

We first discretize (2.1) in time by the following first-order semi-discrete scheme:

[TABLE]

$f^{n+1}$ can be rewritten as follows:

[TABLE]

where $\beta$ is some constant that depends on the spectral radius of the linearized collision operator of $\mathcal{Q}$ around the Maxwellian $M$. We refer to [26, Section 2] on the intuition and justification of choosing $\beta$. For example, one can set

[TABLE]

and other choices $\beta$ are also available [26, 28]. We numerically evaluate the collision term $\mathcal{Q}(f^{n},f^{n})$ in (3.5) by applying the fast spectral method developed in [52].

• 2.

Though the above equation appears to be implicit due to $M^{n+1}$, it can be solved explicitly, thanks to the conservation property of $\mathcal{Q}(f,f)$ and the BGK operator $P(f)$. By multiplying the equation (3.5) with the vector $m(v)$ in (2.5), we can get the following equation:

[TABLE]

where $W:=(\rho,\rho u,E)$ that consists of the macroscopic quantities (mass, momentum and energy). With $W^{n+1}$, one computes $M^{n+1}$ from the Maxwellian (2.7). Finally, we can update $f^{n+1}$ explicitly from (3.6).

For the spatial discretization in (3.6), we employ a second order upwind MUSCL scheme as in [26], and a second order minmod slope limiter is used to suppress possible spurious oscillations near discontinuities or sharp gradients [46]. In addition to (3.7), a second order TVD scheme with a minmod slope limiter is also applied, see [5, 28] for details of implementation.

3.3 A low-fidelity solver

For the low-fidelity model, instead of solving the Euler system (2.14) directly, we shall semi-discretize its equivalent form (3.7) with $f$ replaced by the Maxwellian $M(v)_{\rho,u,T}$, i.e.,

[TABLE]

where the relation between $W:=(\rho,\rho u,E)$ and $M$ is given in (2.13) and (2.7). The initial data of $\rho$, $u$ and $E$ are obtained from the initial distribution $f_{\text{in}}$ for the Boltzmann equation, by using (2.8) and (2.13). That is, the initial data for the low- and high-fidelity models are consistent. The scheme (3.8) is numerically solved in the same way as equation (3.7) in the AP solver for the Boltzmann equation. Since Euler system is marching the macroscopic quantities, instead of marching the distribution solution $f$ to the Boltzmann equation, the scheme (3.8) can be solved with a much reduced computational cost and memory consumption.

Remark 3.1**.**

We remark that instead of taking the solution $f$ to the Boltzmann equation via the scheme (3.6) as the high-fidelity solutions, we consider its corresponding macroscopic quantities of interest $\mathcal{U}^{H}(z)=[\rho^{H}(z),u^{H}(z),T^{H}(z)]^{\top}$ as the high-fidelity snapshot solutions in order to connect the macroscopic quantities computed from the low-fidelity models. The low-fidelity solutions $\mathcal{U}^{L}(z)=[\rho^{L}(z),u^{L}(z),T^{L}(z)]^{\top}$ we considered are computed from the Euler system by using (3.8). During the point selection step to construct $\gamma_{N}$ in Algorithm 1, we shall select the important parameter points based on the concatenated macroscopic quantity snapshot, namely $\mathcal{U}^{L}(z)=[\rho^{L}(z),u^{L}(z),T^{L}(z)]^{\top}$.

Remark 3.2**.**

We acknowledge that there could be other choices of low-fidelity models that lead to more accurate bi-fidelity approximation, e.g., the compressible Navier-Stokes equations. To estimate if the low-fidelity model would be useful for constructing a reasonable accurate bi-fidelity approximation, one can explore an a priori estimate developed in the recent work [33].

4 Accuracy and convergence analysis

To establish the accuracy and convergence results, we first give a summary of the hypocoercivity framework and notations used in [50], then introduce the relation between the solutions to the Boltzmann and compressible Euler system in suitable norms. To study the difference between the high- and bi-fidelity solutions, one can split it into two parts: the projection error and the remainder. In section 4.1, we show the estimate for the projection error in Theorem 4.1. In section 4.2, we give the regularity of high-fidelity solution, then prove the accuracy and convergence results of our bi-fidelity method adapted to the Boltzmann equation in Theorem 4.2.

4.1 The projection error

The subject of hydrodynamic limits and rigorous derivations of macroscopic models such as the fundamental PDEs of fluid mechanics from the kinetic theory of gases is a challenging task and has been studied for decades, see for example [24, 51, 10, 4, 3, 37, 32]. We shall show that for each fixed $z\in I_{z}$, the error between solutions to the Euler system and the macroscopic quantities (2.8) obtained from the Boltzmann equation (with consistent initial data) is small and of order $\varepsilon$, which will be described in (4.3).

Hypocoercivity framework. First, we review the hypocoercivity framework and notations for the norms used in [50]. Let $f$ be the solution to the Boltzmann equation (2.1). Consider a linearization around the global equilibrium and perturbation of $f$:

[TABLE]

with

[TABLE]

then $h$ satisfies the perturbed equation

[TABLE]

where the linearized operator $\mathcal{L}$ and the nonlinear operator $\mathcal{F}$ are given by

[TABLE]

Denote $\partial_{l}^{j}:=\partial/\partial v_{j}\,\partial/\partial x_{l}$ for multi-indices $j$ and $l$. Introduce the following Sobolev norms:

[TABLE]

Refer to [9, Theorem 2.5], we extend its analysis to our case of the Boltzmann equation in the acoustic regime. Let $h_{\varepsilon}$ be the perturbed solution to the linearized equation (4.1). Suppose the initial data for (4.1) and (2.14) are consistent for each $z$. If the initial distribution $h_{in}\in\text{Null}(\mathcal{L})$ and $h_{in}\in H_{x,v}^{s}$, then for each $z$, $\left(h_{\varepsilon}\right)_{\varepsilon>0}$ converges strongly to

[TABLE]

in $L_{[0,T]}^{2}H_{x}^{s}L_{v}^{2}$ as the Knudsen number $\varepsilon\to 0$, where $\rho$, $u$, $T$ (with $E$ obtained by (2.13)) satisfy the Euler system (2.14). We adapt our acoustic scaling to [9, Theorem 2.5] and get the follows: For all $z\in I_{z}$, if $h_{in}$ belongs to $H_{x}^{s}L_{v}^{2}$, then

[TABLE]

where $\forall T>0$, $V_{T}(\varepsilon)$ is defined as

[TABLE]

Error splitting. Let $\left\langle\,\cdot\,\right\rangle^{H}$ be an inner product space corresponding to the high-fidelity solution and $\left|\left|\,\cdot\,\right|\right|^{H}$ be the corresponding induced norm, see [53]. For each $z$, to study the total error $\left|\left|u^{H}(z)-u^{B}(z)\right|\right|^{H}$, one can split it into two parts:

[TABLE]

[53, Lemma 4.3] shows the estimate for the second term:

[TABLE]

where we refer $c$ to $\varepsilon_{1}+\varepsilon_{2}+\varepsilon_{1}\varepsilon_{2}$ in [53, Lemma 4.3], which is small, based on the reasonable assumptions made there. The last term above is related to the non-invertibility of high-fidelity Gramian matrix and usually negligible. Here ${\bf G}^{L}$ (or ${\bf G}^{H}$) is the Gramian matrix of $u^{L}(\gamma)$ (or $u^{H}(\gamma)$) given by (3.3), the vector ${\bf f}^{L}$ has entries

[TABLE]

and ${\bf Q}:={\bf I}-{\bf P}$ is the orthogonal projection onto its kernel (with ${\bf P}$ the orthogonal projection matrix onto its range), see [53]. In addition, since $\left|\left|P_{U^{H}(\gamma_{N}^{L})}u^{H}(z)\right|\right|^{H}\leq||u^{H}(z)||^{H}$, thus

[TABLE]

Projection error. The rest of this section will study the estimate for the projection error (the first term on the right-hand-side of (4.4)) and conclude it in Theorem 4.1. We now adapt the analysis in [53, subsection 4.1] and incorporate our high-fidelity (Boltzmann) and low-fidelity (Euler) models, by utilizing the knowledge of (4.3). Denote

[TABLE]

The “best” achievable distance for approximation from a general $N$-dimensional subspace is the Kolmogorov $N$-width, defined by

[TABLE]

The following is similar to [53, Lemma 4.1], except that we use the inequality (4.3). Since it is lengthy, we present as Lemma 1.1 with its proof in the Appendix.

Theorem 4.1**.**

If all the assumptions in Lemma 1.1 are satisfied for each $n=0,1,\cdots,N-1$, then

[TABLE]

where $d_{N}(u^{H}(I_{z}))$ is the $N$-width of the functional manifold $u^{H}(I_{z})$, and the constant $C=\sqrt{2}(\delta_{1}-\delta_{2}\,\varepsilon)^{-1}$ with $0<\delta_{1}-\delta_{2}\,\varepsilon<1$.

Proof.

Lemma 1.1 shows that it is a weak greedy procedure to use the nodal choices of $z_{n}^{L}$: $\exists\,0<\delta_{1}-\delta_{2}\,\varepsilon<1$ such that for all $n=0,1,\cdots,N-1$,

[TABLE]

and [19, Corollary 3.3] indicates that (4.9) holds with $C=\sqrt{2}(\delta_{1}-\delta_{2}\,\varepsilon)^{-1}$. ∎

4.2 Smoothness of the solution and convergence analysis

In this section, we first study the smoothness of the high-fidelity solution $u^{H}:I_{z}\mapsto V^{H}$, where $z$ is a multivariate random parameter and $V^{H}$ is a Hilbert space with an inner product $\left\langle\,\cdot,\cdot\right\rangle^{H}$. Then we shall establish an estimate bound for the Kolmogorov width in Theorem 4.1, finally combine all the arguments and show the convergence result for our bi-fidelity procedure in Theorem 4.2.

We introduce the standard multivariate notation. Denote the countable set of “finitely supported” sequences of nonnegative integers by

[TABLE]

with $|\nu|:=\sum_{j\geq 1}|\nu_{j}|$. For $\nu\in\mathcal{F}$ supported in $\{1,\cdots,J\}$, one defines the partial derivative in $z$

[TABLE]

and the multi-factorial $\nu!:=\prod_{j\geq 1}\nu_{j}!$, where $0!:=0$.

Regularity of $u^{H}$. It is not hard to extend the sensitivity analysis in [50] to the multi-dimensional random variable case, see Remark 2.8 in [65] on the extension. [50] tells us that 1) the uncertainties from the initial data and collision kernel (under suitable assumptions) will eventually diminish and the solution will exponentially decay in time to the deterministic global equilibrium $\mathcal{M}$; 2) the regularity of the initial data in the random space is preserved at later time.

Let $h_{\varepsilon}$ be a perturbed solution to the equation (4.1). If its random initial data satisfies for all $r$,

[TABLE]

then at time $t>0$,

[TABLE]

Moreover, $h$ is analytic in the random space, meaning that

[TABLE]

See Appendix B for a discussion on the linearized Boltzmann case.

We now use a weaker version by letting $s=1$ in (4.11) and (4.12). If the initial data $||h_{\varepsilon}^{in}(z)||_{H_{x,v}^{1,r}L_{z}^{\infty}}\leq C_{I}$, then

[TABLE]

By the definition of $u$ containing perturbed macroscopic quantities (see [9, section 2.2.4]) by the Cauchy-Schwarz inequality, one easily gets

[TABLE]

for $|\nu|\leq r$, where $C$ and $\tau>0$ are all generic constants independent of $\varepsilon$. We assume that the high-fidelity solution $u^{H}$ follows a similar behavior as the analytic solution $u$ (computed from $h_{\varepsilon}$) to the Boltzmann equation, with an error that depends on the numerical scheme used in the high-fidelity model. If the initial distribution of the high-fidelity model satisfies

[TABLE]

then for a fixed time $t>0$ and $|\nu|\leq r$,

[TABLE]

where $\xi$ depends on the order and discretization parameters $\Delta t$, $\Delta x$, $\Delta v$ used in the high-fidelity solver. Thus for all $z\in I_{z}$, one gets

[TABLE]

Note that $C_{I}$, $C$, $C^{\prime}$ and $\tau$ in the inequalities (4.11)–(4.15) are all positive generic constants independent of $\varepsilon$.

We make the following assumption on the random collision kernel:

Assumption 1**.**

Assume the collision kernel take the form

[TABLE]

and $\left(\psi_{j}\right)_{j\geq 1}$ be an affine representer (see definition in [15]) of the cross section $b$, that is,

[TABLE]

where the sequence $\left(||\psi_{j}||_{L^{\infty}(\eta)}\right)_{j\geq 1}\in\ell^{p}$ for $0<p<1$ (see [15]). One also assumes that

[TABLE]

for all $\eta\in[-1,1]$, $|\nu|\leq r$. Here $C$, $C_{0}$, $C_{1}$, $C_{2}$ are all positive constants.

The above (4.16)–(4.18) extend the conditions in [50, Theorem 4.4 (ii)] from one-dimensional random space to the multi-dimensional random space.

An upper bound for $N$-width. We can now utilize the result in [15, Corollary 3.11], which works for general parametric PDEs. Define the norm

[TABLE]

where $V$ is the physical space considered, in our case $V=L^{2}_{x}$. Under Assumption 1, by the analyticity of $u^{H}$ given in (4.15), for $z\in I_{z}$, one obtains

[TABLE]

where $\left(P_{k}\right)_{k\geq 0}$ is the sequence of renormalized Legendre polynomials on $[-1,1]$, $\Lambda_{N}$ is the set of indices that corresponds to the $N$ largest $||w_{\nu}||_{V}$, and the constant $C:=\left|\left|\left(||w_{\nu}||_{V}\right)_{\nu\in\mathcal{F}}\right|\right|_{\ell^{p}}<\infty$. By (4.15), one formally gets $C=c\,e^{-\varepsilon\tau t}+\xi$.

Recall (4.8) and the definition $d^{H}(v,V_{N}):=\min_{w\in V_{N}}||v-w||_{V}$ (refer to [15, equation (8.3)]), the best achievable error in $L^{\infty}(I_{z},V)$ is described by the $N$-width of the solution manifold $\mathcal{M}=u^{H}(I_{z})$:

[TABLE]

One can use the polynomial approximation bound in $L^{\infty}(I_{z},V)$ to estimate an upper bound for $d_{N}(\mathcal{M})_{V}$, as studied in [15, Section 4]. Let the $N$-dimensional subspace

[TABLE]

For $N\geq 1$, one observes that

[TABLE]

where $\sum_{\nu\in\Lambda_{N}}w_{\nu}P_{\nu}$ is the truncated Legendre expansion and (4.19) is used in the last inequality.

We now conclude with our main result on convergence analysis:

Theorem 4.2**.**

If the assumptions for the random initial data, random collision kernel, namely (4.14) and Assumption 1 are satisfied, for fixed time $t>0$ and fixed numerical discretization parameters $\Delta t$, $\Delta x$ and $\Delta v$, then for all $z\in I_{z}$,

[TABLE]

where $N$ is the size of the subspace $\gamma_{N}$ in Algorithm 1, and $q$ is given in (4.19) with $p$ depending on the $\ell^{p}$-summability assumption of $(\psi_{j})_{j\geq 1}$, $C_{1}=\mathcal{O}\left(\frac{1}{\delta_{1}-\delta_{2}\varepsilon}\right)$, $C_{2}$ and $\tau$ are constants that depend on the initial data $u^{in}$ and Assumption 1 on the collision kernel. $\delta_{1}$, $\delta_{2}$ are all sufficiently small with $0<\delta_{1}-\delta_{2}\,\varepsilon<1$. Definitions of ${\bf G}^{L}$, ${\bf G}^{H}$, ${\bf Q}$ and $f^{L}$ are given below (4.5). $\chi$ is associated to the order and discretization mesh in the high-fidelity solver. $\varepsilon$ is the Knudsen number in the Boltzmann equation (2.1).

Proof.

According to the inequalities (4.4), (4.6), (4.21), (4.15) and Theorem 4.1, one gets for all $z\in I_{z}$,

[TABLE]

here $C$, $C_{1}$, $C_{2}$, $\tau$ are all generic positive constants independent of $\varepsilon$, $C_{2}=c\,C^{\prime}$ which is small, and $\chi=C_{1}\sqrt{\xi}+\xi$. ∎

Theorem 4.2 indicates that the error between the bi- and high-fidelity solutions decays algebraically with respect to the number of high-fidelity runs $N$. The convergence rate $q/2$ is independent of the dimension of the random space and the regularity of the initial data; it only relates to the $\ell^{p}$ summability of the affine representer $(\psi_{j})_{j\geq 1}$ in Assumption 1.

Remark 4.3**.**

We make the following remarks:

The estimate in Theorem 4.2 may not be sharp. Deriving a sharper estimate requires a better understanding on the role of the Knudsen number $\varepsilon$ in the accuracy analysis.

• 2.

It is not our goal of the current work to establish stability and error analysis for the deterministic AP method for the multiscale Boltzmann equation in **[26]**, which is difficult due to the penalization used in the scheme thus has not been studied to our best knowledge. Thus deriving (4.15) from (4.13) rigorously remains challenging.

We hope to report more results regarding the above two issues in future researches.

5 Numerical Tests

To examine the performance of the method, we shall compute numerical errors in the following way: we choose a fixed set of points $\{\hat{z}_{i}\}_{i=1}^{n}\subset I_{z}$ that is independent of the point sets $\Gamma$, and evaluate the following error between the bi-fidelity and high fidelity solutions at a final time $t$:

[TABLE]

where $||\cdot||_{L^{2}(D)}$ is the $L^{2}$ norm in the physical domain $D=\Omega\times\mathbb{R}^{2}$. The error can be considered as an approximation to the average $L^{2}$ error in the whole space of $D\times I_{z}$.

Since our goal is to numerically solve the multiscale Boltzmann equation with random inputs, we solve the Boltzmann equation (2.1) as the high-fidelity AP solver discussed in Section 3.2. We assume the random collision kernel in the form of

[TABLE]

and consider Maxwell molecules, i.e., $\lambda=0$ in (5.2). The low-fidelity model is chosen as the Euler equation solved by the forward Euler in time and second-order MUSCL scheme in space, by using the same spatial and temporal resolutions as the high-fidelity model, but with a different number of quadrature points $N_{v}^{l}$ in velocity discretization.

In all the examples, the spatial domain is chosen to be $[0,1]$ with $N_{x}$ grid points, and periodic boundary condition is assumed except for the shock tube tests. The velocity domain is chosen as $[-L_{v},L_{v}]^{2}$ with $L_{v}=8.4$ and $N_{v}$ grid points in each dimension. Without loss of generality, the $d$-dimensional random variable ${\bf z}$ is assumed to follow the uniform distribution on $[-1,1]^{d}$. The training set $\Gamma$ is chosen to be $M=1000$ random samples of ${\bf z}$. We examine the error of bi-fidelity approximation with respect to the number of high-fidelity runs by computing the norm defined in (5.1) (evaluated over an independent set of $n=1000$ Monte Carlo samples). All numerical experiments are conducted by MATLAB R2018b on macOS Mojave system with 2.4 GHz Intel Core i5 processor and 8GB DDR3 memory),

5.1 A double-peak initial data test

We first consider the following initial data to mimic the Karhunen-Loeve expansion of the random field:

[TABLE]

The uncertain collisional cross section is given by

[TABLE]

Here ${\bf z}^{\rho}=\left(z_{1}^{\rho},\cdots,z_{d_{1}}^{\rho}\right)$, ${\bf z}^{T}=\left(z_{1}^{T},\cdots,z_{d_{1}}^{T}\right)$, and ${\bf z}^{b}=z_{1}^{b}$ represent the random variables in the collision kernel, initial density and temperature. Let the initial distribution $f_{0}$ follow a double-peak non-equilibrium initial data [26]. Set $d_{1}=7$, thus this is a $d=15$ dimensional problem in the random space. We use the Boltzmann equation as the high-fidelity model and the Euler system as the low-fidelity model, set $\Delta x=0.01$, $\Delta t=8\times 10^{-4}$ (in both the high- and low-fidelity models), $N_{v}^{h}=16$, and the final time $t=0.1$.

In Figure 1, we consider the fluid regime with $\varepsilon=10^{-4}$. This figure shows the mean $L^{2}$ errors of $\rho$, $u_{1}$, $T$ between the high- and bi-fidelity solutions with different quadrature points in velocity space. Here $\bf{u}_{1}$ in the figures below stands for the first component of the two-dimensional bulk velocity ${\bf u}$. It is clear that the error decays fast with the number of high-fidelity runs. In addition, when $N_{v}^{l}$ increases, the error between the high- and bi-fidelity solutions decreases. This is expected because the Euler equation solved by more quadrature points in velocity space can capture more information about the high-fidelity model.

In Figure 2, fluid regime is considered and we vary $\varepsilon$ from $\varepsilon=10^{-2}$ to $\varepsilon=10^{-4}$. The Euler equation is chosen as the low-fidelity model, solved by the same forward Euler in time and second-order MUSCL scheme in space, and the same spatial and temporal meshes as the Boltzmann equation in the high-fidelity model, and with $N_{v}^{l}=8$ velocity quadrature points. One observes that the smaller $\varepsilon$ is, the lower level the errors saturate. This is expected, because when the Knudsen number $\varepsilon$ approaches to zero, the Euler equation as the low-fidelity model commits less modeling error and can capture more information of the high-fidelity model.

In Figure 3, we investigate the performance of the bi-fidelity approximation for the kinetic regime with $\varepsilon=1$. Fast convergence of the mean $L^{2}$ errors with respect to the number of high-fidelity runs is observed. Even though $\varepsilon$ is relatively large compared to the previous two tests, a satisfactory accuracy in characterizing the behaviors of the solution in the random space is achieved in both cases: $N_{v}^{l}=8$ and $N_{v}^{l}=4$; and the errors with $N_{v}^{l}=8$ is smaller than that of $N_{v}^{l}=4$. On the right column of Figure 3, we plot the high-, low- and the corresponding bi-fidelity solutions (with $r=20$, $N_{v}^{l}=8$) for a particular sample point $z$. One observes that the high- and bi-fidelity solutions match quite well, whereas the low-fidelity solutions appear inaccurate at some spatial points. This example seems to indicate that although in the kinetic regime, the fluid description breaks down in the physical space, the bi-fidelity solution can still capture important variations of the high-fidelity model (Boltzman equation) in the random space.

5.2 Sod shock tube test

We next consider a more challenging problem where the initial data is discontinuous. Assume the random collision kernel in the form

[TABLE]

and the random initial distribution

[TABLE]

where the initial data for $\rho^{0}$, $u^{0}$ and $T^{0}$ is given by

[TABLE]

Here ${\bf z}^{b}=\left(z_{1}^{b},\cdots,z_{d_{1}+1}^{b}\right)$ and ${\bf z}^{T}=\left(z_{1}^{T},\cdots,z_{d_{1}}^{T}\right)$ represent the random variables in the collision kernel and initial temperature. Set $d_{1}=7$, then the total dimension $d$ of the random space is $15$. We use the Boltzmann equation as the high-fidelity model, and solve it by $\Delta x=0.01$, $\Delta t=8\times 10^{-4}$, and $N_{v}^{h}=24$, until the final time $t=0.15$. We shall employ the Euler equation as the low-fidelity model, and solve it with the same spacial and temporal resolution with the high-fidelity model but with $N_{v}^{l}=12$. We consider the fluid regime with $\varepsilon=10^{-4}$ in this test.

From the left column of Figure 4, we see a fast convergence of $L^{2}$ errors between the high- and bi-fidelity solutions. With only 10 high-fidelity runs, the bi-fidelity approximation can reach an accuracy level $\mathcal{O}(10^{-3})$ for a 15-dimensional problem in random space, while the low-fidelity approximation is quite poor with an accuracy level $\mathcal{O}(10^{-1})$. To further illustrate the performance of our bi-fidelity method, we compared the high-, low- and the corresponding bi-fidelity solutions (with $r=10$) for a particular sample point $z$. One observes that the high- and bi-fidelity solutions match really well, whereas the low-fidelity solutions seem to be quite inaccurate at some points in the spatial domain. Even in this case, the bi-fidelity solutions can approximate the high-fidelity solutions very well.

Figure 5 shows clearly that the mean and standard deviation of the bi-fidelity approximation of $\rho$, $u_{1}$ and $T$ agree well with the high-fidelity solutions by using only $10$ high-fidelity runs. The result is a bit surprising yet reasonable, suggesting that even though the Euler model may be inaccurate in the physical space, it still can capture the behaviors and characteristics of the solution to the Boltzmann equation in the random space. Moreover, since the high-fidelity model (Boltzmann) with $N_{v}^{h}=24$ costs approximately $43$ times of the low-fidelity solver (Euler) with $N_{v}^{l}=12$ (the former takes $30.6$ seconds, the latter takes $0.7$ seconds for one single run; a significant speedup is quite noticeable in this case.

5.3 Mixed regime test

The next test we shall consider is more challenging than the previous two tests. Because various scales are involved, good accuracy of the AP scheme for the Boltzmann equation is required for all ranges of $\varepsilon$. We consider a mixed regime with the Knudsen number $\varepsilon$ varying in space show in Figure. 6 and given by

[TABLE]

The random initial data and collision kernel are given by (5.7) and (5.8). The total dimension of the random space is $d=15$.

All the numerical parameters used in temporal and spatial discretizations are the same as that in Section 5.1. We solve the Boltzmann equation (2.1) for the high-fidelity solution with $N_{v}^{h}=24$, and the Euler system for the low-fidelity solution with $N_{v}^{l}=8$.

From the left column of Figure 7, we observe a fast convergence of $L^{2}$ errors between the high- and bi-fidelity solutions, where they saturate quickly when $r$ reaches about $25$. It is worth noting that the dotted lines that represent the errors between the high- and low-fidelity solutions are much larger $\mathcal{O}(10^{-1})$ compared to that between the high- and bi-fidelity solutions. This indicates that even though the low-fidelity solutions alone are relatively not accurate in the spatial domain, it might be still able to behave similarly in the random space, therefore the resulted bi-fidelity approximation based on a small number of high-fidelity runs (say $r=25$) can reach a reasonable accuracy level up to $\mathcal{O}(10^{-3})$.

The right column of Figure 7 shows the high-, low- and bi-fidelity solutions at a randomly chosen sample point $z$. One can see that the high- and bi-fidelity solutions match really well, whereas the low-fidelity solutions are not accurate. In addition, with $N=1000$ low-fidelity runs of the Euler model, together with only $25$ runs of the AP solver to the Boltzmann model, one can get the bi-fidelity solutions which are able to capture behavior of the solutions to the Boltzmann equation in the random space, up to an accuracy of $10^{-3}$; on the other hand, using the low-fidelity model (Euler equation) alone can not achieve this result, especially under the multiple scalings where $\varepsilon$ ranges from $10^{-3}$ to 1 (since the errors between macroscopic quantities calculated from the Boltzmann and Euler equation deteriorate when $\varepsilon$ becomes large). This observation certainly highlights the merits of our bi-fidelity method. Figure 8 presents the mean and standard deviation of $\rho$, $u_{1}$, $T$ by using $25$ high-fidelity runs. One can see that the high- and bi-fidelity solutions match well, indicating that the bi-fidelity solutions have captured well the characteristics of the macroscopic quantities in the random space.

Once we construct the bi-fidelity model, the online computational cost can be significantly reduced. In this example, the high-fidelity model (Boltzmann) with $N_{v}^{h}=16$ takes about $50$ times computation time of the low-fidelity model (the former takes $12.3$ seconds, while the latter takes $0.23$ seconds for a single run). Since the dominant cost of the bi-fidelity reconstruction for a high-fidelity solution lies in the corresponding low-fidelity run, a significant amount of computational cost is saved in our method.

6 Conclusion

In this work, we study the multiscale Boltzmann equation with multi-dimensional random parameters by a bi-fidelity collocation method [53, 72]. By choosing the low-fidelity solution as the solution of the corresponding first order macro-model – the compressible Euler equations with a consistent initial data, our bi-fidelity approximation can capture well the variations of macroscopic quantities computed from the high-fidelity AP solver of the Boltzmann equation with multiple scales, at a much reduced computational cost and memory footprint. An error analysis developed in [53] has been extended by incorporating the knowledge of the regularity of our high-fidelity and low-fidelity solutions. The computational accuracy and efficiency are demonstrated in various numerical examples and holds promise to accelerate the computation for more complex problems in multiscale kinetic equations with uncertainty.

Appendix A The Lemma proof

Lemma 1.1**.**

Let $z_{\ast}^{H}$, $z_{\ast}^{L}$ be the maximizers of (4.7a) and (4.7b), respectively. One assumes that

1. $\exists$ positive constants $K^{H}$, $K^{L}$, $C_{1}$, $C_{2}$ such that

[TABLE]

[TABLE]

2. $\exists\,\eta^{L}\leq 1$ and $\eta^{H}\leq 1$ such that

[TABLE]

[TABLE]

3. $\exists$ constants $A^{H}$ and $A^{L}$ satisfying

[TABLE]

such that

[TABLE]

[TABLE]

Then

[TABLE]

with $0<\delta_{1}-\delta_{2}\,\varepsilon<1$.

Proof.

Using the above assumptions, one gets

[TABLE]

Similarly, we have

[TABLE]

Therefore,

[TABLE]

where we assume that

[TABLE]

and $\exists\,\delta_{2}$ such that $0<\delta_{1}-\delta_{2}\,\varepsilon<1$ with

[TABLE]

∎

Appendix B Proof of analyticity for the linearized equation

For simplicity, consider the linearized Boltzmann equation under the acoustic scaling. The perturbative solution satisfies

[TABLE]

For simplicity, we write out the proof for one-dimensional random variable $z$. Assume the collision cross-section has the form $b(\cos\theta)=b_{0}(\cos\theta)+b_{1}(\cos\theta)\,z$, and we require $|b_{1}|\leq C_{b}$. Define the operator $\mathcal{T}_{\varepsilon}:=\frac{1}{\varepsilon}\mathcal{L}-v\cdot\nabla_{x}$, and the notation $h_{l}:=\partial_{z}^{l}h$ for all $l\geq 0$. Take $\partial_{z}^{l}$ of the equation (B.1), then

[TABLE]

where $\mathcal{L}_{b_{1}}$ is defined by substituting $b(\cos\theta)$ by $b_{1}(\cos\theta)$ in the linearized collision operator $\mathcal{L}$. Thus

[TABLE]

By the hypocoercivity theory and $\mathcal{L}$ being a bounded operator,

[TABLE]

where $\|h\|_{H_{\Lambda}^{k}}^{2}=\sum_{|j|+|l|=k}\left\|\frac{\partial^{2}}{\partial v_{j}\partial x_{l}}h\right\|_{\Lambda}^{2}$ and $\|h\|_{\Lambda}=\left\|h(1+|v|)^{\gamma/2}\right\|_{L^{2}_{x,v}}$, with $\gamma$ shown in (2.4). Since $||h_{l}||_{H_{\Lambda}^{k}}$ controls $||h_{l}||_{H^{k}}$, then by setting $g_{l}=||h_{l}||_{H^{k}}$,

[TABLE]

Using [47, Lemma 6] in a similar way, one gets

[TABLE]

Following [47, Theorem 5], if the initial data satisfies $||\partial_{z}^{l}h_{0}(z)||_{H^{k}}\leq R^{l}$, $\forall l\geq 0$, then at time $t$,

[TABLE]

where $C$, $C_{0}$, $C_{1}$ are independent of $l$ and $\varepsilon$. The convergence radius for $h$ at any point $z_{0}$ is defined by

[TABLE]

which is independent of $z_{0}$, thus $h$ (or $f$) is analytic in $z$.