Uniform high-order convergence of multiscale finite element computation on a graded recursion for singular perturbation

1.   Introduction

Singular perturbation problems are ordinary/partial differential equations with perturbed parameters. For cases those are of a small parameter $\varepsilon$, the convective term dominates in the convection-diffusion model, and the boundary layer phenomenon involves. A remarkable property is that the equation order reduces as $\varepsilon$ approaches 0. At the same time, the solutions demonstrate either boundary or interior layer(s) in which their derivatives are extremely large. The singularly perturbed problems appear from mathematical modellings of scientific engineerings, such as flows with large Reynolds numbers, heat transfer with large Peclet numbers, turbulent interactions of waves and currents, mechanical and electrical system, chemical reaction, biotic species, etc.

Basically, the presence of boundary/interior layer(s) would create computational difficulties, and it is hard to reach the closed-form solutions with regard to elementary functions. It is known that traditional methods for conducting singular perturbations are inefficient because to resolve layers precisely an extremely fine mesh covering the global domain is required. This sufficiently fine discretization will lead to systems that are prohibitively expensive for solving, especially in high dimensions. The a priori knowledge of the layer location is need to take into account. Stability of the numerical approach is frequently another concern. In general, two popular approaches are processed for designing efficient algorithms: one is based on appropriately refined layer-adapted meshes, and the other is based on locally exact or fitted schemes. A characteristic property of those potential methods is a uniform convergence in terms of the perturbed parameter; it is independent of whatever the parameter value is to be validated the stability of numerical solutions.

Many kinds of literature have focused on fitted numerical methods for the singularly perturbed problems [14], such as finite difference method (FDM) [6,10,16,17], finite element method (FEM) [2,3,4,5,7,18], finite volume method (FVM) [19], discontinuous Galerkin method [8,20], reproducing kernel method [1], collocation method, spectral method, etc. Shishkin [16] proposes a finite difference method on the a priori adapted meshes, it converges almost $\varepsilon$-uniformly for a singularly perturbed parabolic equation. In [6] a one dimensional parabolic systems of reaction-diffusion are studied, components splitting are used by the central finite difference scheme to discretize in space and by the implicit Euler scheme to discretize in time. Chen and Xu [5] present a streamline diffusion finite element method and obtain an optimal maximum norm with stability and accuracy on the adapted grid. [3] gains the existence and uniqueness of the solution of a nonlinear elliptic problem, it is verified to be a parameter robust method for solving nonlinear cases. A finite volume discretization is provided for the perturbed equation in [19], the robust a posteriori error estimate is based on flux reconstruction with anisotropic meshes. Du and Chung [8] investigate an adaptive staggered discontinuous Galerkin method for a steady state convection-diffusion equation, and the a posteriori error estimator is derived under the boundedness assumption on $h/\varepsilon$.

The multiscale finite element method (MsFEM) is an updated approach for this type of problems. We [11] consider a two dimensional reaction-diffusion model and propose an adapted Petrov-Galerkin MsFEM, it provides more flexibility in functional space and eliminates the multiscale resonance automatically to reach a stable convergence. In [12] the MsFEM is applied for a singularly perturbed convection-diffusion problem with the Robin boundary, through the reduced mapping matrix and special interpolations we resolve the rapid oscillation of boundary layers.

Motivated by these works, the goal of this paper is to present a high accuracy and a high-order convergence of the multiscale finite element computation for singular perturbations, it will resolve the solution's singularity with an underlying mesh adapted to its singular nature. It is accomplished by a reduced multiscale strategy, with the help of multiscale basis functions to capture the microscopic information locally and to save the macroscopic capacity globally.

The paper proceeds as follow. In Section 2, a singularly perturbed convection-diffusion equation with a small parameter is introduced, and the classical finite element approximation is used for it first. Then in Section 3 we propose a reduced multiscale strategy for solving the model. A recursion of graded mesh is constructed, and the numerical analysis of its multiscale error estimate is provided. Numerical results are offered to testify the behaviors of the multiscale simulation in Section 4. We end up with conclusions in Section 5.

2.   Model problem and Galerkin method

We consider a convection-diffusion equation

where $0<\varepsilon<<1$ is a small parameter. $u(x)$ is a solution, and $b(x),c(x)$ are sufficiently smooth coefficients with properties

It is known that the solution has a bound

its $k$-th order of derivative $k = 0,1,\cdots,q$ (for a prescribed $q\in \mathbb{N}$) and $exp(x,\beta,\varepsilon) = \exp^{-{\beta x\over\varepsilon}}+\exp^{-{\beta (1-x)\over\varepsilon}}$. Throughout in this context, $C$ is a generic independent constant. Note that in this form the solution may have boundary layers at $x = 0$ and/or $x = 1$. Convection-diffusion equation of this type arises in linearised versions of the Navier-Stokes equation.

