Modified Bakhvalov–Shishkin Mesh Analysis

• Modified Bakhvalov–Shishkin meshes are layer-adapted meshes that merge Shishkin partitioning with Bakhvalov grading to resolve exponentially thin layers.
• They use analytic inversions of exponential decay profiles to achieve uniform step sizes and eliminate logarithmic degradation in fine-region discretizations.
• These meshes are applicable across various methods, ensuring parameter-uniform convergence for problems with boundary, interior, or overlapping layers.

Modified Bakhvalov–Shishkin meshes are layer-adapted meshes for singularly perturbed differential equations in which a Shishkin-style decomposition of the domain is combined with Bakhvalov-style grading in one or more layer subregions. Across the cited literature, the term does not denote a single universal formula. In some papers it refers to a classical Bakhvalov–Shishkin mesh used on a standard boundary-layer region; in others it denotes a genuine modification of the classical construction to accommodate interior layers, multiple perturbation parameters, overlapping layers, or Robin boundary conditions. The common purpose is parameter-uniform resolution of exponentially thin layers while reducing or removing the logarithmic degradation characteristic of piecewise-uniform Shishkin meshes (Cheng et al., 2022, Babu, 2021, Ramanujam et al., 1 Sep 2025).

1. Concept and defining features

The core idea of the modified Bakhvalov–Shishkin mesh is to preserve the explicit layer localization of a Shishkin mesh while replacing piecewise-uniform fine-region spacing by a graded map derived from the layer profile. In the literature surveyed here, this is expressed in several closely related forms.

For one-dimensional single-boundary-layer problems, a standard layer-adapted framework uses a generating function $\varphi$ on $[0,1/2]$, a transition point

$$\tau=\min\left\{\frac12,\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(\frac12\right)\right\},$$

a coarse region $[0,1-\tau]$, and a fine region $[1-\tau,1]$, with mesh nodes

$x_i= \begin{cases} \frac{2i}{N}(1-\tau), & i=0,1,\dots,N/2-1,\[0.4em] 1-\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(1-\frac{i}{N}\right), & i=N/2,\dots,N. \end{cases}$

For the classical Bakhvalov–Shishkin choice,

$$\varphi(t)=-\ln\!\big(1-2(1-N^{-1})t\big), \qquad \psi(t)=e^{-\varphi(t)}=1-2(1-N^{-1})t,$$

so that $\tau\sim (\sigma\varepsilon/\alpha)\ln N$ and $\max|\psi'|=2$ (Yan et al., 2022, Cheng et al., 2022).

This construction already exhibits the structural hallmark of the MBS family: the transition is Shishkin-like, but the layer region is graded by inversion of an exponential-type profile rather than meshed uniformly. A closely related description appears in the finite element analysis of Bakhvalov-type meshes, where an explicit transition is combined with a logarithmic grading branch and an essentially uniform coarse region; that framework is presented as encompassing MBS-type choices when the breakpoint is selected to reproduce a Shishkin-style transition in physical space (Zhang et al., 2020).

A further strand of the literature uses the term “Shishkin–Bakhvalov mesh” for the same principle: transition points are chosen as in a Shishkin mesh, whereas the layer subintervals are graded by inverting the relevant exponential layer functions. This terminology appears for two-parameter elliptic and parabolic problems with discontinuous data (Roy et al., 2022, Roy et al., 2022).

2. Terminological variability and scope

The phrase “Modified Bakhvalov–Shishkin mesh” is not used uniformly. That variability is itself a substantive feature of the literature.

In the LDG analysis of a two-dimensional convection–diffusion problem, the authors explicitly identify the Bakhvalov–Shishkin-type mesh used in the paper as the “modified” Bakhvalov–Shishkin mesh. There the 1D generator is

$$\phi(t)=-\ln\big(1-2(1-N^{-1})t\big), \qquad t\in[0,1/2],$$

