Parameter-uniform convergence analysis of a domain decomposition method for singularly perturbed parabolic problems with Robin boundary conditions

We construct and analyze a domain decomposition method to solve a class of singularly perturbed parabolic problems of reaction-diffusion type having Robin boundary conditions. The method considers three subdomains, of which two are finely meshed, and the other is coarsely meshed. The partial differential equation associated with the problem is discretized using the finite difference scheme on each subdomain, while the Robin boundary conditions associated with the problem are approximated using a special finite difference scheme to maintain the accuracy. Then, an iterative algorithm is introduced, where the transmission of information to the neighbours is done using a piecewise linear interpolation. It is proved that the resulting numerical approximations are parameter-uniform and, more interestingly, that the convergence of the iterates is optimal for small values of the perturbation parameters. The numerical results support the theoretical results about convergence.


Introduction
Singularly perturbed problems (SPPs) arise in the mathematical modeling of several practical problems in engineering and applied mathematics, for example, in describing the theory of gyroscopes [25], in studying linear spring-mass system without damping but with forcing and a small spring constant [31], in variational problems in control theory [10], etc., where the higher order derivative appears multiplied by a small parameter. The fascinating aspect of SPPs is that their solutions have boundary and/or interior layers (the regions in which there are steep gradients). Because of this reason, classical numerical schemes are not appropriate to accurately and efficiently solve these problems. This has led to the development of special numerical approaches for SPPs. The developed numerical methods for SPPs are called as parameter-robust or parameter-uniform or uniformly convergent, meaning that the convergence of the method is independent of the involved parameter. They are developed considering fitted operators, fitted meshes, and domain decomposition approaches [1, 2, 9, 11, 14-16, 26, 28-31, 35]. Further, some recent advancements in finite difference methods can be seen in [3][4][5][6].
Consider the following SPP with the Robin boundary conditions where D = (0, 1), a(x, t) ≥ α > 0 on , and 0 < ε ≤ 1 is known as the perturbation parameter. It is known that problem (1)-(3) has a unique solution exhibiting boundary layers near x = 0 and x = 1 [12,22,37]. Such problems arise in the modelling of certain types of bioswitches [41]. Singularly perturbed problems similar to (1)-(3) with Dirichlet-type boundary conditions have been studied extensively in the literature (see [7,18,20,24,27,34,38] and the references therein). However, there are only a very few studies of such problems with Robin boundary conditions (RBCs) [12,13,22,37]. Note that all these studies are based on the fitted mesh approach.
Domain decomposition is shown to be a very versatile and effective approach for solving problems in ordinary and partial differential equations [8,17,32]. The original approach dates back to 19th century [36]. However, the approach flourished more in the last three decades. Schwarz waveform relaxation (SWR) is a special approach that is very popular for solving time dependent problems [32]. In this technique, the spacetime domain is split into several subdomains. Then, on each subdomain across the whole time interval the solution is computed and the exchange of time-space boundary values takes place. We can treat the subdomains differently using this approach. Also, we have the flexibility to locally address any singularity in the solution. Further, a non uniform mesh can be avoided with this approach. Lastly, one can speed-up the computations by implementing it on parallel computers. Domain decomposition methods for SPPs have been developed in [19-21, 33, 34, 39, 40] and the references therein. Note that in all these works model problems with Dirichlet boundary conditions are considered. Moreover, we are not aware of any studies involving the domain decomposition approach for SPPs with Robin boundary conditions.
Hence, the objective of this paper is two-fold: first, to introduce a domain decomposition method of SWR type to numerically solve problem (1)-(3), and second, to present a parameter-uniform convergence analysis of the introduced method. We consider a space-time partitioning of the original domain using the Shishkin transition point. The PDE associated with problem (1)-(3) in each subdomain is discretized by the finite difference scheme, while the Robin boundary conditions of problem (1)-(3) are approximated by a special finite difference scheme to maintain the accuracy. Then an iterative algorithm is introduced, where transmission of the information to the neighbours is done using the piecewise-linear interpolation. We discuss parameter-uniform convergence of the constructed method using auxiliary problems that separates the discretization and iteration errors. Note that the convergence analysis of domain decomposition methods for problems with Dirichlet boundary conditions cannot be straightforwardly applied to the present method due to the presence of the Robin boundary conditions. Firstly, we require a different definition of auxiliary problems to handle the Robin boundary conditions. The boundary and initial conditions of the auxiliary problems in earlier works are simply defined using the exact solution of the problem, but that is not entirely true for the present problem. Secondly, we require a more complex constant coefficient problem (in Theorem 3.2) and altogether different barrier functions (in Theorems 3.2 and 3.3) for bounding the error between the auxiliary solution and the Schwarz iterates. It is proved that the resulting numerical approximations are parameter-uniform and, more interestingly, that the iteration convergence is optimal for small ε. Finally, numerical results are provided to support the convergence theory.
The rest of the article is organized in different sections as follows. Section 2 provides a priori bounds for the solution derivatives, and Sect. 3 includes the development of the SWR method for problem (1)-(3). Section 4 includes the convergence analysis of the developed method. Finally, in Sect. 5, we present the numerical results for two test examples confirming our convergence theory.
Notation: C denotes a generic positive constant which is independent of the parameters N , M, ε, and k. For v ∈ C( ), let us define v i, j = v(x i , t j ). Here, . Q is used to define the maximum norm on a closed and bounded set Q, and . Q N ,M is used to define the corresponding discrete maximum norm.

