Publisher main menu

Uniform spectral convergence of the stochastic Galerkin method for the linear transport equations with random inputs in diffusive regime and a micro–macro decomposition-based asymptotic-preserving method

Abstract

In this paper we study the stochastic Galerkin approximation for the linear transport equation with random inputs and diffusive scaling. We first establish uniform (in the Knudsen number) stability results in the random space for the transport equation with uncertain scattering coefficients and then prove the uniform spectral convergence (and consequently the sharp stochastic asymptotic-preserving property) of the stochastic Galerkin method. A micro–macro decomposition-based fully discrete scheme is adopted for the problem and proved to have a uniform stability. Numerical experiments are conducted to demonstrate the stability and asymptotic properties of the method.

This equation arises in neutron transport, radiative transfer, etc., that describes particles (for example neutrons) transport in a background media (for example nuclei), in which f(t, x, v, z) is the density distribution of particles at time \(t\ge 0\), position \(x \in (0,1)\). \(v=\Omega \cdot e_x = \cos \theta \in [-1, 1]\) where \(\theta \) is the angle between the moving direction and x-axis. \(\sigma (x,z)\), \(\sigma ^a(x,z)\) are total and absorption cross section, respectively. S(x, z) is the source term. \(\varepsilon \) is the dimensionless Knudsen number, the ratio between particle mean free path and the characteristic length (such as the length of the domain). The Dirichlet boundary conditions are given in the incoming direction by

We are interested in the problem that contains uncertainty in the collision cross section, source, initial or boundary data. The uncertainty is characterized by the random variable \(z\in {\mathbb {R}}^d\) with probability density function \(\omega (z)\). Thus in our problem \(f, \sigma , \sigma ^a\) and S all depend on z.

In recent years, there have been extensive activities to study partial differential equations or engineering problems with uncertainties. Many numerical methods have been introduced. In this article, we are interested in the polynomial chaos (originally introduced in Wiener’s work [26])-based stochastic Galerkin method which has been shown to be competitive in many applications, see [4, 27, 28]. The stochastic Galerkin method has been used for linear transport equation with uncertain coefficients [25]. Here we are interested in the problem that contains both uncertainty and multiscale. The latter is characterized by the Knudsen number \(\varepsilon \), which, in the so-called optically thin region (\(\varepsilon \ll 1\)), due to high scattering rate of particles, leads the linear transport equation to a diffusion equation, known as the diffusion limit [1, 3, 18]. For the past decades, developing asymptotic-preserving (AP) schemes for (deterministic) linear transport equation with diffusive scaling has seen many activities, see for examples [5, 7–9, 17, 19, 20, 22]. Only recently AP scheme for linear transport equation with both uncertainty and diffusive scaling was introduced in [15] (in the framework of stochastic Galerkin method, coined as s-AP method). See more related recent works along this line in [6, 12, 13]. A scheme is s-AP if the stochastic Galerkin method for the linear transport equation, as \(\varepsilon \rightarrow 0\), becomes a stochastic Galerkin method for the limiting diffusion equation. It was realized in [15] that the deterministic AP framework can be easily adopted to study linear transport equations with uncertain coefficients. Moreover, as shown in [6, 12], kinetic equations, linear or nonlinear, could preserve the regularity in random space of the initial data at later time, which naturally leads to spectral accuracy of the stochastic Galerkin method.

When \(\varepsilon \ll 1\), however, the energy estimates and consequently the convergence rates given in [6, 12] depend on the reciprocal of \(\varepsilon \), which implies that one needs the degree of the polynomials used in the stochastic Galerkin method to grow as \(\varepsilon \) decreases. In fact, this is typical of a standard numerical method for problems that contain small or multiple scales. While AP schemes can be used with numerical parameters independent of \(\varepsilon \), to prove this rigorously is not so easy and has been done only in a few occasions [5, 14]. A standard approach to prove a uniform convergence is to use the diffusion limit, as was done first in [5] in the deterministic case and then in [11] for the uncertain transport equation. See also the review article [7]. However such approaches might not give the sharp convergence rate.

In this paper, we provide a sharp error estimate for the stochastic Galerkin method for problem (1). This requires a sharp (\(\varepsilon \)-independent) energy estimate on high-order derivatives in the random space for f, as well as \([f]-f\) where [f] is the velocity average of f defined in (5) which is shown to be bounded even if \(\varepsilon \rightarrow 0\). Then the uniform in \(\varepsilon \) spectral convergence naturally follows, without using the diffusion limit.

