A high-order convergence analysis for semi-Lagrangian scheme of the Burgers' equation

1.   Introduction

One of the most popular techniques for solving advection-type equations is the backward semi-Lagrangian method (BSLM), which has the advantage over conventional Eulerian approaches of using a larger temporal step size while avoiding the mesh deformation of pure Lagrangian methods. Due to these attributes, BSLMs have been much studied in a vast range of models since their introduction in the early 1980s ^[29] such as weather forecasting ^[31,33], oceanography ^[29,38], engineering control ^[16], and various important fluid dynamics applications (^{[3,4,5,6,15,23,30,32,36,37]} and the references therein). In addition, numerous research efforts have been recently made to develop more efficient BSLM algorithms in computational aspects ^[10,21,24] as well as in mass conservation and multi-dimension ^[11,12,13].

Their convergence analyses have been conducted in various problems including a second-order BSLM in linear advection–diffusion equations ^[19] and nonlinear advection-type equations ^[27] using the conservation property; second-order finite element BSLMs in the Navier–Stokes equations given a priori estimates ^[8,9], and second-order finite difference BSLMs in advection–diffusion equations with the Dirichlet boundary condition ^[22]. Despite these convergence analysis studies, there remained unaddressed nonlinear advection-type problems of finite difference BSLMs with a high-order spatial discretization due to the complexity and difficulties of understanding the structure for the Cauchy problem represented by a characteristic curve.

The BSLM offers a way of dealing with the nonlinear advection term by splitting the model equation into a Cauchy problem and a total differential equation along characteristic curves. Even though this nonlinearity appears to vanish in the Lagrangian procedure, it is expressed as the velocity of the characteristic curve, which creates the uncertainty of its departure points. This uncertainty makes the convergence analysis more challenging because the truncation error of the departure points of Cauchy problem contains the errors of the solution from the previous temporal steps in the time discretization (see (2.16) and (2.17) for example). Analyzing the error arising from the diffusion term along the characteristic curves is another difficulty in high-order spatial discretizations. The second-order spatial discretization gives a positive-definite matrix for the diffusion term satisfying the Dirichlet boundary condition ^[22], which allows the Cholesky decomposition to split the errors into consecutive temporal steps. The other issue of the convergence analysis is in the temporal evolution of the truncation errors and bounds of the error from each time step. For the Dirichlet problem with the second-order spatial discretization, the discrete Poincaré inequality allows to evaluate the error at the final step using the spatial finite difference of the error. However, the mentioned techniques are no longer possible to analyze the BSLM with high-order spatial discretizations in the Dirichlet problem, which is the concern of this paper.

This work is a part of the sequential research project on the convergence analysis of various proposed finite difference BSLMs due to ^{[25,26,27,28]} as well as a subsequent result of ^[22]. The periodic problem is permitted to discretize a diffusion term with a symmetric high-order spatial approximation. From this symmetry, as the second-order analysis, the Cholesky decomposition allows that the errors related the diffusion term can be split into the ones at two consecutive temporal steps, (see Lemma 4). For the analysis of the finite difference BSLM with a high-order spatial discretization, we apply the Cholesky decomposition of the positive definite matrix for a discretization of the Laplacian in the periodic problem. Instead of using the discrete Poincaré inequality which is no longer allowed in the periodic problem, we apply the telescoping summation technique in order to obtain the bound of the error at the final step using the finite difference of error in time. This not only allows a high-order approximation analysis for spatial discretization but also provides the general convergence analysis without boundary restriction for nonlinear advection–diffusion-type equations. More precisely, we consider Burgers' equation with the periodic boundary condition, which is a prototype of nonlinear advection-diffusion equations, as follows:

where $u(t, x)$ represents the particle velocity, $\nu$ stands for the kinematic viscosity, and $u^0(x)$ is a given smooth function. To illustrate the BSLM for the model problem (1.1), we assume that the temporal steps are uniformly distributed throughout the domain $[0, T]$; that is, the grid points are formed by calculating

with the temporal step size $h = T/N$. For each temporal step $t_n$, let ${\chi}({x}, t_{n+1};t)$ be the characteristic curve arriving the spatial point $x$ at time $t_{n+1}$ that satisfies the nonlinear Cauchy problem described as

Combining (1.1) and (1.2), the total derivative of $u (t, { \chi}({x}, t_{n+1};t))$ along the characteristic curve ${ \chi}({x}, t_{n+1};t)$ provides the following diffusion equation given by

Notice that the characteristic curve ${ \chi}({x}, t_{n+1};t)$ satisfying the Cauchy problem (1.2) is used to evolve the solution $u$ of (1.3) along the characteristic curve, while the solution $u$ is used as a velocity for the characteristic curves with a given initial value at the fixed grid points. This self-consistency requirement of the system generates the main aforementioned hindrance to analyzing the constructed method. The nonlinear Cauchy problem can be solved by several approaches, including iteration methods ^[1,34,36], iteration-free methods ^[26,27], and so on. The solution of (1.3) is usually obtained by combining three procedures: a backward finite difference method for the total derivative, an interpolation technique for the solution at off-grid departure points, and a solver for the diffusion term, for instance, a finite difference method ^[14], a finite element method ^[3,4,14], a spectral finite element method ^[36], or a discontinuous Galerkin finite element method ^[35].

In this paper, we present a concrete convergence analysis of the finite difference BSLM that adopts an iteration-free error correction method (ECM) ^[26,27] to solve the Cauchy problem (1.2), the second-order backward difference formula (BDF2) for the total derivative, and the fourth-order central finite difference scheme for the diffusion term in the discretization of the Eq (1.3). The substantial contribution of this study is to provide not only the comprehensive analysis of convergence and stability of a high-order spatial discretization BSLM but also the experiments to support the theoretical results. For the prerequisite process, we review the BSLM, which is modified from ^[26,27], to determine the numerical solution of the model problem while providing analysis of the local truncation errors of departure points obtained by the ECM, the temporal and spatial discretizations from the finite difference methods, BDF2 and the fourth-order central difference formula. Using the relations between these truncation errors and the boundedness properties of the interpolation, we provide the error equation of the proposed scheme at each time step with the help of the symmetric discretization for the diffusion term. Finally, after taking the telescoping summation of the error equations through all the temporal steps, we obtain the error of the numerical solution at the final time (see Theorem 3.6). The resulting order of convergence $\mathcal{O}(h^2+\triangle x^4+ \triangle x^{p+1}/h)$ with respect to the discrete $\ell_2$-norm supports convergence analyses for linear advection-type equations ^[19] and finite element BSLMs ^[8,9], where $p$ represents the degree of the interpolation. The result also validates the claim shown by the numerical results, including the non-monotonic property ^[17] of $\mathcal{O}(\triangle x^{p+1}/h)$. We further investigate the stability of the proposed scheme by presenting that the discrete norm of the numerical solution is uniformly bounded without any restriction of the temporal step size $h$, which means that the proposed BSLM is unconditionally stable (see Theorem 3.7).