with $[0,1/2]$0, and the 2D grid is the tensor product of the corresponding $[0,1/2]$1- and $[0,1/2]$2-meshes (Cheng et al., 2022).

By contrast, the third-order LDG paper on convection–diffusion type problems explicitly states that no modified BS mesh is introduced or analyzed. It treats only the Shishkin-type mesh, the classical Bakhvalov–Shishkin mesh, and a Bakhvalov-type mesh. Its results for BS therefore refer to the classical choice

$[0,1/2]$3

not to a separate modified variant (Yan et al., 2022).

For problems with interior layers or multiple small parameters, “modified” denotes a genuine extension of the classical BS mesh. In the SDFEM study of a weakly coupled $[0,1/2]$4 convection–diffusion system with a discontinuous source at $[0,1/2]$5, the modification consists of introducing two graded subregions, one on $[0,1/2]$6 for the weak interior layer and one on $[0,1/2]$7 for the boundary layer, with separate generating functions $[0,1/2]$8 and $[0,1/2]$9 and transition points at $$\tau=\min\left\{\frac12,\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(\frac12\right)\right\},$$0, $$\tau=\min\left\{\frac12,\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(\frac12\right)\right\},$$1, and $$\tau=\min\left\{\frac12,\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(\frac12\right)\right\},$$2 (Babu, 2021).

For coupled reaction–diffusion systems with Robin boundary conditions and two perturbation parameters $$\tau=\min\left\{\frac12,\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(\frac12\right)\right\},$$3, the modified BS mesh retains a five-subinterval Shishkin partition but replaces uniform layer spacing by explicit exponential-inversion formulas in four layer subintervals. The paper stresses that this yields the global bound $$\tau=\min\left\{\frac12,\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(\frac12\right)\right\},$$4 and exact second-order convergence in the maximum norm (Ramanujam et al., 1 Sep 2025).

This suggests that “MBS mesh” should be understood as a family of constructions rather than a single canonical mesh.

3. Canonical construction patterns

Despite the terminological differences, the papers exhibit a small number of recurrent design patterns.

Boundary-layer-only modification

In the simplest setting, modification means smooth Bakhvalov grading inside a single fine region attached to one boundary. The BS-type generator

$$\tau=\min\left\{\frac12,\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(\frac12\right)\right\},$$5

gives $$\tau=\min\left\{\frac12,\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(\frac12\right)\right\},$$6, $$\tau=\min\left\{\frac12,\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(\frac12\right)\right\},$$7, and $$\tau=\min\left\{\frac12,\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(\frac12\right)\right\},$$8, in contrast to the Shishkin choice $$\tau=\min\left\{\frac12,\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(\frac12\right)\right\},$$9 with $[0,1-\tau]$0 (Cheng et al., 2022, Yan et al., 2022). In this form, the modification is analytic rather than topological: the domain split is unchanged, but the fine-region distribution is smoother and more strongly graded.

Interior-layer and boundary-layer modification

When a discontinuity induces a weak interior layer, the mesh requires an additional graded segment. In the coupled convection–diffusion system with a discontinuous source at $[0,1-\tau]$1, the modified BS mesh uses four subregions with $[0,1-\tau]$2 points each:

1. $[0,1-\tau]$3, equidistant;
2. $[0,1-\tau]$4, graded via $[0,1-\tau]$5;
3. $[0,1-\tau]$6, equidistant;
4. $[0,1-\tau]$7, graded via $[0,1-\tau]$8.

The transition widths are

$[0,1-\tau]$9

and the modified BS functions are

$[1-\tau,1]$0

The interior-layer grading on $[1-\tau,1]$1 is precisely the feature absent from the classical single-boundary BS mesh (Babu, 2021).

A related six-subinterval design appears for two-parameter one-dimensional problems with discontinuous convection and source terms. There the interval is partitioned as

$[1-\tau,1]$2

with transition points