The s-AP scheme in [15] uses the AP framework of [8] that relies on the even- and odd-parity formulation of the transport equation. In this paper, we use the micro–macro decomposition-based approach (see [22]) to develop a fully discrete s-AP method. The advantage of this approach is that it allows us to prove a uniform (in\(\varepsilon \)) stability condition, as was done in the deterministic counterpart in [23]. In fact, we will show that one can easily adopt the proof of [23] for the s-AP scheme.

The paper is organized as follows. In Sect. 2 we summarize the diffusion limit of the linear transport equation. The generalized polynomial chaos-based stochastic Galerkin method for the problem is introduced in Sect. 3 and shown formally to be s-AP. The uniform in \(\varepsilon \) regularity of the stochastic Galerkin scheme is proven in Sect. 4, which leads to a uniform spectral convergence proof. The micro–macro decomposition-based fully discrete scheme is given in Sect. 5, and its uniform stability is established in Sect. 6. Numerical experiments are carried out in Sect. 7. The paper is concluded in Sect. 8.

The micro–macro decomposition, a useful tool for the study of the Boltzmann equation and its fluid dynamics limit [24], and for the design of asymptotic-preserving numerical schemes for kinetic equations [2, 16, 22], takes the form

3 The gPC-stochastic Galerkin approximation

We assume the complete orthogonal polynomial basis in the Hilbert space \(H({\mathbb {R}}^d; \omega (z)\,\mathrm {d}z)\) corresponding to the weight \(\omega (z)\) is \(\{\phi _i(z), i=0,1, \ldots ,\}\), where \(\phi _i(z)\) is a polynomial of degree i and satisfies

4 The regularity in the random space and a uniform spectral convergence analysis of gPC-SG method

In this section, we assume \(\sigma ^a=S=0\) for clarity and periodic boundary condition

$$\begin{aligned} f(t, 0, v, z) = f(t, 1, v, z) \end{aligned}$$

(19)

for simplicity. We prove that, under some suitable assumptions on \(\sigma (z)\), the solution to the linear transport equation with random inputs preserves the regularity in the random space of the initial data uniformly in\(\varepsilon \). Then based on the regularity result, we conduct the spectral convergence analysis and error estimates for the gPC-SG method and will also obtain error bounds uniformly in\(\varepsilon \).

4.1 Notations

We first recall the Hilbert space of the random variable introduced in Sect. 3,

where \(Q=[0,1]\times [-1,1]\) denotes the domain in the phase space. For simplicity, we will suppress the dependence of t and just use \(\Vert f\Vert _\Gamma , \Vert f\Vert _{\Gamma ^k}\) in the following proof.

4.2 Regularity in the random space

We will study the regularity of f with respect to the random variable z. To this aim, we first prove the following lemma. For simplicity, we state and prove the following lemma only for one-dimensional case. Proof for high-dimensional case is identical except the change of coefficient.

Lemma 4.1

Assume \(\sigma (z)\ge \sigma _{\mathrm {min}}>0\), then for any integer k and \(\sigma \in W^{k,\infty }, g\in H^k\) we have

where C is clearly independent of \(\varepsilon \). This completes the proof of the theorem. \(\square \)

Theorem 4.1 shows the derivatives of the solution with respect to z can be bounded by the derivatives of initial data. In particular, the \(\Vert D^k f\Vert _{\Gamma }\) bound is independent of \(\varepsilon \)! This is crucial for our later proof that our scheme is s-AP. However, this estimate alone is not sufficient to guarantee that the whole gPC-SG method has a spectral convergence uniform in \(\varepsilon \) (since there is \(O(1/\varepsilon ^2)\) coefficient in front of the projection error, such we need \(O(\varepsilon ^2)\) estimation of \([ f] - f\) to cancel this coefficient). To this aim, we first provide the following lemma.

where \(C_{p+1}\) is a constant independent of \(\varepsilon \). So by mathematical induction, we complete the proof of the theorem. \(\square \)

Remark 4.1

We remark that all the above lemma and theorems are proved for \(z\in {{\mathbb {R}}}\) and \(\sigma \) depending only on z. However, our conclusions and techniques are not limited to these cases. For \(z\in {{\mathbb {R}}}^d\), it is straightforward to prove and for \(\sigma (x,z)\) also a function of x, we only need to modify the proof of Lemma 4.2 by using the same technique as in the proof of Theorem 4.1.

4.3 A spectral convergence uniformly in \(\varvec{\varepsilon }\)

Let f be the solution to linear transport equation (1)–(2). We define the Mth-order projection operator

Remark 4.2