Asymptotic behavior
We discuss a priori bounds on the solution derivatives of (1)- (3). The existence of a unique solution of problem (1)-(3) is established in the following lemma. 1), and the compatibility conditions up to the first order are true. Then, problem (1)-(3) has a unique solution u ∈ C (2+β,1+β/2) ( ).
In the following lemma, crude bounds for the derivatives of u are given. 1), and the compatibility conditions up to the second order are true. Then, problem (1)-(3) has a unique solution u ∈ C (4+β,2+β/2) ( ). Moreover, it holds Proof The proof follows from the arguments in [37,Theorem 5].
The above bounds are not enough for the convergence analysis in Sect. 4. We require to split u into regular and singular parts, and bounds for their derivatives as given in the following lemma.

Proof
The lemma can be proved following the arguments in [37,Theorem 6].

The domain decomposition method
We partition into p = D p × (0, T ], p = , m, r , where D = (0, 2ρ), D m = (ρ, 1 − ρ), D r = (1 − 2ρ, 1) with the parameter ρ defined by Here, ρ is chosen as the Shishkin transition parameter [9,26], which enables us to have fine meshes (for small ε) on and r (the subdomains corresponding to the boundary layers). Further, we discretize each subdomain p with a uniform mesh in both the spatial and temporal directions. On each spatial subdomain When the domain is obvious, for convenience, the term p will be dropped from the subscript in Here, M and N are discretization parameters in time and space directions respectively. Defining We discretize (2) as follows Note that we have approximated the boundary conditions using a special discretization scheme. If we had used the standard upwind scheme for the discretization of the boundary conditions we would have obtained only first order accuracy. Therefore, we use a special discretization scheme for the boundary conditions, which is based on improving the truncation error to maintain second order accuracy.
To compute as the approximate solution of problem (1)-(3) the algorithm is defined as follows. We start the iterative process with U [0] (x i , t j ) as the initial approximation, defined as follows: We then compute U [k] by The iterations are performed until t . Now using the arguments in [37,Theorem 7] we can prove the following maximum principle for The stability of the numerical scheme is given by the following lemma.

Lemma 5 For any mesh function Z p , we have
Then, it is easy to verify that Hence, from Lemma 4, the proof follows.

Convergence analysis
We now establish that the method gives parameter-uniform approximations to the solution of problem (1)-(3). Further, we prove that the iterative process converges optimally for small ε. We consider the following auxiliary mesh function where U p , p = , m, r , are such that Here and (1 − ρ, t j ). By the triangle inequality, we get The following lemma gives the bound of the first term on the LHS of the inequality in (11).