Standard Sobolev spaces $W^{k,p}(\Omega), L^p(\Omega) = W^{0,p}(\Omega), H^k(\Omega) = W^{k,2}(\Omega), H_0^k(\Omega)$ are used, and $(\cdot,\cdot)$ is used for the $L^2(\Omega)$ inner product. Related seminorm and norm are $|\cdot|_{k,p,\Omega}$ and $\|\cdot\|_{k,p,\Omega}$ on $W^{k,p}(\Omega)$. In the case $p = 2$, we simplify $|\cdot|_{k,\Omega} = |\cdot|_{k,2,\Omega}$ and $\|\cdot\|_{k,\Omega} = \|\cdot\|_{k,2,\Omega}$. A $\varepsilon$-weighted energy norm is introduced as

The original motivation of this paper is to practice the multi-scale nature through a scale decomposition of $u(x)$ into two parts, a smooth $S(x)$ and a singular $E(x)$, i.e.,

Set $LS = f$, $LE = 0$ to satisfy $Lu = f$, and their $k$-th order of derivatives have

The weak formulation of (1) is based on a bilinear form

and an inner product

to satisfy

Denote $K$ as an element and $\mathcal{K}^h$ as a mesh partition (for any $K \in \mathcal{K}^h$), and $V^h$ is denoted as

where $P_m(K)$ denotes a space of polynomials of degree less than or equal to $m\in \mathbb{N}$. The classical Galerkin finite element method is to find $u_g\in V^h$ such that

here $u_g$ is the Galerkin FEM solution.

It is known that for cases of a large parameter $\varepsilon$ in (1), there is no singular perturbation so that traditional methods such as FEM works well. However, for cases of a small parameter, the boundary layers phenomena are quite troublesome, the solution or its derivative would exhibit jumps or huge steepness. Our goal is to obtain the most accurate and convergent solution (with possibly least computations), and the result does not rely on the perturbed parameter $\varepsilon$. In this way, we are committed to finding powerful methods to address this situation. In the paper, we would like to present a novel strategy of the multiscale finite element method on an adaptive coarse mesh, which is capable of providing enough accuracy with just moderate computer resources, and it ensures uniform convergence for singular perturbations.

3.   Multiscale adaption and error estimate

3.1.   Multiscale finite element method

Being different from the Galerkin FEM, the idea of a multiscale approach is to split the multiple solution space into resolvable and unresolvable scales. It is processed by using the piecewise polynomials space $V^h$ to represent the resolvable scale, and using a projection operator $P:U\rightarrow V^h$ to represent the unresolvable scale,

We reformulate the weak form (8) as to find $u_1\in V^h$, $u_2\in U^h$ such that

In order to eliminate the unresolvable scales in (14), $B(u_1)$ and $F(f)$ are defined as the following solutions to find $B(u_1)\in U^h$, $F(f)\in U^h$ such that

Thus the solution of (14) becomes $u_2 = F(f)-u_1 = F(f)+B(u_1)$. And with the elimination of resolvable scales in (13) we find $u_1\in V^h$ such that

As we can see that space $U^h$ is to represent the unresolvable scales, it will be enriched with the multiscale basis functions $\varphi_i$ (in addition to the finite element basis functions $\psi_i$), that is

These multiscale bases $\varphi_i$ have abilities to reflect the microscopic information of the macroscopic problem (1), and they may capture the local perturbation in boundary layers automatically. It is fulfilled through conducting following local problems in the finite element scheme.

We solve the multiscale basis functions on each coarse element $K$,

The boundary condition Kronecker $\delta_{ij}$ is refined as $\varphi_i(x_j) = 0$ when $i\neq j$ and $\varphi_i(x_j) = 1$ when $i = j$. As for the local problem (19) and the original problem (1) they have the identical differential operator $L$, we solve (19) by applying the FE scheme (with a fine partition $M$) in coarse elements. Consequently on the coarse level, the microscopic information of boundary layers is captured through the multiscale basis functions $\varphi_i$. We should point out that it is quite different from the classical FEM, which is computed on the fine level and would consume more computational costs.

As a result, we are driven to find $u_h\in U$ in a weak form such that

here $u_h$ is the MsFEM solution. The strength of multiscale computation comes from its enrichment of the multiscale space.

