An efficient adaptive grid method for a system of singularly perturbed convection-diffusion problems with Robin boundary conditions

A system of singularly perturbed convection-diffusion equations with Robin boundary conditions is considered on the interval [0,1]\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$[0,1]$\end{document}. It is shown that any solution of such a problem can be expressed to a system of first-order singularly perturbed initial value problem, which is discretized by the backward Euler formula on an arbitrary nonuniform mesh. An a posteriori error estimation in maximum norm is derived to design an adaptive grid generation algorithm. Besides, in order to establish the initial values of the original problems, we construct a nonlinear optimization problem, which is solved by the Nelder–Mead simplex method. Numerical results are given to demonstrate the performance of the presented method.


2)
© The Author(s) 2020. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. (1.4) where α = min{α 1 , α 2 }, k = 1, 2. Under these conditions, there exists a unique solution − → u (t) of problem (1.1). Such type of problems appear in optimal control problems and in certain resistor-capacitor circuits (see [1]).
It is well known that many numerical methods have been developed for the solution of singularly perturbed problems (see [2][3][4] for example). To capture the singularly perturbed nature of these problems, one effective approach is based on the use of meshes that are designed to be very fine where sharp layer appears in the solution. Such meshes contain two classes: special meshes (e.g., Shishkin meshes) that are chosen by a priori information and adaptive meshes that are generated by a grid iterative algorithm. Over the past few decades, the special mesh approach [5][6][7][8][9][10] was used to solve some system of singularly perturbed problems. The authors in [11,12] proposed adaptive grid methods to solve the system of singularly perturbed convection-diffusion problems with Dirichlet boundary conditions. As far as we know, the adaptive grid methods for the system of singularly perturbed convection-diffusion Robin boundary problems are not found in the literature.
In this paper, we propose an adaptive grid method to solve the above problem (1.1).
An upwind finite difference scheme is developed to approximate the system of first-order singularly perturbed differential equations transformed from problem (1.1). Then, an a posteriori error estimate for the presented finite difference scheme is derived, which is used to design an adaptive grid algorithm(see [13]). In addition, in order to establish the approximation values of u 1 (0) and u 2 (0), we design a nonlinear optimization problem, which is solved by the Nelder-Mead simplex method [14]. Finally, linear and nonlinear numerical examples are used to verify the effectiveness and practicability of the proposed adaptive grid method.
Notations. Throughout this paper, let C be a generic positive constant that is independent of all perturbation parameters ε i (i = 1, 2) and mesh parameter N . It may take different values in different place. Besides, in our estimates, we use the L ∞ norm defined by v(x) ∞ = ess sup

Reformulation of the boundary value problem
Our numerical method for solving (1.1) is based on reformulating it as a system of firstorder singularly perturbed differential equations, for which we will design an adaptive grid method. Now, we integrate the two sides of the first two equations of problem (1.1) from 0 to x and write it into the following matrix form: Then, from the Robin boundary conditions of problem (1.1), we have It should be pointed out that the constant vector − → η is an initial condition which will be adjusted so that the solution satisfies the right Robin boundary condition of problem (1.1). In recent years, the numerical methods of this type of equations were discussed in [15,16]. Our numerical method for solving (1.1) is based on discretizing (2.1). That is, we solve our original Robin boundary problems by applying the adaptive grid approach to solve problem (2.1). This numerical method is cheap: a direct method to settle (1.1) would entail solving a linear system of equations with 2N + 2 unknowns, but when using (2.1) instead, we only solve N second-order linear systems.
For the later analysis of our numerical method, further information about the structure of the continuous solution − → u (x) of (2.1) is needed. We shall mainly consider the following stability result for problem (2.1).

Lemma 2.1 (Stability result) The solution
− → u (x) of the initial value problem (2.1) satisfies the following inequality: Proof The proof can be seen in Theorem 2.4 of [15].

Finite difference discretization
. . , N , be the local mesh size. Then the finite difference scheme of (2.1) can be constructed as follows: dt. Next, we give the stability property for the discrete scheme (3.1). The discrete maximum principle implies the following discrete stability result.
Proof The proof can be seen in Lemma 3.2 of [14].