The paper is organized as follows. The BSLM, comprising the ECM and finite difference schemes for the discretization of Eqs (1.2) and (1.3), is reviewed and local truncation errors are analyzed in Section 2. Based on the error equation, the main result of the paper, the convergence and stability analysis of the BSLM for Burgers' equation, are presented in Section 3. In Section 4, numerical experiments are presented to support the theoretical analysis. Finally, Section 5 presents our conclusions and further discussion. Additionally, in Appendix, we review the Lagrange interpolation formula and its properties.

2.   Calculation of departure points and discretization of BSLM

For the derivation of the BSLM, we review the ECM for solving the Cauchy problem in Subsection 2.1 and the fourth-order spatial discretization for the Helmholtz equation in Subsection 2.3. Especially, Subsection 2.2 is devoted to the estimation of a concrete bound for the truncation errors of its departure points.

2.1.   Review of ECM

We start this subsection with a concise review of the ECM for solving the nonlinear self-consistent Cauchy problem given by

where ${ \chi}_{j}(t): = { \chi}_{j} ({x}_{j}, t_{n+1};t)$ is the characteristic curve that arrives at the grid point ${x}_{j} \; (1\leq j \leq {J})$ at a fixed temporal time step $t_{n+1}$ and $u$ is the solution of the problem (1.1). Here, we assume that spatial grids ${x}_{j} \; (1\leq j \leq {J})$ are uniformly spaced in the domain $I: = [a, b]$ as follows:

The ECM then focuses on finding the departure points $\chi_j(t_{n-k}) = \chi(x_j, t_{n+1};t_{n-k})$, $k = 0, 1$, which will be used for the BSLM. For notational simplicity, we denote ${u}^m({x}): = {u}(t_m, {x})$, ${u}^{m}_{j}: = {u}^{m}(x_{j})$, and $\mathbf{u}^m: = \left[{u}^{m}_{1}, \ldots, {u}^{m}_{J}\right]^{\mathtt T}$. Furthermore, we assume that approximations $\mathbf{U}^m: = \left[{U}^{m}_{1}, \ldots, {U}^{m}_{J}\right]^{\mathtt T}$ for the solutions $\mathbf{u}^m$ ($m\leq n$) are already calculated at all spatial grids ${x}_{j} \; (1\leq j \leq {J})$ in the domain $I$.

Let ${{\phi}} (t)$ be a perturbed characteristic curve given by an Euler polygon $y(t)$ such that

Then, by the Taylor expansion of $\chi_j(t)$ about $t = t_{n+1}$ one can see that

for some $\xi_1\in (t, t_{n+1})$. Again, the Taylor expansion of $u^{n+1}_j$ about $t = t_n$ implies that

where $\xi_2$ belongs to the interval between $t_{n}$ and $t_{n+1}$ and a constant $C$ depends only on the bounds of $u$, $u_t$ and $u_x$. Also, by (2.1)–(2.4), the derivative $\phi'(t)$ satisfies

where $T_1$ is bounded by

for a constant $C$ depending only on $u, u_t, u_x, u_{xx},$ and $u_{xt}$. We now use the bound (2.6) and the mid-point quadrature rule for integrating (2.5) over $[t_{n-1}, t_{n+1}]$ to obtain

The property that ${\phi}(t_{n}) = \frac{1}{2} \left({\phi}(t_{n+1})+{ {\phi}}(t_{n-1}) \right) +h^2\phi''(\xi_3)$ for some $\xi_3\in(t_{n-1}, t_{n+1})$ with ${\phi}(t_{n+1}) = 0$ yields

Since ${y}(t_n)$ is not a grid point in the usual sense, after truncating the asymptotic term of (2.7) the ECM uses the Lagrange interpolation to evaluate $u^n_x\left({y}(t_n)\right)$ and ${u}^n\left(y(t_n)\right)$. We then define the approximation $\phi^{n-1}$ of $\phi(t_{n-1})$ by

where ${\mathfrak L}$ and ${\widetilde D}_{x}{\mathfrak L}$ represent the Lagrangian interpolation polynomial of degree $p$ and its piecewise average derivative at the grid points, respectively (see (A.1) for details). Here, the approximation $\mathrm y^n$ of the Euler polygon $y(t)$ at $t = t_n$ satisfies

where $e^n_j$ is the error of $\mathbf{u}^n$ at $x_j$ such as

Instead of solving (2.1) directly, the iteration-free ECM uses the approximation $\phi^{n-1}$ of (2.8) and the relation (2.3) to obtain an approximation ${ \chi}_{j}^{n-1}$ of $\chi_j(t_{n-1})$ defined by

In addition, to find an approximation ${ \chi}^{n}_{j}$ for $\chi_j(t_n)$ we use the Taylor expansion of ${ \chi}_{j}(t)$ at $t_{n-1}$, given by

By using ${ \chi}^{n-1}_{j}$ and ${\mathfrak L} { \mathbf{U}}^{n-1}({ \chi}_{j}^{n-1})$ and (2.12), we finally define ${ \chi}^{n}_{j}$ by

2.2.   Estimation of truncation errors for departure points

In this subsection, we focus on estimating the truncation errors ${\rm{ \mathsf{ τ} }}^{n-k}_j$ of the approximations $\chi_j^{n-k}$ for two departure points ${ \chi}_j(t_{n-k})$ defined by

From (2.7) and (2.8), the approximation error of the perturbed characteristic curve $\phi(t)$ given by (2.4) satisfies

Therefore, Eqs (2.3), (2.9), (2.11), and (2.15) make up the truncation error for the departure points

where ${\epsilon}_{\rm{ \mathsf{ τ} }}$ is given by

Moreover, Eqs (2.12) and (2.13) imply that the truncation error at $t = t_n$ satisfies

To estimate the truncation errors of ${\rm{ \mathsf{ τ} }}^{n-k}_j$, we introduce ${\delta}_x$ that is a first-order backward spatial finite difference operator given by

Lemma 2.1. Assume that the Euler polygon and its approximation ${ \mathrm y}^{n}$ and $y(t_{n})$ belong to the same subinterval $I_j$ at $t = t_{n}$ for each $j$. Then we have