In our codes the available multiscale basis functions are deposited in a mapping matrix $R$, whose every row and column in each coarse element $K$ has local microscopic information. In result, a (reduced) global matrix $A_{\text {ms}} = R*A*R^T$ and a (reduced) global vector $F_{\text {ms}} = R*F$ are stored on the coarse level. Here $A, F$ are the FEM (full) stiffness matrix and (full) right vector on the fine level. Through the multiscale finite element scheme, we assemble the global equations to solve $A_{\text {ms}}u_h = F_{\text {ms}}$ for $u_h$ on the coarse level.

It should be noted that when we apply the FEM on a very fine mesh, the global system $Au_g = F$ is massive. On the contrary, if we apply the MsFEM on a relatively coarse mesh to obtain the MsFEM solution (with enough accuracy), the corresponding system turns to $A_{\text {ms}}u_h = F_{\text {ms}}$, which are quite small-scaled. This merit shows the advantage of our reduced MsFEM, especially for high dimension problems. As a consequence, this multiscale computation is simply of limited costs to achieve an accurate solution. We are pleased to show that the MsFEM on an adaptive coarse mesh is an optimal strategy in the following.

3.2.   Mesh adaption

We provide sorts of mesh partitions $\mathcal{K}^h$ in this subsection. Suppose the model (1) has a large parameter of $\varepsilon$, there is no boundary layer and traditional methods behave well. It is validated that usage of Uniform mesh is sufficient (whose size $h$ is a constant ${1/N}$, here $N$ is a partition number). However, with a small $\varepsilon$ the model will produce boundary layers whose width is $O(\tau) = O(\varepsilon\text{ln} N)$. As for this situation, even though with a very fine partition $N$ the traditional method would like to perform a fake result.

For a small parameter $\varepsilon$, the domain $\Omega$ is divided into smooth and singular parts concerning a transition point $\tau = \text{min}\{{1\over2},{\varepsilon\text{ln} N\over \beta}\}$. We take the a priori estimate for the model problem and figure out the appropriate locations of boundary layers; then several non-Uniform meshes are to be employed. Here the total partition number $N$ is kept fixed, as for in a sub-domain it is refined to approach the singular part, while in another sub-domain it is coarsened to approach the smooth part. These sorts of $h$-adapted strategy are offered as follows.

First we provide Shishkin's idea [13,17]. Suppose that boundary layers appear near the right side $x = 1$, we divide $\bar\Omega = [0,1]$ into $\Omega_1 = [0,1-\tau]$ and $\Omega_2 = [1-\tau,1]$, both sub-domains are partitioned by $N/2$ elements. Shishkin grid node is

Next, a Bakhvalov mesh in [12] is presented to be another adapted layer. Suppose the boundary layer is near $x = 1$, its mesh distribution is determined by $e^{-{\beta (1-x_i)\over\varepsilon}} = Ai+B$ and $x_{N\over 2} = 1-\tau$, $x_N = 1$. In this way, Bakhvalov grid node is

Finally, a more flexible improvement is a Graded mesh provided in [9,15], which is a highly anisotropic non-Uniform. It is based on a recursive iteration, the Graded grid node is

where $0<h<1, \sigma>0$ are initial constants for the mesh generation. It is taken $h = 0.5, \sigma = 1$ in the following experiments. Be aware that $N$ is no longer fixed number any more and it satisfies $x_{N-1}<1$ and $(1+\sigma h)x_{N-1}\geq 1$. In this case plenty of grid nodes will be concentrated on the left boundary $x = 0$.

However, if the boundary layer appears on the right boundary $x = 1$, in our Matlab codes a command $ones-fliplr$ is used for flipping arrays from left to right. As a result, these adaptive nodes will be concentrated on the other side contrarily. Another thing to remark is about the generated partition number $N$, it could be even or odd positive integer according to the recursive algorithm. These outputs would be testified in the following numerical result section 4.

For more general situations, suppose that singular perturbations dominate on both sides $x = 0$ and $x = 1$, we provide a mesh modification as

here $N_0$ is the size of the generated $N$ in (23). It will produce a quasi-symmetric mesh partition. The Graded mesh adaption could be applied to efficiently address the solutions that have both boundary layers.

In the paper we set $\alpha = k+1$. A transition point is supposed to define as $\tau = \alpha\varepsilon\text{ln} N$ here. Near the left point $x = 0$ we set $x_{N\over 2} = \tau$ and in the sub-domain $[0,\tau]$ a grid generating function $g$ is introduced, and it has $g(0) = 0,g({1\over 2}) = \log N$. Then for grid nodes we have $x_i = \alpha\varepsilon g({i\over N})$. Also, another auxiliary grid characterising function $\omega$ is defined as