$[1-\tau,1]$3

and Bakhvalov-type inversion used in each layer subinterval (Roy et al., 2022). The parabolic analogue uses the same spatial partition in combination with Crank–Nicolson time stepping (Roy et al., 2022).

Overlapping-layer modification for two small parameters

For weakly coupled reaction–diffusion systems with $[1-\tau,1]$4, the modified BS mesh uses five subintervals:

• $[1-\tau,1]$5,
• $[1-\tau,1]$6,
• $[1-\tau,1]$7,
• $[1-\tau,1]$8,
• $[1-\tau,1]$9,

with

$x_i= \begin{cases} \frac{2i}{N}(1-\tau), & i=0,1,\dots,N/2-1,\[0.4em] 1-\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(1-\frac{i}{N}\right), & i=N/2,\dots,N. \end{cases}$0

and, in the regime $x_i= \begin{cases} \frac{2i}{N}(1-\tau), & i=0,1,\dots,N/2-1,\[0.4em] 1-\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(1-\frac{i}{N}\right), & i=N/2,\dots,N. \end{cases}$1,

$x_i= \begin{cases} \frac{2i}{N}(1-\tau), & i=0,1,\dots,N/2-1,\[0.4em] 1-\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(1-\frac{i}{N}\right), & i=N/2,\dots,N. \end{cases}$2

The node formulas are explicit exponential inversions tailored separately to the $x_i= \begin{cases} \frac{2i}{N}(1-\tau), & i=0,1,\dots,N/2-1,\[0.4em] 1-\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(1-\frac{i}{N}\right), & i=N/2,\dots,N. \end{cases}$3- and $x_i= \begin{cases} \frac{2i}{N}(1-\tau), & i=0,1,\dots,N/2-1,\[0.4em] 1-\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(1-\frac{i}{N}\right), & i=N/2,\dots,N. \end{cases}$4-layers on the left and right, with a uniform middle region (Ramanujam et al., 1 Sep 2025).

4. Mesh-generating functions and analytic properties

A recurrent analytic objective of MBS constructions is the replacement of Shishkin’s $x_i= \begin{cases} \frac{2i}{N}(1-\tau), & i=0,1,\dots,N/2-1,\[0.4em] 1-\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(1-\frac{i}{N}\right), & i=N/2,\dots,N. \end{cases}$5 fine-region steps by $x_i= \begin{cases} \frac{2i}{N}(1-\tau), & i=0,1,\dots,N/2-1,\[0.4em] 1-\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(1-\frac{i}{N}\right), & i=N/2,\dots,N. \end{cases}$6 or otherwise sharpened bounds.

In the reaction–diffusion system with Robin boundary conditions, the modified BS mesh satisfies

$x_i= \begin{cases} \frac{2i}{N}(1-\tau), & i=0,1,\dots,N/2-1,\[0.4em] 1-\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(1-\frac{i}{N}\right), & i=N/2,\dots,N. \end{cases}$7

which is singled out as the key mesh-size lemma underlying the exact second-order maximum-norm bound (Ramanujam et al., 1 Sep 2025). By contrast, the standard Shishkin mesh on the same problem has step sizes $x_i= \begin{cases} \frac{2i}{N}(1-\tau), & i=0,1,\dots,N/2-1,\[0.4em] 1-\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(1-\frac{i}{N}\right), & i=N/2,\dots,N. \end{cases}$8 in the layer subintervals.