where $C$ and $C'$ depend only on the degree $p$ of the interpolation and the bounds of the first and second partial derivatives of $u$.

Proof. Using $\mathbf{e}^n$ in (2.10), we split the quantity $|{\mathfrak L} \mathbf{U}^n({ \mathrm y}^n) - u^n({y}(t_n))|$ into three terms such as

The definition (A.1) of interpolation and its associated property (A.3), respectively, provide that

where $C$ is independent of the indices $j_\ell$ and $j$. Also, the mean value theorem and the property (2.9) give $I_3\leq Ch|e^n_j|$ for some constant $C$ that depends only on the bounds of the partial derivatives of $u$, which implies that the first inequality holds. To estimate $\widetilde{D}_x {\mathfrak L} \mathbf{U}^n({ \mathrm y}^n) - u^n_x\left({y}(t_n)\right)$, we assume that $\mathrm y^n$ is in the interior of $I_j$ and split it into three terms. By virtue of Lemma A.1, (2.9), the mean value theorem, and the interpolation error (A.5), we have

where $C$ depends only on the bounds of the partial derivatives of $u$. This establishes the second inequality. The estimates of the previous lemma give the invertibility of $1+h{\widetilde D}_{x}({\mathfrak L} \mathbf{U}^n)({ \mathrm y}^n)$ as follows.

Corollary 2.2. Assume that for $n < N$, $\mathbf{e}^n$ and $\delta_x \mathbf{e}^n$ are sufficiently small. Then, there exists a temporal step size $h < 1$ such that $1+h{\widetilde D}_{x}{\mathfrak L} \mathbf{U}^n({ \mathrm y}^n)$ is invertible and

Proof. To show the invertibility of $\big|{\widetilde D}_{x}{\mathfrak L} \mathbf{U}^n({ \mathrm y}^n)\big|$ for a small $h < 1$, it is enough to show that it is uniformly bounded. This follows easily from the second result in Lemma 2.1 by considering

Combining the outcomes of Lemma 2.1 and Corollary 2.2 yields the estimates for the truncation errors of the departure points at the previous two steps:

Corollary 2.3. Assume that both ${ y}(t_{n-k})$ and ${ \mathrm y}^{n-k}$, the Euler polygon and its approximation as defined by (2.3), (2.8), and (2.11), belong to the same $j^{th}$ subinterval $I_{j}$ for $k = 0, 1$. Then the truncation error ${{\rm{ \mathsf{ τ} }}}^{n-k}_{j}$ $(k = 0, 1)$ of the departure points in (2.14) can be estimated as follows:

for some constant $C$ depending on $p$ and on the bounds of $u, u_t$ and $u_x$.

Proof. Note that the quantity $\epsilon_{\rm{ \mathsf{ τ} }}$ defined in (2.17) can be estimated from the bounds of (2.4) and Lemma 2.1. Thus, using the uniform bound of Corollary 2.2 yields the desired order of the estimate for ${{\rm{ \mathsf{ τ} }}}^{n-1}_{j}$. By applying the result for ${{\rm{ \mathsf{ τ} }}}^{n-1}_{j}$ in (2.18) and Lemma 2.1, we can also have the bound of ${{\rm{ \mathsf{ τ} }}}^{n}_{j}$ from the estimation

which can be obtained by splitting into three terms as performed in the proof of Lemma 2.1.

2.3.   Fourth-order spatial discretization of the Helmholtz equation

In this subsection, we review the BSLM ^[25,26,27] for solving (1.3). To do this, the fourth-order central finite difference method is chosen to discretize the diffusion term in space and the second-order BDF (BDF2) is used to discretize the total derivative using the departure points obtained in Subsection 2.1 with the interpolation ${\mathfrak L}$ which is reviewed in Appendix.

To this end, we first evaluate the diffusion equation (1.3) at time $t = t_{n+1}$ with the setting $s = t_{n+1}$ and then apply BDF2 to the total derivative, resulting the asymptotic one-dimensional Helmholtz equation

For the discretization of the diffusion term $u^{n}_{xx}$, we employ the finite difference operator ${\delta}_{x}^4$ with the fourth-order accuracy defined by

Using the difference operator ${\delta}_{x}^4$ and the approximations $\chi_{j}^{n-k}\; (k = 0, 1)$ for the departure points $\chi_j(t_{n-k}): = \chi(x_j, t_{n+1};t_{n-k})$ developed in the previous subsection, the Eq (2.20) can be discretized at each grid point ${x}_{j}$, as follows:

for some $\xi^k_j\; (k = 0, 1)$ between $\chi_j(t_{n-k})$ and $\chi_j^{n-k}$. Using the interpolation polynomials ${\mathfrak L} { \mathbf{U} }^{n-k}(\chi_j^{n-k})$ for ${u}^{n-k} (\chi_j^{n-k})$ ($k = 0, 1$) and truncating the asymptotic term ${\mathfrak T}^{n+1}_j$ in (2.22), which is bounded by the result of Corollary 2.3, we introduce a discrete system for the approximation $U^{n+1}_j$ of the solution $u^{n+1}_j$ given by

For convenience, let us define a vector ${\mathfrak L} { \mathbf{U} }^{n-k} (\mathit{\boldsymbol{\chi}}^{n-k})$ by

and a symmetric circulant matrix ${\mathcal A}$ for the periodic model problem by

Using these vectors and the matrix ${\mathcal A}$, the system of Eq (2.23) can be simplified as

where ${\mathcal I}$ denotes the identity matrix of size $J\times J$.

Remark 2.4. The eigenvalue computations in [20,Theorem 3.1] for ${\mathcal A}$ involving the periodic boundary condition can be used as an efficient solver for the system (2.25). Since the matrix $\mathcal{A}$ defined by (2.24) is symmetric circulant its eigenvalue decomposition is given by

where $\mathbf{q}^{j}$ is the eigenvector corresponding to the eigenvalue $\lambda_j$. With the decomposition in (2.26), the system (2.25) can be solved by the following procedure:

Remark 2.4. Note that since $\cos(2\theta)-16\cos (\theta)+15\geq 0$ for $0\leq \theta \leq 2 \pi$, all eigenvalues $\lambda_{j}$ in (2.26) are non-positive real values. Thus, the entries of the diagonal matrix $\frac{3}{2h}{\mathcal I}- \nu\Lambda$ in (2.27) are strictly positive. Therefore, the matrix $\frac{3}{2h}{\mathcal I}- \nu\Lambda$ is invertible and hence the discretized system (2.27) is uniquely solvable.

3.   Convergence and stability analysis of BSLM