As we can see that the Graded is not a pre-determined mesh, it relies on the recursive iteration, an initial parameter $h$ and a perturbed parameter $\varepsilon$ in (23) and (24). Once the Graded nodes $x_i$ are available, we denote $h_i = x_i-x_{i-1}$ as a varying mesh size in the $i$-th coarse element $K$. We confirm that the mesh size is monotone increasing/decreasing, and this grid is prepared for the multiscale finite element computation on the coarse level. The advantage of Graded mesh is that it may adaptively adjust to the local quantity of original perturbed problem, and it is shown to be an optimal mesh for capturing boundary layer phenomena accurately and economically in the multiscale computation.

3.3.   Multiscale error estimate

Now we are ready to prove the error estimate for the multiscale approximation on the Graded mesh, an error estimate is to be derived between the exact solution $u$ and the multiscale solution $u_h$. Associated with the above coarse mesh $\mathcal{K}^h$, we investigate the enriched multiscale space $U$ in (11).

Lemma 3.1. [9] Let $u$ and $u_I$ be the solution and its interpolation of (8), respectively. Then

Lemma 3.2. Let $u = S+E$ in (4) be a solution decomposition into a smooth part and a singular part, $S_I$ and $E_I$ their interpolants, repectively. Then

Proof. Recall that $\alpha = k+1$, $\tau = \alpha\varepsilon\text{ln} N$ and $\Omega_1 = [0,1-\tau]$, $\Omega_2 = [1-\tau,1]$, we prove (28) first.

As for (30), since

then we get $\|E-E_I\|_{0,\Omega_2} \leq C\varepsilon^{1\over 2}N^{-\alpha} \max|\omega'|^\alpha$. The proofs for (27) and (29) are omitted.

Theorem 3.3. Let $u$ and $u_h$ be the solution of (8) and the multiscale solution of (20), respectively. Then

Proof. From the triangular inequality, it is known that

For the first term on the right, from Lemma 3.2 we have

then $\|u-u_I\|_{0,\Omega} \leq C\varepsilon^{1\over 2}N^{-\alpha}\max|\omega'|^\alpha$. With respect to the definition of the energy norm in (3) and from Lemma 3.1, we get

For the second term, let $e = u_I-u_h \in U^h$, using the error $a(u-u_h,e) = 0$ and the bilinear form in (7) we consider

in the last step, the fourth term attracts a concern. Since