Nonuniform mesh generation algorithm
It is well known that many researchers developed an adaptive grid method to solve the single singularly perturbed differential equation (see, e.g., [15][16][17][18][19]). In these works, the authors used the arc-length monitor function to design a grid generation algorithm. Recently, Liu and Chen [13] proposed an adaptive grid method for a system of singularly perturbed differential equations. In this paper, similar to [13,17], we also choose the following monitor function: where U j (x) ∈ C(0, 1] is a piecewise linear interpolation function through the knots (x i , U j,i ), j = 1, 2, i = 0, 1, . . . , N . Thus, the key technologies of the adaptive mesh method is to find ( where j = 1, 2, i = 0, 1, . . . , N -1. Here, in order to solve the above equidistribution problem (3.2), we construct the following grid iteration algorithm [17]. Step . . , N} as an initial mesh.
Step 2. For given k, assume that the mesh¯ (k) be the maximum arc-length between the points ( Then Step 3. Choose a user-chosen constant C 0 > 1. If C 0 satisfies the following stopping criterion: then go to Step 5. Otherwise, continue to Step 4. Step 4. Generate a new mesh¯ (k+1) where i = 1, 2, . . . , N . k = k + 1, return to Step 2.

A posterior error analysis
Then, based on this interpolation function (4.1), we can derive the following a posteriori error estimate for the discrete scheme (3.1).
Proof By the definition of the differential operator L, we have For the first term of (4.3), we have where condition (1.4) is used.
Obviously, the second term of (4.3) satisfies the following estimation: For the third term of (4.3), we have where we have used the assumptions of a ij (x) and Theorem 3.1.

Numerical results and discussion
In Sect. 5.1, we first design a numerical method to obtain the parameter vector − → θ defined in (3.1). Then, numerical experiments are presented in Sect. 5.2 to demonstrate the validity and efficiency of the presented adaptive grid method. In the numerical experiments below, we shall take C 0 = 1.1.
For comparison purpose, the presented finite difference scheme (3.1) is also studied on a piecewise-uniform Shishkin mesh, which is constructed as follows [20]: Let N be divisible by 3 and λ > 0 be a mesh parameter. Fixed mesh transition points σ 1 and σ 2 are as follows: Then the mesh is obtained by dividing each of the intervals [0, σ 1 ], [σ 1 , σ 2 ], and [σ 2 , 1] into N/3 equal subintervals.

A calculating method of initial values
At first, from the first Robin boundary conditions of problem (1.1), we have Then, for given − → θ , combining with (5.1), we can use the presented adaptive grid method to calculate the numerical solution − → U i of (2.1). Obviously, the numerical solution − → U i depends on the value of − → θ . Finally, in order to obtain the approximation value of − → θ , we construct the following optimization problem: Since the above objective function Z( − → θ ) is an implicit function above variable − → θ , we choose the Nelder-Mead simplex method [14] to solve the above nonlinear optimization problem (5.2).

Example 1: a linear problem
In this section, we give a linear test problem to illustrate the theoretical result of the presented adaptive grid method with Robin boundary conditions Here, we choose f 1 (x) and f 2 (x) to agree with the following exact solutions: Recalling that − → U N i is the solution of the discrete scheme (3.1), we calculate the errors where x i are the points in the final mesh generated by the above iteration algorithm. Rates of convergence are calculated by For different values of ε i (i = 1, 2) and N , we use the Nelder-Mead simplex method to calculate the above nonlinear optimization problem (5.2) with (2, 2) T as the initial value of − → θ . Next, for ε 1 = 2 -9 , ε 2 = 2 -2k and ε 2 = 2 -9 , ε 1 = 2 -2k , where k = 0, 1, . . . , 10, we use the presented adaptive grid method to calculate the test problem with N = 32, 64, 128, 256, 512, 1024, and the errors in the maximum norm, the rates of convergence, and the number of iterations K are given in Tables 1-2. Besides, to illustrate the advantages of the adaptive grid method, we also compare the numerical results computed on an adaptive grid to those obtained on a Shishkin mesh, see Table 3. It is shown from these results that the convergence order of the presented adaptive grid method is first order and the accuracy of adaptive grid method is higher than that of Shishkin mesh method. For ε 1 = 10 -6 , ε 2 = 10 -8 , and N = 64, Fig. 1 shows the changing process of mesh movement with the number of iterations K . The right-hand part of this figure is labeled with the value of C 0 for which the stopping criterion (3.3) becomes an equation. It is shown from Fig. 1 that the solution of Example 1 has a boundary layer at x = 0.

A nonlinear problem
To demonstrate that our presented adaptive grid method can be successfully applied in a nonlinear setting, consider the following system of quasi-linear convection-diffusion equations: with Robin boundary conditions u 1 (0) + ε 1 β 1 u 1 (0) = s 11 , u 1 (1) + ε 1 γ 1 u 1 (1) = s 21 , (5.9) u 2 (0) + ε 2 β 2 u 2 (0) = s 12 , u 2 (1) + ε 2 γ 2 u 2 (1) = s 22 . (5.10)  As far as we know, the layer-adapted meshes and adaptive grid approaches for a single singularly perturbed quasi-linear two-point boundary value problem have attracted much attention; see [19,21] and the references therein. However, the adaptive grid method for the system of quasi-linear convection-diffusion equations (5.7)-(5.10) is not found in the literature to the best of our knowledge. Thus, in this paper, we will also use the presented

Figure 2
Mesh movement with ε 1 = 10 -4 , ε 2 = 10 -3 , and N = 64 for Example 2 gularly perturbed initial value problems. It also led to the development of an efficient adaptive grid method for solving the original problem.
Our numerical experiment has shown that the presented method can achieve first-order convergence. It should be pointed out that the ideas presented in this paper can be used to deal with the other high-order singularly perturbed differential equations with Robin boundary conditions.