This section mainly aims to present the convergence analysis of the BSLM by manipulating the bounds of the truncation errors along the discrete time evolution. To do this, we first introduce several definitions and hypotheses to be used in the subsequent analysis. Let ${\delta}_{t}$ be the first-order backward temporal finite difference operator for the vector ${\mathbf v}^{n+1}: = [v^{n+1}_1, \ldots, v^{n+1}_J]^{\mathtt T}\in {\mathbb R}^{J}$ defined by

Furthermore, let $\langle {\mathbf f}, {\mathbf g} \rangle$ be the discrete $\ell_2$-inner product defined by

and let $\|\mathbf{g}\|_{2}$ be the corresponding discrete $\ell_2$-norm for the vector ${\mathbf g}$ defined by

The convergence analysis for the BSLM presented here is based on an induction technique under mesh restriction (see ^[8,9,22]); we assume the induction hypothesis for errors at each temporal step and mesh size restriction as follows:

($A_1$) For $N\geq2$, there exists a constant $C_{A_1} > 0$ independent of $h$ and $\triangle x$ satisfying

for efficiently small $\triangle x\ll 1$.

($A_2$) We assume that the ratio between the spatial mesh size and temporal mesh size is satisfied by

for some constant $C_{A_2} > 0$.

Remark 3.1. Applying the previous hypotheses to the result of Corollary 2.3, we have

for $C_{\mathtt 1}$ depending only on $p$, $C_{A_1}$ and $C_{A_2}$, since

3.1.   Convergence analysis

We first derive an error equation for the proposed method using (2.22) and (2.23) as follows. From (2.22) and the fact that $u^m(x) - {{\mathfrak L} \mathbf{u}^m}(x) = \mathcal{O}(\triangle x^{p+1})$ for a sufficiently smooth solution $u$, subtracting (2.23) from (2.22) and multiplying the resulting equation by $\frac{2}{3}$ lead to an equation for the error ${\mathbf {e}}^{n+1}$:

where

and

Note that from Remark 3.1, each component ${r}_{j}^{n+1}$ can be estimated by

where $C_{\mathtt 2}: = \max\left\{ \frac{4}{3}C_{\mathtt 1}\|u_x^n\|_\infty, \frac{1}{3}C_{\mathtt 1}\|u^{n-1}_x\|_\infty\right\}$. Taking the discrete $\ell_2$-inner product with ${\delta}_t {{ \mathbf{e}}}^{n+1}$ after some manipulation of (3.4), one can obtain

The basic idea for estimating $\|{\mathbf {e}}^{n+1}\|_2$ is to use an induction hypothesis with a bound of $\left\| {\delta}_x {{ \mathbf{e}}}^{n+1} \right\|_{2}$ obtained by estimating each term of the Eq (3.7). We begin with an estimation of the last term $\langle\mathbf{r}^{n+1}, {\delta}_{t} {\mathbf {e}}^{n+1}\rangle$ of (3.7) in the following lemma.

Lemma 3.2. Suppose that all assumptions in Lemma 2.1 are satisfied. Then $\langle\mathbf{r}^{n+1}, {\delta}_{t} {\mathbf {e}}^{n+1}\rangle$ can be estimated by

for a constant $C_{\mathtt 3}: = {6}C_{\mathtt 2}^2 (p+1)^2$, where ${\mathcal O}_{\Gamma }$ is defined by (3.5).

Proof. Applying the triangle inequality to (3.6) shows that

By the periodicity of $e^{n-k}_j$ and the definition of discrete $\ell_2$-norm, we have

Using Young's inequality for the inner product $\langle\mathbf{r}^{n+1}, {\delta}_{t} {\mathbf {e}}^{n+1}\rangle$ completes the proof.

Next, for the estimation of $\langle {{\mathfrak L}} \mathbf{e}^{n-k}({ \mathit{\boldsymbol{\chi}}}^{n-k})- \mathbf{e}^{n-k}, \delta_t \mathbf{e}^{n+1}\rangle$ ($k = 0, 1$) in (3.7), let us define ${\widetilde \chi}_{j}(t)$ to be an arbitrary approximation of the characteristic curve ${ \chi}_{j}(t)$ over $[t_{n-1}, t_{n+1}]$ satisfying

where ${ \chi}_{j}^{n+1} = {x}_{j}$ and ${ \chi}_{j}^{n-k}$ ($k = 0, 1$) are the approximations of the departure points ${ \chi}_{j}(t_{n-k})$ and $M$ is a fixed constant. Then the property (3.8) and the fundamental theorem for a given vector $\mathbf{e}^{n-k}$ imply that

Thus, if we assume ${\widetilde \chi}_{j}(t) \in {I}_{j}$ for all $t\in [t_{n-1}, t_{n+1}]$, then by applying Hölder's inequality and using the results of Lemma A.1 on the right-hand side of (3.9), we obtain

where $d_1$ is a constant defined in Lemma A.1. Using this bound, we estimate the middle two terms of the right-hand side in (3.7) in the following lemma.

Lemma 3.3. Suppose that the approximation ${\widetilde \chi}_{j}(t)$ is included in the subinterval $I_{j}$ for $1\leq j \leq {J}$ and the derivatives $\frac{d {\widetilde \chi}_{j} }{dt}(t)$ satisfy the conditions of (3.8). Then there exist constants $C_{\mathtt 4}^k > 0, \; k = 1, 2$ depending on $p$ and $M$ for vectors ${\delta}_{x} \mathbf{e}^{n-k}: = \left[\delta_{x} e_{1}^{n-k}, \delta_{x} e_{2}^{n-k}, \dots, \delta_{x} e_{J}^{n-k} \right]^{\mathtt T}\; (k = 0, 1)$ such that

Proof. By the definition of the discrete $\ell_2$-norm and the periodicity of $\delta_{x}e^{n-k}_{j}$, the inequality (3.10) shows

Thus, applying Young's inequality and choosing $C_{\mathtt 4}^k: = \frac{1}{9}pd_1^2M^2b_k^2 (k+1)^2$ complete the proof.

We now estimate the term $\langle \mathcal A{ \mathbf{e}}^{n+1}, {\delta}_{t}{ \mathbf{e}}^{n+1} \rangle$ in (3.7). For this, let ${\delta}^2_x$ be the second-order central difference operator given by

Then the finite difference operator $\delta_{x}^4 { \mathbf{e}}^{n+1}_{j}$ defined by (2.21) can be written in terms of the operators $\delta_x$ and $\delta^2_x$ as follows

Using the expression in (3.12), we have the following result.

Lemma 3.4. For $n\geq0$, the symmetric circulant matrix ${\mathcal A}$ defined in (2.24) satisfies