Theorem 4.3 gives a uniformly in \(\varepsilon \) spectral convergence rate; thus, one can choose M independent of \(\varepsilon \), a very strong s-AP property. If the scattering is anisotropic, namely \(\sigma \) depends on \(\nu \), then one usually obtains a convergence rate that requires \(M\gg \varepsilon \) (see for example [12]). In such cases the proof of s-AP property is much harder, and one usually needs to use the diffusion limit, see [5] in the case of deterministic case and [11] in the random case.

5 The Full discretization

As pointed out in [15], by using the gPC-SG formulation, one obtains a vector version of the original deterministic transport equation. This enables one to use the deterministic AP scheme. In this paper, we adopt the AP scheme developed in [22] for gPC-SG system (16).

One of the most important and challenge problems for linear transport equation is the treatment of boundary conditions; here, we refer to the early work by Jin and Levermore [10] and more recent work by Lemou and Méhats [21] for their study of AP property and numerical treatment of physical boundary conditions.

We take a uniform grid \(x_i = ih, i = 0, 1, \ldots N\), where \(h=1/N\) is the grid size, and time steps \(t^n=n \Delta t\). \(\rho ^n_{i}\) is the approximation of \(\rho \) at the grid point \((x_i, t^n)\), while \(g^{n+1}_{i+\frac{1}{2}}\) is defined at a staggered grid \(x_{i+1/2} = (i+1/2)h, i = 0, \ldots N-1\).

where \(K =\frac{1}{3} \Sigma ^{-1}\). This is the fully discrete scheme for (17). Thus the scheme is stochastically AP as defined in [15].

We also notice that \([ \hat{g}^n_{i+\frac{1}{2}}]=0\) for every n which will be used later.

6 The uniform stability

One important property for an AP scheme is to have a stability condition independent of \(\varepsilon \), so one can take \(\Delta t \gg O(\varepsilon )\) when \(\varepsilon \) becomes small. In this section we prove such a result. The proof basically follows that of [23] for the deterministic problem.

For clarity in this section we assume \(\sigma ^a =S=0\). The main theoretical result about the stability is the following theorem:

Remark 6.1

Since the right-hand side of (87) has a lower bound when \(\varepsilon \rightarrow 0\) (and the lower bound being that of the stability condition of discrete diffusion equation (86)), the scheme is asymptotically stable and \(\Delta t\) remains finite even if \(\varepsilon \rightarrow 0\).

6.1 Notations and useful lemma

We give some useful notations for norms and inner products that are used in our analysis. For every grid function \(\mu =(\mu _i)_{i=0}^{N-1}\) define:

Now we give some notations for the finite difference operators that are used in scheme (85). For every grid function \(\phi =(\phi _{i+\frac{1}{2}})_{i\in {{\mathbb {Z}}}}\), we define the following one-sided difference operators:

We first recall some basic facts. For every grid functions \(\phi =(\phi _{i+\frac{1}{2}})_{i=0}^{N-1}, \psi =(\psi _{i+\frac{1}{2}})_{i=0}^{N-1}\), and \(\mu =(\mu _{i})_{i=0}^{N-1}\), one has (see [23]):

6.2 Energy estimates

Now we provide the details of the energy estimate. The proof is similar to that for deterministic problem in [23].

First, multiplying (85a) and (85b) by \({{\hat{\rho }}}^{n+1}\) and \(\varepsilon ^2 {\hat{g}}^{n+1}\), respectively. With the assumption that \(\sigma _i^a=0, {\hat{S}}_i=0\), and using the fact that \(\Sigma \ge \sigma _{\text {min}} Id\), one has

When \(\varepsilon \) is small, we use this as the reference solution, as it is accurate with an error of \(O(\varepsilon ^2)\). Hereafter we set \(\varepsilon =10^{-8}\). For large \(\varepsilon \) or in the case we cannot get an analytic solution, we will use the collocation method (see [27]) with the same time and spatial discretization to micro–macro system (11) as a comparison in the following examples. For more details about the collocation method, especially the s-AP property, we refer to the discussion in the work of Jin and Liu [12]. In addition, the standard 30-points Gauss–Legendre quadrature set is used for the velocity space to compute \(\rho \) in the following example.

To examine the accuracy, we use two error norms: the differences in the mean solutions and in the corresponding standard deviation, with \(\ell ^2\) norm in x:

In Fig. 1, we plot the errors in mean and standard deviation of the gPC numerical solutions at \(t = 0.01\) with different gPC orders. Three sets of results are included: solutions with \(\Delta x = 0.04\) (squares), \(\Delta x = 0.02\) (circles), \(\Delta x = 0.01\) (stars). We always use \(\Delta t = 0.0002/3\). One observes that the errors become smaller with finer mesh. One can see that the solutions decay rapidly in N and then saturate where spatial discretization error dominates. It is then obvious that the errors due to gPC expansion can be neglected at order \(M = 4\) even for \(\varepsilon = 10^{-8}\). The solution profiles of the mean and standard deviation are shown on the left and right of Fig. 2, respectively.

We also plot the profiles of the mean and standard deviation of the flux vf in Fig. 3. Here we observe good agreement among the gPC-Galerkin method, stochastic collocation method with 20 Gauss–Legendre quadrature points and analytical solution (103).

In Fig. 4, we examine the difference between the solution \(t = 0.01\) obtained by the 4th-order gPC method with \(\Delta x = 0.01, \Delta t = \Delta x^2/12\) and limiting analytical solution (103). As expected, we observe the differences become smaller as \(\varepsilon \) is smaller in a quadratic fashion, before the numerical errors become dominant.

7.2 Example 2: mixing regime

In this test, we still set \(\sigma = 2 + z\). We consider \(\varepsilon >0\) depending on the space variable in a wide range of mixing scales:

The reference solution is obtained using collocation method with 30 points. The parameters are set up as the following: the mesh size is \(\Delta x = 0.01\), and the corresponding t direction mesh size is \(\Delta t = \Delta x^2/3\). And we use the 5th-order gPC-Galerkin method to evolve the equation to different time \(t=0.005, t=0.01, t=0.05, t=0.1\). For the v integral, we use Legendre quadrature of 30 points.

Figure 6 shows the \(\ell ^2\) error of the mean and standard deviation with respect to the gPC order. We also see the fast (spectral) convergence of the method.

Fig. 6

Example 2 with initial data (111)–(112). The \(\ell ^2\) error of mean and standard deviation (dash line) with respect to gPC order

7.3 Example 3: random initial data

We then add randomness on the initial data (\(\sigma =2+z\) still random).

$$\begin{aligned} f(0,x,v,z) = f(0,x,v)+0.2z \end{aligned}$$

(113)

where f(x, v, 0) is the same as in (111). This time we set \(\Delta x=0.01\) and \(\Delta t = \Delta x^2/12\) and final time \(T=0.01\). First we test in the fluid limit regime \(\varepsilon =10^{-8}\) as shown in Fig. 7.

7.4 Example 4: random boundary data

We also test when \(\varepsilon =10^{-8}\) and \(\varepsilon =10\) as shown in Figs. 9 and 10, again, good agreements are observed between the gPC-SG solutions and the solutions by the collocation method.

where we set \(\sigma =4\) and \(z_1, z_2\) are both uniformly distributed in \((-1,1)\). The mean and standard deviation of the solution \(\rho \) at \(t=0.01\) obtained by the 5th-order gPC Galerkin with \(\Delta x = 0.025, \Delta t = 0.0002/3\) are shown in Fig. 11. We then use the high-order stochastic collocation method over 40\(\times \)40 Gauss–Legendre quadrature points to compute the reference mean and standard deviation of the solutions. In Fig. 12, we show the errors of the mean (solid lines) and standard deviation (dash lines) of \(\rho \) with respect to the order of gPC expansion. The fast spectral convergence of the errors can be clearly seen.

Errors of the mean (solid line) and standard deviation (dash line) of \(\rho \) with respect to gPC order, with the \(d=2\)-dimensional random input

8 Conclusions

In this paper we establish the uniform spectral accuracy in terms of the Knudsen number, which consequently allows us to justify the stochastic asymptotic-preserving property of the stochastic Galerkin method for the linear transport equation with random scattering coefficients. For the micro–macro decomposition-based fully discrete scheme we also prove a uniform stability result. These are the first uniform accuracy and stability results for the underlying problem.

It is expected that our uniform stability proof is useful for more general kinetic or transport equations, which is the subject of our future study.

Acknowledgements

Research was supported by NSF Grants DMS-1522184, DMS-1514826 and the NSF RNMS: KI-Net (DMS-1107291 and DMS-1107444). Shi Jin was also supported by NSFC Grant No. 91330203 and by the Office of the Vice Chancellor for Research and Graduate Education at the University of Wisconsin-Madison with funding from the Wisconsin Alumni Research Foundation.

Notes

In honor of the 70th birthday of Professor Björn Engquist

Declarations

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Authors’ Affiliations

(1)

Department of Mathematics, University of Wisconsin-Madison

(2)

Institute of Natural Sciences, School of Mathematical Science, MOE-LSEC and SHL-MAC, Shanghai Jiao Tong University