In the SDFEM study with an interior layer at $x_i= \begin{cases} \frac{2i}{N}(1-\tau), & i=0,1,\dots,N/2-1,\[0.4em] 1-\frac{\sigma\varepsilon}{\alpha}\varphi\!\left(1-\frac{i}{N}\right), & i=N/2,\dots,N. \end{cases}$9, the modified BS mesh has bounded $$\varphi(t)=-\ln\!\big(1-2(1-N^{-1})t\big), \qquad \psi(t)=e^{-\varphi(t)}=1-2(1-N^{-1})t,$$0, whereas the Shishkin comparator satisfies $$\varphi(t)=-\ln\!\big(1-2(1-N^{-1})t\big), \qquad \psi(t)=e^{-\varphi(t)}=1-2(1-N^{-1})t,$$1. This boundedness is directly linked to the improvement from $$\varphi(t)=-\ln\!\big(1-2(1-N^{-1})t\big), \qquad \psi(t)=e^{-\varphi(t)}=1-2(1-N^{-1})t,$$2 to $$\varphi(t)=-\ln\!\big(1-2(1-N^{-1})t\big), \qquad \psi(t)=e^{-\varphi(t)}=1-2(1-N^{-1})t,$$3 in the discrete energy norm (Babu, 2021).

For the two-dimensional LDG supercloseness analysis, the BS-type mesh yields $$\varphi(t)=-\ln\!\big(1-2(1-N^{-1})t\big), \qquad \psi(t)=e^{-\varphi(t)}=1-2(1-N^{-1})t,$$4 and $$\varphi(t)=-\ln\!\big(1-2(1-N^{-1})t\big), \qquad \psi(t)=e^{-\varphi(t)}=1-2(1-N^{-1})t,$$5, while the Shishkin mesh has $$\varphi(t)=-\ln\!\big(1-2(1-N^{-1})t\big), \qquad \psi(t)=e^{-\varphi(t)}=1-2(1-N^{-1})t,$$6. The fine-region sizes satisfy $$\varphi(t)=-\ln\!\big(1-2(1-N^{-1})t\big), \qquad \psi(t)=e^{-\varphi(t)}=1-2(1-N^{-1})t,$$7, and the transition thickness is $$\varphi(t)=-\ln\!\big(1-2(1-N^{-1})t\big), \qquad \psi(t)=e^{-\varphi(t)}=1-2(1-N^{-1})t,$$8 under the singularly perturbed assumption (Cheng et al., 2022).

In the analysis of Bakhvalov-type finite element meshes for two-parameter elliptic problems, the logarithmic generating functions transform exponential layer decay into algebraic behavior in the mesh parameter. For example, on a left exponential layer subinterval,

$$\varphi(t)=-\ln\!\big(1-2(1-N^{-1})t\big), \qquad \psi(t)=e^{-\varphi(t)}=1-2(1-N^{-1})t,$$9

and one obtains bounds such as

$\tau\sim (\sigma\varepsilon/\alpha)\ln N$0

This analytic mechanism is presented as characteristic of the MBS philosophy: Shishkin-style subdomain decomposition combined with Bakhvalov-style grading inside each fine subdomain (Zhang et al., 2020).

A different modified Bakhvalov construction for semilinear reaction–diffusion problems emphasizes smoothness of the mesh generator itself. There the mesh-generating function $\tau\sim (\sigma\varepsilon/\alpha)\ln N$1 is piecewise defined, includes a cubic transition segment, is symmetric with respect to $\tau\sim (\sigma\varepsilon/\alpha)\ln N$2, and satisfies

$\tau\sim (\sigma\varepsilon/\alpha)\ln N$3

These second-order mesh-regularity bounds are central to the proof of $\tau\sim (\sigma\varepsilon/\alpha)\ln N$4-uniform $\tau\sim (\sigma\varepsilon/\alpha)\ln N$5 convergence for the fitted difference schemes considered there (Karasuljić et al., 2018).

5. Discretization frameworks and convergence results

Modified Bakhvalov–Shishkin meshes are not tied to a single discretization. The literature associates them with finite elements, SDFEM, LDG, and fitted finite differences, and the mesh often determines whether logarithmic factors remain in the final estimate.

For a weakly coupled convection–diffusion system with discontinuous source term, continuous piecewise linear finite elements with streamline diffusion on the modified BS mesh satisfy