Proof. From discrete summation by parts with the periodicity of $\mathbf{e}^{n+1}$, the term $\langle \mathcal A{ \mathbf{e}}^{n+1}, {\delta}_{t}{ \mathbf{e}}^{n+1} \rangle$ can be expressed as

Also, from the definition of ${\mathcal E}^{n+1}_j$ and the property ${\delta}_{x}{\delta}_t {{ \mathbf{e}}}^{n+1} = {\delta}_t {\delta}_{x}{{ \mathbf{e}}}^{n+1}$, the term $\langle {\mathcal A} {\mathbf {e}}^{n+1}, {\delta}_{t}{\mathbf {e}}^{n+1} \rangle$ in (3.13) can be split into two terms as follows:

Using the periodicity of $\mathbf{e}^{n+1}$ and discrete summation by parts again, the term ${\mathcal A_1}$ can be estimated as follows:

Similarly, from the definition of the operator $\delta_t$, the term ${\mathcal A_2}$ can be estimated as

Finally, combining (3.15) and (3.16), the term $\langle {\mathcal A} {{ \mathbf{e}}}^{n+1}, {\delta}_{t}{\mathbf {e}}^{n+1} \rangle$ can be bounded as follows:

with the help of Young's inequality. This completes the proof.

Remark 3.5. The negative of ${\mathcal A}$ is positive-definite, which allows the Choleksy decomposition to split the errors as shown in the first equation of (3.17). This is the discrete version of a common integration by parts for the Laplacian with the periodic boundary condition.

Combining the results of Lemmas 3.2–3.4, we can establish our main theorem, namely the following convergence theorem for the proposed scheme.

Theorem 3.6. Suppose that ${\delta}_t {{ \mathbf{e}}}^{1}$, ${\delta}_x{{ \mathbf{e}}}^{1}$, and ${\delta}_x^2 {{ \mathbf{e}}}^{1}$ are bounded and the hypotheses in Lemmas 3.2 and 3.3 are satisfied. Then there exists constant $C_{\mathtt{5}^*}$ such that

where $C_{\mathtt{5}^*}$ depends only on $C_{\mathtt {3}}, {C_{\mathtt 4}^k}, C_{A_1}, C_{A_2}, p, \nu, M,$ and $T$.

Proof. We apply Young's inequality on the first term of (3.7), i.e., $\left|\langle {\delta_{t}}{ \mathbf{e}}^{n}, {\delta}_{t}{{ \mathbf{e}}}^{n+1} \rangle \right| \leq \frac{1}{2}\left\|{\delta}_{t}{{ \mathbf{e}}}^{n}\right\|_{2}^2+\frac{1}{2}\left\|{\delta}_{t}{{ \mathbf{e}}}^{n+1} \right\|_{2}^2$ and combine the results of Lemmas 3.2–3.4 into (3.7) to obtain

where $C_{\mathtt 5}: = \max \left\{C_{\mathtt 3}, \frac{3}{2}C^k_{\mathtt 4}, k = 0, 1 \right\}$. By adding $\frac{1}{6}\left\|{\delta}_t {{ \mathbf{e}}}^{1} \right\|^2_{2}$ on both sides after summing from $n = 0$ to $N-1$ in (3.18) under the induction hypothesis ($A_1$) on the term $\|\delta_x \mathbf{e}^{m}\|_{2}$ for $m < N$, we have

in which the asymptotic order $h^2$ from the induction hypothesis is incorporated in ${\mathcal O}_{\Gamma}^2$. Since $\frac{1}{hT}\left\| {{ \mathbf{e}}}^{N} \right\|^2_{2} = \frac{1}{h^2N}\left\| {{ \mathbf{e}}}^{N}- {{ \mathbf{e}}}^{0} \right\|^2_{2} = \frac{1}{N}\left\|\sum\limits_{n = 1}^{N} {\delta}_t {{ \mathbf{e}}}^{n } \right\|^2_{2} \leq \sum\limits_{n = 1}^{N} \left\| {\delta}_t {{ \mathbf{e}}}^{n } \right\|^2_{2}$ from the telescoping summation, we can obtain that multiplying by $6hT$

which yields the desired estimate of $\| \mathbf{e}^{N}\|_{2}$ with the generic constant $C_{\mathtt 5^*}$ from $T$ and ${\mathcal O}_{\Gamma}$.

Remark 3.7. $(1)$ In (3.19), the estimations of $\|\delta_x \mathbf{e}^{m}\|_2$ for $m < N$ from the induction hypothesis are applied to obtain the bound of $\| \mathbf{e}^{N}\|_2$. Equations (3.19) and (3.20) show that the bound of $\|\delta_x \mathbf{e}^{N}\|_2$ increases when either $T \rightarrow \infty$ or $\nu\rightarrow 0$, and so does $\| \mathbf{e}^{N}\|_2$, cognate to a result shown in ^[1,5,22].

$(2)$ Note that for the Dirichlet problem [22,Lemma 2], it satisfies $e^N(x_j) = 0$ for all $x_j\in \partial I = \{a, b\}$, so that $\| \mathbf{e}^N\|_2\leq C\|\delta_x \mathbf{e}^{N}\|_2$ for constant $C = \frac{1}{2}(b-a)^2$, which is the discrete Poincaré inequality. But, the analysis of the model problem with the periodic boundary condition is no longer allowed to use this inequality as mentioned in the introduction. To overcome this structural drawback, we have instead investigated the error $\delta_t \mathbf{e}^n$ using the temporal finite difference operator in each time step. By applying the property of the telescoping summation of $\delta_t \mathbf{e}^n$ for $n\leq N$, the error $\mathbf{e}^{N}$ at the final step can be estimated without the use of the bound $\|\delta_x \mathbf{e}^{N}\|_2$. This approach gives more general convergence analysis of the BSLMs for solving advection-type equations without restrictions of boundary conditions compared to the strategy in ^[22].

$(3)$ The error expansion presented here is similar to the result for cases of linear advection problems ^[17] and also nonlinear problems ^[22] with Dirichlet boundary condition. In particular, the non-monotonic dependence of the temporal step size $\mathcal{O}\left({\triangle x^{p+1}}/{h}\right)$ is a well-known result, and indicates that the superior accuracy of semi-Lagrangian schemes with smaller temporal step size may not always be achieved. As expected, the other two terms in $\mathcal{O}(h^2+\triangle x^4)$ are derived from global truncation errors in the ECM strategy equipped with BDF2 for the Cauchy problem and fourth-order central finite difference for the diffusion term. Therefore, the precise convergence error of BSLM depends on the scheme of backward integration for solving the Cauchy problem as well as the interpolation formula for estimating the solution at departure points.