Lemma 6
Suppose u is the solution of (1)- (3) and U is the auxiliary mesh function defined in (10). Then On the other hand, from (1)-(2) and (8) it follows that Using Taylor expansions and (4) we get and where we have used Taylor expansions to get (14). Since h ≤ C √ εN −1 ln N , using (4) it follows from (14) that Hence, by Lemma 4 with C( t + (N −1 ln N )
The following notation will be used in the next two theorems.
We will show the parameter-uniform convergence of the method in two cases: when ρ = 2 ε α ln N and ρ = 1/4. In the next theorem, for ρ = 2 ε α ln N , we first obtain a bound of ||U − U [1] || N ,M and then combine it with Lemma 6 to get a bound of ||u − U [1] || N ,M . So, this result demonstrates that only one iteration is enough to achieve the required accuracy. [1] is its approximation obtained after the first iterate of the proposed method. Then, for ρ = 2 ε α ln N , it holds ||u − U [1]

Theorem 1 Suppose u is the exact solution of (1)-(3) and U
For an arbitrary mesh function Y , let us define [1] for t j ∈ M t (15) with and ξ i , i = 1, 2, are as follows The solution of (15) is given by which is monotonically increasing and satisfies χ ≥ 0 in N ,M . Thus, one gets and Hence, applying Lemma 4 to χ ± , we obtain This implies that [1] ϕξ [1] Cϕ Similarly, it is Next, consider U m − U [1] m which satisfies and since the mesh points (ρ, t j ) and (1 − ρ, t j ) belong to N ,M and N ,M r , respectively. Therefore, Lemma 4 gives Combining (17), (18), and (19), we obtain Since , by Lemma 6 one gets μ ρ ≤ C( t + (N −1 ln N ) 2 ). Hence, using (20) and Lemma 6 in (11) the proof is complete.
After repeating the above arguments one gets We simplify this to get Thus, it is Further, since (2ρ, t j ), ( N −1 ln N ) 2 ). Since μ [1] ≤ C, using (26) and Lemma 6 in (11) the proof is complete.

Numerical results
Numerical results considering a couple of examples are given that verify the convergence theory.
where f , g , and g r are such that The solution plots for various values of ε are given in Fig. 1. Note that the boundary layers are close to x = 0, 1. Taking the tolerance error tol = N −2 , we compute the approximate solution and denote it by U N , t . We then evaluate parameter uniform errors and convergence orders as follows    Table 2 gives the iterations needed for convergence; from this, we observe that one iteration is enough to get the desired results for very small values of ε.

Example 2 Consider
The solution plots for various values of ε are given in Fig. 2. Note that the boundary layers are close to x = 0, 1. Since the solution is unknown, we compute another where U 2N , t/4 is obtained taking in each subdomain 2N + 1 points in space and t/4 mesh size in time, but using the same subdomain parameter ρ as is considered for U N , t . The error plots for various values of ε are given in Figs. 3 and 4 for Examples 1 and 2 respectively.
The computed errors E N , t ε , E N , t , and uniform convergence rates R N , t for Example 2 are tabulated in Table 3, showing the uniform convergence of the method. Table 4 gives iterations that are performed to achieve the convergence for Example 2.
One can see that only one iteration is needed to get the solution up to the desired accuracy for very small values ε.

Conclusions
An SWR technique to solve singularly perturbed parabolic reaction-diffusion equations with Robin boundary conditions is developed and analyzed in this paper. The original domain is initially divided into three overlapping subdomains. The problem is discretized using the backward difference and central difference schemes for the time and space derivatives, respectively. It is shown that the proposed scheme is unconditionally stable. Error analysis is also discussed in this work, and it is demonstrated that the proposed approach is uniformly convergent with order almost two in space and one in time. Furthermore, it is shown that for small values of ε, one iteration is sufficient to achieve the specified accuracy. The idea discussed in this paper can also be extended to singularly perturbed semilinear differential equations having boundary conditions of Robin type. Further, it is our intention in the future to extend the SWR technique to higher dimensional singularly perturbed problems.