$\tau\sim (\sigma\varepsilon/\alpha)\ln N$6

whereas the corresponding Shishkin-mesh bound is

$\tau\sim (\sigma\varepsilon/\alpha)\ln N$7

The constants are independent of $\tau\sim (\sigma\varepsilon/\alpha)\ln N$8 under $\tau\sim (\sigma\varepsilon/\alpha)\ln N$9 (Babu, 2021).

For the coupled reaction–diffusion system with Robin boundary conditions, the interior central difference scheme together with cubic spline approximations of the boundary derivatives yields

$\max|\psi'|=2$0

under the M-matrix conditions and $\max|\psi'|=2$1 (Ramanujam et al., 1 Sep 2025).

For the two-dimensional LDG method on a BS-type mesh identified there as MBS, the main supercloseness result is

$\max|\psi'|=2$2

and the associated $\max|\psi'|=2$3 error estimate is

$\max|\psi'|=2$4

uniformly in $\max|\psi'|=2$5, assuming $\max|\psi'|=2$6 and $\max|\psi'|=2$7 (Cheng et al., 2022).

For the third-order singularly perturbed convection–diffusion-type problem, the LDG analysis on the classical BS mesh—not an MBS mesh—gives

$\max|\psi'|=2$8

again uniform in $\max|\psi'|=2$9 in the convection-dominated regime $$\phi(t)=-\ln\big(1-2(1-N^{-1})t\big), \qquad t\in[0,1/2],$$0 (Yan et al., 2022). This is important because it shows that some benefits often attributed to “modified” BS grading are already present for the classical BS mesh.

For two-parameter one-dimensional problems with discontinuous data, upwind finite difference schemes on Shishkin–Bakhvalov meshes achieve parameter-uniform first-order convergence:

$$\phi(t)=-\ln\big(1-2(1-N^{-1})t\big), \qquad t\in[0,1/2],$$1

for the elliptic case (Roy et al., 2022), and

$$\phi(t)=-\ln\big(1-2(1-N^{-1})t\big), \qquad t\in[0,1/2],$$2

for the parabolic case with Crank–Nicolson in time (Roy et al., 2022).

For semilinear reaction–diffusion problems, fitted schemes on a modified Bakhvalov mesh satisfy the $$\phi(t)=-\ln\big(1-2(1-N^{-1})t\big), \qquad t\in[0,1/2],$$3-uniform estimate

$$\phi(t)=-\ln\big(1-2(1-N^{-1})t\big), \qquad t\in[0,1/2],$$4

with numerical evidence supporting the theoretical rate (Karasuljić et al., 2018). A different semilinear study proves $$\phi(t)=-\ln\big(1-2(1-N^{-1})t\big), \qquad t\in[0,1/2],$$5 convergence on a modified Shishkin mesh of order $$\phi(t)=-\ln\big(1-2(1-N^{-1})t\big), \qquad t\in[0,1/2],$$6, while using a modified Bakhvalov mesh only in numerical experiments, where observed rates are near $$\phi(t)=-\ln\big(1-2(1-N^{-1})t\big), \qquad t\in[0,1/2],$$7 and absolute errors are smaller than on the Shishkin meshes considered there (Karasuljić et al., 2020).

6. Comparison with Shishkin, BS, and Bakhvalov meshes

The role of the modified Bakhvalov–Shishkin mesh is best understood relative to neighboring mesh families.

Mesh type	Typical transition philosophy	Fine-region distribution
Shishkin	Explicit layer width, usually $$\phi(t)=-\ln\big(1-2(1-N^{-1})t\big), \qquad t\in[0,1/2],$$8	Piecewise uniform
Classical BS	Shishkin-style transition with BS generator	Bakhvalov–Shishkin grading
Modified BS / Shishkin–Bakhvalov	Shishkin-style transition, often extended to multiple layer regions	Bakhvalov-type grading in each layer subregion
Bakhvalov-type	Layer-width control tied directly to perturbation parameter	Fully parameter-driven grading