$(4)$ The above convergence analysis can be naturally extended in higher-dimensional problems: For example, the truncation errors of the characteristic curves on departure points in the two-dimensional problem satisfy consistent results with Lemma 2.1, Corollaries 2.2 and 2.3, (see ^[22] for details). With the help of induction hypotheses, one can show that the error equation corresponding to (3.4) satisfies the same asymptotic behavior in (3.5). Furthermore, the positive definite matrix, the negative of ${\mathcal A}$ for diffusion term (2.24) in a higher-dimensional problem can be expressed as a sum of the tensor product of the identity matrix and a positive definite matrix, so that the result of Lemma 3.4 can be easily extended. The rest of the process to obtain the error at the final step is achieved by the temporal evolution as Theorem 3.6 in a similar way. Due to space limitations, we only provide numerical results, (see Example 4.2 in Section 4).

3.2.   Stability analysis

In this subsection, we analyze the numerical stability of the proposed semi-Lagrangian scheme by showing the boundedness of the numerical solution for the model problem (1.1) with respect to the discrete ${\ell}_2$-norm without any restriction of temporal step size $h$ and spatial step size $\triangle x$. As mentioned in the Introduction, the stability analysis of the Dirichlet Problem can be obtained by using the discrete Poincaré inequality from the boundary condition, which is not permissible for the periodic problem. For the stability of the BSLM with periodic boundary, we instead begin by applying the discrete Grönwall's lemma to achieve the bound of $\delta_x {{ \mathbf{U}}}^{m}$ for each $m < N$ and then obtain the bound of ${{ \mathbf{U}}}^{N}$ by summation of $\left\| \delta_x {{ \mathbf{U}}}^{m} \right\|_{2}$ from $m = 0$ to $N-1.$

Theorem 3.8. Assume that the trajectories of the characteristic curves satisfy that $\chi^{n-k}_j\in I_j$, $(k = 0, 1)$ and the condition (3.8) for all $j = 1, \cdots, J$. If ${{ \mathbf{U}}}^{0}$ is bounded, so that ${\delta}_x {{ \mathbf{U}}}^{0}$ and ${\delta}^2_x {{ \mathbf{U}}}^{0}$ are also bounded, then the discrete solution $\mathbf{U}^N$ of the proposed BSLM is unconditionally stable in the following sense:

where $C_{s}$ depends only on $\nu$, $p$, $M$ and $T$ and is independent of $h$ and $\triangle x$.

Proof. We multiply the numerical system (2.25) by ${2}/{3}$ and add $-{ \mathbf{U}}^n/h$ to obtain that

Applying the inner product with $\delta_t{ \mathbf{U}}^{n+1}$ in (3.21) yields

We apply $\mathbf{U}^n$ instead of using $\mathbf{e}^n$ as done in Lemma 3.3 and (3.17) respectively, using similar techniques to obtain

and

where $\left\|{\delta}_x {{ \mathbf{U}}}^{m}\right\|^2_{2^*}: = \left\|{\delta}_x {{ \mathbf{U}}}^{m}\right\|^2_2+ \frac{\triangle x^2}{12}\left\|{\delta}_x^2 {{ \mathbf{U}}}^{m}\right\|^2_2$ and $C_{\mathtt 4}^k (k = 0, 1)$ are the constants in Lemma 3.3. Young's inequality applying on the first term of the right-hand side of (3.22) provides

one can rearrange Eq (3.22) using the bounds of (3.23) and (3.24), so that

where ${C_{\mathtt{4}}: = \frac{3}{2} \max \left\{C^0_{\mathtt 4}, C^1_{\mathtt 4}\right\}. }$ Summing on both sides of (3.26) from $n = 1$ to $N-1$ yields

Note that $\left\| {{ \mathbf{U}}}^{N} \right\|_{2}\leq \left\| {{ \mathbf{U}}}^{0} \right\|_{2}+h\sum\limits_{n = 1}^{N} \left\| {\delta}_t {{ \mathbf{U}}}^{n } \right\|_{2}$, so that $\left\| {{ \mathbf{U}}}^{N} \right\|^2_{2}\leq 2\left\|{{ \mathbf{U}}}^{0} \right\|^2_{2}+2Nh^2\sum_{n = 1}^{N} \left\| {\delta}_t {{ \mathbf{U}}}^{n } \right\|^2_{2}$. Multiplying by $6hT$ in (3.27) and applying the previous lower bound for $\sum_{n = 1}^{N} \left\| {\delta}_t {{ \mathbf{U}}}^{n } \right\|^2_{2}$, we obtain

Dividing (3.28) by $2\nu T$ induces that

Hence, from the discrete Grönwall's lemma we conclude that

which implies

where $C_{\mathtt{4}^*}$ satisfies

Again, applying (3.29) to (3.28), we have

We finally complete the proof by choosing $C_{s}^2 = \max\left\{ 2+\frac{12}{\nu T} C_{\mathtt{4}}T^2 e^{\frac{6C_{\mathtt{4}}T}{\nu}}, 2\nu T+{12} C_{\mathtt{4}}T^2 e^{\frac{6C_{\mathtt{4}}T}{\nu}} \right\}$.

Remark 3.9. $(1)$ The stability analysis is performed to obtain the generic constant $C_{s}$, which is independent of $\triangle x$ and $h$, using neither the induction hypothesis nor mesh restrictions. Therefore, the result says that the numerical algorithm is unconditionally stable, which is often numerically asserted. We will present numerical experiments to support this stability analysis.

$(2)$ This proof procedure offers a more general approach to achieving the stability of the BSLM than the previous stability presentation in ^[22], so that the obtained results can be applied to the BSLM in Dirichlet problems.

4.   Numerical experiments

This section provides several experiments with two simple problems to support the theoretical convergence result and the numerical stability in Theorems 3.6 and 3.8, respectively. To verify the order of convergence of the proposed scheme, we exhibit rate-varying spatial mesh sizes and temporal step sizes for varying viscosity $\nu$: the temporal convergence order, the spatial convergence order, and the reciprocal term. For the second-order temporal convergence, we perform numerical tests with different sizes of $h$ for a fixed sufficiently small size of $\triangle x$, so that the error is not much affected by the terms $\triangle x^2$ and ${\triangle x^{p+1}}/{h}$, and calculate the rate of temporal convergence (RTC) via the formula

Here, $\|\mathbf{e}^N(h)\|_2$ denotes the discrete $\ell_2$-norm error corresponding to the temporal grid size $h$ at the end time $t_N$. To investigate the order of spatial convergence from the two terms $\triangle x^4$ and $\triangle x^{p+1}/h$, we restrict the relation between two step sizes such as $h = \triangle x^2$ using the quintic Lagrange interpolation, $p = 5$, in order that

