A second order numerical method for a class of parameterized singular perturbation problems on adaptive grid

1 Introduction

The study of singularly perturbed problems (SPPs) is exceptionally useful because they describe many problems of practical interest. These types of problems are identified by the differential equations in which a small parameter (ε) known as perturbation parameter multiplies the highest order derivative. It is well known that the classical numerical methods for solving SPPs are unstable and fail to give accurate results when ‘ε’ is very small, 0 < ε ≪ 1. Therefore, special approaches are required for obtaining a good approximation. So, it is important to develop such methods, whose accuracy do not depend on ε. The layer-adapted meshes: Shishkin type meshes (standard Shishkin mesh (S-mesh), Bakhvalov-Shishkin mesh (B-S-mesh)) and newly developed adaptive grid can be used to obtain parameter uniform numerical methods. In recent years, development of parameter uniform numerical methods based on layer adapted meshes have witnessed substantial progress. For more details on the singular perturbation and layer adapted meshes, one can refer the books ([7, 13, 19]) and the references therein.

Consider the following singularly perturbed Boundary Value Problem(BVP) depending on a parameter:

where 0 < ε ≪ 1 and λ known as the control parameter. Here s₀, s₁ are the given constants. The function f(x,u,λ) is assumed to be sufficiently smooth and satisfies the following bounds:

With these assumptions, the BVP (1.1) possesses a unique solution having a boundary layer of width O(ε) near x = 0 (refer [1, 17, 18]). The parameter λ has no connection with the eigenvalue of the nonlinear differential equation. Since there are two unknowns, two boundary conditions are given in (1.1) to determine it exactly.

Parameterized BVP have been considered for many years. The existence and uniqueness of the solution for the BVP(1.1) were first considered by Pomentale [17]. Jankowski [8] and Liu [11] constructed monotonic iterative methods for these problems. But the above-mentioned articles were concerned with the regular case. In recent years, many researchers considered the singular perturbation cases for these problems. Amiraliyev et. al. [1] gave a uniform finite difference method on a standard Shishkin mesh [7] for BVP(1.1) and shown that the method is first order convergent up to a logarithmic factor i.e. of the order O(N^{− 1} ln N). A hybrid difference scheme which combines upwind scheme on the fine mesh with the midpoint upwind scheme on the coarse mesh was considered by Cen [3]. Xie et. al. [23] used the boundary layer correction technique to solve the parameterized problem. Turkyilmazoglu [21] constructed a methodology based on the homotopy analysis technique to approximate the analytic solution. A uniform difference scheme was developed by Amiraliyeva et. al. [2] for parameterized delay differential equation. Recently, Das [6] provide a comparison of apriori and posteriori meshes for BVP (1.1).

Here, we have considered two types of numerical schemes: the backward Euler scheme and the post-processing Richardson extrapolation technique on the adaptive grid generated via equidistribution principle. Richardson extrapolation is a technique where two computed solutions are approximated by an average to provide a better approximation. This technique has the advantage that it can be extended to problems in more than one dimensions. Vulanovic et. al. [22] used this technique for singularly perturbed reaction-diffusion problem using Bakhvalov mesh. Its application to the numerical solution of singularly perturbed convection-diffusion problem was examined in [16] by Natividad and Stynes. Using this idea, Mohapatra and Natesan [14] constructed a second order post processing technique for solving singularly perturbed delay differential equation. In this work, we have solved the BVP (1.1) by backward Euler difference scheme on an adaptive grid. Then, we have successfully applied the Richardson extrapolation technique on computed solution to enhance the accuracy from first order to second order.

This paper is organized as follows: In section 3, we describe the numerical scheme for the problem (1.1). In section 4, a brief description of the nonuniform mesh is presented. The error estimates for the approximate solution are obtained for the proposed schemes in section 5. Finally, in section 6, we present two numerical examples, which validate the theoretical estimates. Throughout this paper ‘C’ denotes a generic positive constant independent of both ε and N which can take different values at different places and subscripted C’s are fixed constants. Here, we denote g(x_i) = g_i.

2 Analytical results

3 Numerical schemes

In this section, we describe a difference scheme for BVP (1.1). On a arbitrary nonuniform mesh Ω^N:0 < x₀ < x₁ < … < x_N, define the operator:

where h_i = x_i − x_{i − 1}. The BVP (1.1) is discretized by the following difference scheme:

In order to increase the accuracy of the difference scheme (3.1), we follow the idea described in [16]. We construct a new nonuniform mesh Ω2N:0<x~0<x~1<…<x~2N and h~i=x~i−1−x~i, which is generated by bisecting each interval of the original mesh Ω^N. Now, the difference scheme on the mesh Ω^2N is defined as:

We are interested in the extrapolated solution defined on the mesh Ω^N by Uextp,iN=2U¯2i2N−UiN, for i=0,1,…,N; which we expect will improve the accuracy of the approximation obtained at the points x_i ∈ Ω^N. Also, It will enhance the order of the convergence of the scheme.

4 The nonuniform mesh