The Shishkin mesh remains a common baseline because of its simplicity. In the third-order LDG study, it is described as simple to construct and reliable, but its energy-norm error includes logarithmic factors absent or reduced on BS/B meshes (Yan et al., 2022). In the reaction–diffusion system with Robin data, the Shishkin mesh produces “almost” second-order accuracy but incurs a $$\phi(t)=-\ln\big(1-2(1-N^{-1})t\big), \qquad t\in[0,1/2],$$9 degradation in the maximum norm (Ramanujam et al., 1 Sep 2025).

The classical BS mesh removes much of this degradation while retaining explicit transition points. In the LDG supercloseness study, the BS-type mesh is exactly the mesh termed MBS, and it yields $[0,1/2]$00 rather than $[0,1/2]$01 in the supercloseness factor $[0,1/2]$02 (Cheng et al., 2022). In the third-order LDG study, the classical BS mesh is presented as stronger than the Shishkin mesh in the fine region because its local stretching can exceed the Shishkin grading by a factor up to $[0,1/2]$03 (Yan et al., 2022).

Modified BS constructions become necessary when the layer geometry is more complicated than a single boundary layer. The additional graded segment near an interior discontinuity in (Babu, 2021), the four graded layer subintervals for overlapping $[0,1/2]$04- and $[0,1/2]$05-layers in (Ramanujam et al., 1 Sep 2025), and the six-subinterval Shishkin–Bakhvalov layouts in (Roy et al., 2022) and (Roy et al., 2022) are all examples of this broader use.

A plausible implication is that “modification” is best interpreted structurally: the mesh departs from the classical single-layer BS template whenever the solution decomposition contains multiple distinct layer components that require separate graded subdomains.

7. Misconceptions, limitations, and research context

A common misconception is that every paper mentioning BS-type meshes uses a distinct modified BS definition. The literature here shows otherwise. One paper explicitly states that no modified BS mesh is used (Yan et al., 2022), another uses “modified BS” to mean the ordinary BS-type mesh in its LDG setting (Cheng et al., 2022), while several others reserve the term for constructions with extra graded subregions or multi-parameter overlap handling (Babu, 2021, Ramanujam et al., 1 Sep 2025).

A second misconception is that the sole benefit of MBS meshes is smaller constants. In several works the improvement is qualitative at the level of logarithmic factors or even of the convergence order visible in the stated estimate. Representative examples are the change from $[0,1/2]$06 to $[0,1/2]$07 in the SDFEM discrete energy norm (Babu, 2021), and from $[0,1/2]$08 to $[0,1/2]$09 in the maximum norm for the Robin problem (Ramanujam et al., 1 Sep 2025).

At the same time, the literature does not support a universal claim that MBS is always superior. In the two-dimensional LDG experiments, BS and B meshes are reported as comparable, with BS sometimes giving slightly smaller errors for the tested $[0,1/2]$10 while B shows slightly better asymptotic rates (Cheng et al., 2022). In the semilinear reaction–diffusion experiments comparing Shishkin, modified Shishkin, modified Bakhvalov, and Liseikin meshes, the Liseikin mesh often produces the smallest errors, although the modified Bakhvalov mesh substantially improves on the Shishkin meshes (Karasuljić et al., 2020).

The broader research context is therefore not a single canonical mesh formula, but a design principle: identify the asymptotic layer structure, introduce transition points that isolate each layer component, and replace piecewise-uniform fine meshes by Bakhvalov-type inversions of the relevant exponential profiles. Under that principle, the modified Bakhvalov–Shishkin mesh serves as a flexible intermediary between Shishkin simplicity and Bakhvalov grading efficiency across finite element, discontinuous Galerkin, streamline diffusion, and fitted finite difference discretizations (Zhang et al., 2020, Zhang et al., 2020).