we have $|(b(E-E_I),e')|_{\Omega} \leq CN^{-\alpha}\max|\omega'|^\alpha\|e\|_\varepsilon$. Divided by $\|e\|_\varepsilon$ on both sides, and synthesising the estimates the error for $\|u_I-u_h\|_\varepsilon$ is obtained.

Therefore the proof completes.

Remark 1. It is verified in Theorem 3.3 with $\alpha = k+1$ when $k = 1$ linear polynomials are used in our multiscale finite element computation, the second-order of convergence is guaranteed. Besides, with the contribution of the adaptive grid function $|\omega'|^\alpha$, the novel method may have the potency to acquire higher convergent order.

4.   Numerical experiment

In this section, singularly perturbed convection-diffusion models are solved on different meshes by the finite element method and by our reduced multiscale method, respectively. The ability and superiority of the multiscale computation combined with the Graded mesh will be presented.

The finite element method and the multiscale finite element method are applied on different meshes, to get the numerical solutions $u_g$ and $u_h$ respectively, which are compared with the exact solution $u$ to testify their accuracy. For convenience, an initial abbreviation of meshes is represented Uniform as (U), Shishkin as (S), Bakhvalov as (B) and Graded as (G). Computations are executed on the mesh partition number $N$. As a result, two methods introduce their discrete system $Au_g = F$ and $A_{\text {ms}}u_h = F_{\text {ms}}$ for solving, which are of large scale and small scale, respectively.

Numerical errors and convergence orders are tested by

They are criteria to judge the corresponding abilities of numerical methods. Here $u_N$ may be the classical Galerkin solution $u_g$ on the fine level, or the multiscale solution $u_h$ on the coarse level. In addition, $E(2N)$ on the Graded mesh is not from a fixed $2N$ exactly but is from its close approximation.

Example 1. Given in (1) with $b(x) = 2$, $c(x) = 1$ and an exact solution is

We have $u_0 = 0$, $u_1 = 1$ and an identical right side $f(x) = u(x)$. Note that with a small parameter $\varepsilon$ in this example, $u(x)$ presents an abrupt jump near the right side $x = 1$.

We apply the FEM and MsFEM on four types of meshes shown in Figure 1, and detailed results for a small perturbed parameter $\varepsilon = 2^{-20}$ are listed in Table 1. It is known that even though with a very fine partition number $N = 10240$, the FEM(U) behaves poorly and it has a slow convergence. With the help of Shishkin and Bakhvalov meshes, the FEM provides a much better accuracy in the refinement history of $N$. FEM(S) realizes a clear 1.75-order convergence. The accuracy of FEM(B) is based on large partition numbers $N$ too, but we are not quite sure about its somehow fast or slow convergence in this example. That we concern most is the MsFEM on a Graded coarse mesh, whose coarse partition is an automatical product, it is even or odd at sometime. Its convergent order of the energy norm is high and is kept stable from the second to the third, which is in conformity with Theorem 3.3. It should be pointed out that the MsFEM(G) with coarse partitions is matching up or even better than the FEM(B) with fine partitions, consequently the advantage and superiority of MsFEM(G) are quite promising.

It is observed in Figure 2 with partition number $N = 160$ the numerical solution of FEM(U) is far wrong from the exact solution. While those solutions of FEM(S) and FEM(B) are well approximated, on a relatively fine mesh $N = 160$. What is more encouraging is that even with a coarse partition $N = 37$, the MsFEM(G) captures well. We also find that the grid partition is concentrated to the location of the boundary layer, while in the smooth part it is sparse partitioned.

In Figure 3, corresponding discrete errors are plotted. The error (green line) from FEM(U) is quite annoying, which is from 0 to 3.5. Those maximum errors of FEM(S) and FEM(B) are reduced to $6*10^{-4}$ and $1*10^{-3}$, respectively. Moreover, we observe the shape of the boundary layer on the right side in FEM(S), while FEM(B) has almost an equal-distributed error. Besides, with a coarse partition $N = 37$, the maximum error of MsFEM(G) is $2*10^{-3}$. The large error (red line) is detected only at the near location of the boundary layer. Except this layer, there is an almost negligible error on this adaptive coarse mesh.

To further illustrate the uniform convergence of MsFEM on the Graded mesh, for different small $\varepsilon$ we employ this approach to reach for the same magnitude of energy norm error, see Table 2 and Figure 4. The high-order (between the second and the third) convergence is validated again, and the coarse partition $N$ is shown as a linear growth, which is independent of whatever $\varepsilon$. It is reliable that for cases of smaller $\varepsilon$ there is a relatively greater $N$. At the same time, the coarse partition number in the MsFEM is still much less than the fine one in the FEM.

Example 2. Given $b(x) = 1$, $c(x) = 1+\varepsilon$, the exact solution is available as

Then $u_0 = 1+\exp^{-{1+\varepsilon\over\varepsilon}}$, $u_1 = 1+\exp^{-1}$ and a homogeneous right side $f(x) = 0$. In this way, the solution with a small $\varepsilon$ will arise boundary layers on both sides. We apply the FEM on the different fine meshes and the MsFEM on the Graded coarse mesh (24), respectively, to solve this model with $\varepsilon = 2^{-30}$.

The poor result from the FEM(U) is omitted in Table 3. With the refinement history of a fine partition $N$, the accuracy of FEM(S) and FEM(B) is acceptable, their convergence order is less than the second. We find that the coarse partition $N$ in MsFEM(G) is almost one quarter to the fine partition in FEM, while the reduced multiscale method preserves the highest accuracy and it has a superconvergence greater than the second-order in Theorem 3.3.

In Figure 5 and Figure 6, FEM results on a fine partition $N = 320$ and MsFEM results on a coarse partition $N = 106$ are plotted. It is distinct that 106 coarse Graded nodes are automatically distributed to both sides densely, and the boundary layer errors are greatly reduced to the magnitude of $4.5*10^{-4}$. The advantage of the reduced multiscale approach is validated again.

In the experiment results, it is demonstrated that by our reduced multiscale computation on the Graded mesh, the numerical solutions are independent of the perturbed parameters $\varepsilon$, and they are numerically stable and uniformly superconvergent.

5.   Conclusion

In this paper, a novel strategy of the reduced multiscale is proposed to solve a convection-diffusion model with small perturbed parameters. Making use of the multiscale nature of the problem, we decompose it into singular scales and smooth ones. A Graded mesh is built for capturing the local microscopic information among scales. It is a type of adaptive grid, and its unfixed partition number is generated according to the recursive iteration. We present the mathematical foundation on error estimates, and in numerical experiments, the strength of our approach has been shown. The multiscale computation is efficient and robust for no matter of the perturbed parameters; it offers high accuracy and uniform superconvergence with merely moderate computation resources for the singularly perturbed problem.

Acknowledgments

We appreciate referees and editors for their insightful comments and suggestions.