Reducing spatial and temporal step sizes by half under the restricted relation, we compute the rate of spatial convergence (RSC) via the formula

where $\mathbf{e}^N(\triangle x)$ represents the error corresponding to the spatial grid size $\triangle x$ at the end time $t_N$. Finally, to investigate the non-monotonic property of the term $\triangle x^{p+1}/h$, we apply three different orders $p = 1, 2, 3$ of Lagrangian interpolation. For each interpolation order, we demonstrate the error and its ratio by decreasing the temporal step sizes by half for a fixed spatial step size to determine the behavior of the fractional term involving temporal step size $h$ based on the analytical results already obtained.

Example 4.1. As the first experiment, we consider the Burgers' equation of form

with the initial condition given by

With the help of a Cole–Hopf transform, the periodic analytic solution ^[39] is given by

We have numerically integrated to estimate the analytic solution using Gauss–Hermite quadrature with $200$ nodes for the integral form. For the first test, we experiment with the errors and the ratio of the changes while varying temporal step sizes $h = 1/2^k, k = 6, 7, \cdots, 10$ with a fixed small spatial step size $\triangle x = 1/2^{11} \leq h$ in the different viscosity coefficients $\nu = 0.01, 0.1, 1.0$. It is easily noted that $h^2$ is the dominant term of the error, since $h^2 > \triangle x^4/h > \triangle x^4$ for $\triangle x < h \ll 1$, so the error $\| \mathbf{e}^N\|_2$ is dominated by $\mathcal{O}(h^2)$. Thus, we can check the temporal order of convergence for the method by halving temporal step sizes under the condition $\triangle x < h \ll 1$. The numerical results for the errors and their ratios are displayed for the different viscosities in Table 1. The figures for RTC in the table guarantee that the method has the required second-order temporal convergence rate.

To check the order of the spatial convergence we perform a test problem with three different viscosity coefficients $\nu = 0.01, 0.5, 1.0$, by varying step sizes under the relation $h = \triangle x^2$ such as $(h, \triangle x) = (2^{-2k}, 2^{-k})$, $k = 1, 2, \cdots, 9$. The numerical results are displayed in Table 2. The figures show that the method has an almost fourth-order spatial convergence since each term in the error expression has fourth-order convergence under the relation $h = \triangle x^2$ as seen in (4.1). This supports the theoretical order of convergence analyzed in Theorem 3.6.

Finally, in order to check the reciprocal term of the temporal step size $\mathcal{O}({\triangle x^{p+1}}/{h})$, we apply Lagrange interpolations with different orders $p = 1, 2, 3$. It can be seen that the term ${\triangle x^{p+1}}/{h}$ dominates the error as soon as $h < \triangle x^{(p+1)/3} \ll 1$, since

We display the behavior of the error for each different interpolation order $p = 1, 2, 3$ in Figure 1 while varying $h = 1/2^k, k = 3, 4, \cdots, 11$ for a fixed spatial grid size $\triangle x = 1/2^6, 1/2^{7}$ with viscosity $\nu = 0.1, 0.01$. The horizontal and vertical axes are represented by $-\ln(h)$ and $\ln(\|\mathbf{e}^N(h)\|_2)$, respectively. The behavior of the error confirms to a slop of approximately $2$ when $h \geq \triangle x^{(p+1)/3}$, since $\mathcal{O}(h^2)$ dominates the error, as seen in Figure 1. However, the error stops decreasing after the point at which the ${\triangle x^{p+1}}/{h}$ term begins to dominate it. Moreover, we can see that the error is slowly increasing to certain levels as $h$ is decreasing with $h \leq \triangle x^{(p+1)/3}$. This behavior comes from the term $\mathcal{O}(1/h)$ and thus it is not useful for the accuracy of the scheme to choose a temporal step size such as $h < \triangle x^{(p+1)/3}$.

Example 4.2. This example is the unsteady two-dimensional Burgers' equation with periodic boundary condition described as

where $\vec{u}(t, x, y) = [{u}(t, x, y), {u}(t, x, y)]^{\mathtt T}$ and the initial value is given as $u(0, x, y) = -\sin\left(\pi y \right).$

In this example, we uniformly discretize the spatial domain given by $\triangle x = \triangle y$ and present the resulting convergence orders in Theorem 3.6 as extended in the two-dimensional problem. With the lack of an exact solution for the problem, we compute reference solutions by setting $(h, \triangle x) = ({1}/{2^{8}}, {1}/{2^{8}})$ with cubic Lagrangian interpolation and $(h, \triangle x) = ({1}/{2^{17}}, {1}/{2^{8}})$ with quintic Lagrangian interpolation to estimate the temporal convergence and spatial convergence rates, respectively. To estimate the temporal order of convergence, we approach the test problem similarly to Example 1 by varying $h = {1}/{ 2^{k}}$, $k = 1, 2, \cdots, 7$ for a fixed small spatial grid size $\triangle x = {1}/{2^{11}}$, which guarantees that the two dependent spatial terms $\triangle x^4$ and $\triangle x^{p+1}/h$ become sufficiently smaller than $h^2$. The numerical results for three different viscosities $\nu = 0.01, 0.1, 1.0$ are displayed in Figure 2(a). The figures indicate that the method has an almost second-order convergence rate in time, which allows us to use a much larger temporal than spatial step size, guaranteeing the theoretical temporal order of convergence shown in Theorem 3.6 with good stability.

To check the spatial order of convergence, we reduce spatial step sizes by half based on the relation $h = \triangle x^2$, for instance $(h, \triangle x) = (1/2^{2k}, 1/2^{k})$, $k = 1, 2, \cdots, 7$ while applying the quintic Lagrange interpolation for different viscosity coefficients $\nu = 0.01, 0.1, 1.0$, as displayed in Figure 2(b). The figure shows that the proposed scheme has a fourth-order spatial convergence rate, as seen in the previous example.

In the next numerical experiment, we support the numerical stability result of Theorem 3.8 for the proposed BSLM by estimating the generic constant $C_{s}$ in the theorem as well as showing the behavior of the numerical solution over a sufficiently long time $T$ with no restriction of step sizes $\triangle x$ and $h$. The constant $C_{s}$ is computed as the ratio of the solution at between the initial time $t_0 = 0$ with its spatial differences and the final time $T = 1$ such as