Numerical methods using the standard finite difference schemes on uniform meshes are inadequate for solving SPPs [7, 13]. Special approaches are required to dealt with these problems. Layer-adapted meshes are of great interest, not only because they capture the layer behaviour, but also because they can result in numerical approximation whose accuracy is guaranteed independent of the width of the layer. These meshes can be divided into two categories: a priori mesh, for which a prior information about the location and the width of the solution is required and a posteriori mesh, for which we do not need such information. In past few years, different approaches are developed to get uniform convergence where the mesh is chosen apriori. Very few literature are available for the posteriori meshes. One of such approach is the adaptive grid. In [10], an adaptive grid based numerical method is developed for solving a quasi-linear one-dimensional convection-diffusion problem. Mackenzie [12] constructed the uniform convergent upwind method for a convection-diffusion problem on an adaptive grid generated via equidistribution principle. Chen [4] provided uniform convergence of finite difference approximation for SPP on an adaptive grid. Mohapatra and Natesan [15] proved the uniform convergence of upwind scheme for approximating the global solution and the global normalized flux using the grid equidistribution. Here, our aim is to construct efficient numerical method on an adaptively generated posteriori mesh.

Adaptive grid is one of the special kind of nonuniform mesh. A commonly-used technique to generate the adaptive grid is based on equidistribution of an arbitrary non-negative function M(u(x),x) defined on [0, 1] known as the monitor function. A grid Ω^N is said to be equidistributed if

The monitor functions are usually depend on the gradient of the solution. The solution of (4.1) along with a discretized version of the BVP (1.1) produces the numerical approximation to the solution of the BVP (1.1).

We shall follow the similar mesh generation algorithm for solving the parameterized BVP (1.1). The following adaptive algorithm is used to generate the appropriate nonuniform mesh (One may refer [4, 10, 12] for more details on the adaptive algorithm).

5 Error estimates

In this section, we wish to establish the uniform convergence of the extrapolation scheme on the layer-adapted mesh. To investigate the error function of the method, the error function is denoted as ziN=UiN−ui,0≤i≤N,μ^N = λ^N − λ which is the solution of the following discrete problem:

For, i = 1,2,…,N, we use the Taylor’s expansion for f around (x_i,u_i,λ) which gives,

where,

It is easy to verify that the matrix associated with L^N is an M-matrix. Hence, we can use the discrete maximum principle, which states “If {vi}0N and {wi}0N are the mesh functions that satisfy v₀ ≥ w₀ and v_N ≥ w_N and L^Nv_i ≥ L^Nw_i for i = 0,…,N − 1 then v_i ≥ w_i for all i.” Here we use the discrete norms ∥v∥∞,Ω=maxi|vi| and ∥v∥∗,Ω=maxi=1,…,N|∑j=iNhjvj|.

6 Numerical Results

We consider two test problems to show the applicability and efficiency of the proposed methods.

Let u(x) be the exact solution. Let the numerical solution before and after extrapolation be UiN and Uextp,iN respectively, then the Eε,uN with respect to the variable u of the numerical solution before and after extrapolation is calculated as: maxi|ui−UiN| and maxi|ui−Uextp,iN| respectively. The corresponding rate of convergence is calculated as rε,uN=log2Eε,uNEε,u2N.

Figure 1(a) and Figure 1(b) represent the numerical solution before and after extrapolation for Example 6.1 and Example 6.2 respectively for N = 40 and ε = 2⁻⁸. The error behaviour of the scheme for Example 6.1 and 6.2 for N = 20 and ε = 10⁻² are displayed in the Figure 2. One can clearly observe that the error in the layer region as well as in the outer region is less after extrapolation. In Figure 3(a) and Figure 3(b), the maximum pointwise error versus the number of mesh interval is plotted in the logarithmic scale for ε = 10⁻⁸. These figures show the effectiveness and accuracy of the extrapolation scheme. In Table 1, we represent the maximum pointwise error and the corresponding rate of convergence for Example 6.1 with ε = 10⁻⁴, 10⁻⁸ which is enough to show the singularly perturbed nature of the problem. The rapid decrease in the error after extrapolation can be observed from these results. Moreover, improvement in the rate of convergence from first order to second order is also evident as claimed by the theoretical finding. The proposed improvement is compared with results on Shishkin type meshes in Table 2 for Example 6.2 with ε=10⁻⁸. The numerical approximation on the B-S-mesh also results in first-order convergence before extrapolation and second order convergence after extrapolation. But the advantage of the adaptive grid is that it does not need any apriori information about the location and width of the layer. The numerical results are the clear illustration of the error estimates.

Figure 1

Numerical solution with N = 40 and ε = 2⁻⁸.

Figure 2

Error behaviour with N = 20 and ε = 10⁻².

Figure 3

Loglog plot of maximum point-wise errors with ε = 10⁻⁸.

7 Concluding remarks

In this work, a higher order accurate method is used to compute the solution of a class of parameterized nonlinear singularly perturbed differential equation. First, we use the backward Euler difference scheme on the adaptive grid, then Richardson extrapolation technique is used on the computed solution. Clearly, after extrapolation the error of the numerical solution is considerably decreased. Theoretical estimates are supported with the help of numerical results. The advantage of the adaptive grid over Shishkin type meshes is evident from the theoretical estimates obtained as well as from the numerical results shown.