We present the estimates of $C_{s}$ while varying $\triangle x = 1/2^k, \; k = 4, 5, \cdots, 10$ for fixed $h = 1/2^4$ and while varying $h = 1/2^k, \; k = 4, 5, \cdots, 10$ for fixed $\triangle x = 1/2^4$ when $\nu = 0.01, 0.05, 0.1, 0.5$, and $1$, which are plotted in Figure 3(a) and 3(b), respectively. The presented figures show that the constant $C_{s}$ is independent of the sizes of $\triangle x$ and $h$ as claimed in Theorem 3.8. We next exhibit the behavior of the discrete $\ell_2$-norm of the solution $\mathbf{U}^n$ versus time $t^n < T = 10$ with different spatial step sizes $\triangle x = 1/2^k, k = 1, 5, 9, 13$ for fixed $h = 1/2$ when $\nu = 0.01$ and $\nu = 1.0$, which are displayed in Figure 4, respectively. The experiments show that the solutions consistently decrease for a long time without blowing up under various spatial step sizes in both cases of $\nu$, even when $h$ is $2^{12}( = 4096)$ times larger than $\triangle x$ in the case $k = 13$; this provides strong evidence for the numerical stability of the scheme. Therefore, the experiment results guarantee that the numerical scheme of the proposed BSLM is stable without any restriction of the relation between mesh and temporal discretizations, as we claimed.

We end the section with a comment about the performance of the solution with respect to the viscosity. It is well known that one of the difficulties for solving Burgers' equation occurs when a small viscosity $\nu$ is applied, which generates solution curves steepen as a shock like discontinuity as time increases. This creates spurious oscillations called Gibbs phenomena near the shock front ^[18]. The design of the schemes to overcome this problem with sharp front solutions is one of the challenges ^[7].

In Figure 5, we present the behaviors of the numerical solution obtained by the proposed BSLM when viscosity coefficients $\nu = 1, 0.1, 0.01, 0.001$. The figures show the proposed scheme does not suffer from the sharp front phenomena up to the viscosity $\nu = 0.001$. For more related results or comparisons of the proposed method, we refer the interested readers to ^{[22,24,26,27,28]} and references therein.

5.   Conclusions

We have proposed a BSLM for solving Burgers' equation with the periodic boundary condition and have proved that the scheme has second-order and fourth-order convergence in time and space, respectively. In particular, the asymptotic expression possesses a reciprocal term of temporal step size that is inversely proportional to interpolation order. Furthermore, we have shown that the proposed BSLM is unconditionally stable, allowing the use of a large size of temporal step compared to spatial step size, with examples of supporting numerical experiments. The analysis of the model problem suggests an extension of linear advection–diffusion problems and high-dimensional model problems with no restriction due to boundary conditions using the same order of discretization, implying the cognate results for order of convergence. The paper contains a theoretical result, namely, the asymptotic order of convergence and numerical stability of the scheme with respect to the discrete $\ell_2$-norm with varying viscosity coefficients. The overall theoretical order of convergence and numerical stability of the scheme are also numerically verified with examples.

Acknowledgments

The first author (P.Kim) was supported by the R & D program through Korea Institute of Fusion Energy (KFE) funded by the Government funds, Republic of Korea (No. EN2241-3) and the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. 2020R1A4A1018190).

Conflict of interest

The authors declare no conflicts of interest.

Appendix

Lagrange interpolation and its properties

In this appendix, we review the Lagrange interpolation polynomial used to analyze the convergence and stability of the proposed BSLM and also introduce its properties. We consider the domain $I$ to be uniformly discretized by a spatial grid size $\triangle x$ and a fixed number $J$ as shown in (2.2).

For a given periodic function $w({x})$ defined in $I$ and each local cell $I_j: = [x_{j-1}, x_j]$, let ${\mathfrak L}_{j}w({x})$ be a $j^{th}$ local Lagrangian interpolation polynomial of degree $p\geq1$ defined by

where $w_k: = w(x_k)$ and $[m]$ denotes the integer satisfying $[m]\leq m < [m]+1$. For the interpolation polynomial near boundary points, we consider $w_{J+ k} = w_{ k}$ and $w_{-k} = w_{J- k}$ for the periodicity when $J+k$, $-k \in {\Lambda}_j$. Using the local interpolation ${\mathfrak L}_{j}w(x)$ of (A.1), a sufficiently smooth function $w(x)$ can be approximated as

and for $x\in I_j$, its truncation error (see ^[2] for detail) satisfies

for some $\xi_j \in I_j$, since $|\prod_{{j_\ell\in {\Lambda}_{j} }} (x-x_{j_\ell})| \leq\frac{1}{4} {\triangle x ^{p+1}\left[\frac{p+1}{2}\right]! \left[\frac{p+2}{2}\right]!}$ for uniformly spaced. Here, $\mathbf{1}_{{I}_{j}}$ denotes the indicator function and $\|\cdot\|_{\infty}$ is the infinity norm. Note that the interpolation ${\mathfrak L} \mathbf{w}({x})$ is piece-wise continuous in $I$ but in general not differentiable at the grid points. To overcome this non-differentiability of ${\mathfrak L} \mathbf{w}({x})$ for ${x}\in{I}_{j}$ on the interface, we define a derivative $\widetilde{D}_x {\mathfrak L} \mathbf{w}({x})$ such as

We note [27,Lemma 2.3] that for a sufficiently smooth function $w({x})$ defined on $I_j$, the derivative $\widetilde D_x$ defined in (A.4) satisfies

Lastly, we introduce some boundedness properties for the interpolation polynomial.

Lemma A.1. The derivative of ${\mathfrak L}_{j}w({x})$ defined in (A.1) can be bounded by

where ${\Lambda'_j}: = \{j_\ell\in \Lambda_j | 1\leq \ell\leq p \}, \; d_1: = d_0 p(p+1), \; d_0: = \max \limits_{x\in I_j} \Big(\prod\limits_{ \substack{k, \ell \in {\Lambda}_j\\ k\neq \ell} } \left| \frac{ x-x_k}{x_\ell-x_k}\right|\Big)$. Furthermore, for $x\in I_j$, we have

where $d_2: = 2d_0p^2(p+1)$, and $\|\cdot\|$ denotes the standard Euclidean $L_2$-norm.

Proof. First note that the Lagrange interpolation basis $L_{j, k}(x)$ of degree $p$ given by (A.1) satisfies

Hence, using the property (A.6), one can easily check that

From the definition of (A.1), we have that for any ${k \in {\Lambda}_j}$

and Eqs (A.7) and (A.8) imply

which proves the first assertion. For the second assertion, we apply the triangle inequality and the Cauchy-Schwartz inequality to obtain

which completes the proof.