Research Article | Open Access

Ke Li, Dali Guo, Yunxiang Zhao, "Convergence Analysis of Schwarz Waveform Relaxation for Nonlocal Diffusion Problems", *Mathematical Problems in Engineering*, vol. 2021, Article ID 5574403, 25 pages, 2021. https://doi.org/10.1155/2021/5574403

# Convergence Analysis of Schwarz Waveform Relaxation for Nonlocal Diffusion Problems

**Academic Editor:**Qiuye Sun

#### Abstract

Diffusion equations with Riemann–Liouville fractional derivatives are Volterra integro-partial differential equations with weakly singular kernels and present fundamental challenges for numerical computation. In this paper, we make a convergence analysis of the Schwarz waveform relaxation (SWR) algorithms with Robin transmission conditions (TCs) for these problems. We focus on deriving good choice of the parameter involved in the Robin TCs, at the continuous and fully discretized level. Particularly, at the space-time continuous level, we show that the derived Robin parameter is much better than the one predicted by the well-understood equioscillation principle. At the fully discretized level, the problem of determining a good Robin parameter is studied in the convolution quadrature framework, which permits us to precisely capture the effects of different temporal discretization methods on the convergence rate of the SWR algorithms. The results obtained in this paper will be preliminary preparations for our further study of the SWR algorithms for integro-partial differential equations.

#### 1. Introduction

Fractional kinetic equations have been proved particularly useful to model the anomalous slow diffusion (subdiffusion) [1]. Subdiffusive motion is characterized by an asymptotically long-time behavior of the mean square displacement, such aswhere is the anomalous diffusion exponent. The process is usually referred to as subdiffusive when . For the ordinary (or Brownian) diffusion (this corresponds to ), from a continuous point of view the diffusion process is described by the well-known diffusion equation , where represents the probability density of finding a “particle” at space and at time . It turns out that the probability density function that describes anomalous (sub) diffusive particles follows the fractional diffusion equation [1–4]:where is the Riemann–Liouville fractional derivative defined by (6) below.

Subdiffusive equation (2) is the fundamental model for other anomalous diffusion observed in many fields, such as the transport of fluid in porous media [5], diffusion on liquid surfaces [6], and turbulent flow [7]. Other applications of fractional differential equations can be found in [8–17]. Very recently, the so-called fractional cable equation was introduced for modeling the anomalous diffusion in spiny neuronal dendrites [18–20]:where . This equation models spatial-temporal dynamics of the membrane potential along the axial direction of an approximately cylindrical segment of a nerve cell, which takes into account anomalous transport in an external field. In (3), denotes the diameter of the cable, is the equilibrium potential of the membrane, is a membrane constant and is a constant independent of the membrane area, and is an axial material constant, associated to the given cable. By taking , and then by letting , cable equation (3) can be transformed to with . Dropping the tilde superscript, we finally get

Unlike the integer order PDEs, a simple and explicit analytical solution to a fractional PDE is usually unavailable and therefore numerical methods play an important role in investigating the properties of these equations. Moreover, compared to the integer order problems, numerical computation for the fractional PDEs is much more time consuming, because we need to treat convolutions with singular kernels at every time grid. Here, our interest is to apply the Schwarz waveform relaxation (SWR) algorithms to solve the fractional diffusion equations. These algorithms belong to the widely used domain decomposition methods, but with completely new implementation strategy for time-dependent problems: the classical strategy employs the alternating or parallel Schwarz methods to the elliptic PDEs which result from semidiscretizing the time-dependent PDEs in time by some implicit time-integrator [21,22], but in the SWR framework one directly decomposes the space domain into many subdomains and then solve each subproblem simultaneously or alternately. This provides flexibility for treating different subproblems numerically differently with an adapted procedure for each subdomain (see [23,24] for the original idea of the SWR algorithms). There are also other efficient algorithms for handling the fractional PDEs, such as the interesting parallel-in-time algorithms studied in [25–29].

For integer order differential equations, the SWR algorithms have been investigated widely and deeply by Gander and his colleagues [30–35]. However, much less results are known about the algorithms applied to fractional order PDEs; actually, according to our best knowledge, this is the first time that the SWR algorithms are analyzed for fractional problems. In this paper, we analyze the algorithms with Robin transmission conditions (TCs) for the aforementioned fractional cable equation, which, as we will show in Section 2, represents a large class of fractional problems that received considerable attention in the literature. The Robin TCs contain a free parameter, which has a significant influence on the convergence rate of the SWR algorithms; hence, finding a good choice of the Robin parameter is one of the top-priority matters. For the integer order PDEs, optimizing the Robin parameter is deeply studied in [30,35–37] and the main finding is that a good parameter can be computed by solving one or two simple nonlinear equations. This is the so-called equioscillation principle, and nowadays, it is well-understood and widely used in the study of the SWR algorithms.

However, as we will show in Section 3, for the time-fractional diffusion problems, the Robin parameter predicted by the equioscillation principle is unsatisfactory and there exists a much better choice for specified problem/discretization parameters. To make the SWR algorithms practical in real computation, we also perform a convergence analysis at the fully discretized level in Section 4, in the framework of convolution quadratures [38–40]. This framework allows for a convenient treatment of the convergence analysis based on discrete Laplace transform. We validate the theoretical conclusions by numerical results in Section 5 and we finish this paper with some concluding remarks in Section 6.

#### 2. The Model Problem and the SWR Algorithms

Motivated by the fractional cable equation, which includes fundamental equation (2), our model studied throughout this paper is the following initial value problem:where and . This model gets considerable attention in recent years and we can refer to [19,20] for the construction of solutions based on Green’s functions. For any , we denote by the Riemann–Liouville fractional derivative:where is the gamma function . From [41] (Chapter 2), for any it holds that . Hence, the Riemann–Liouville fractional derivative builds a “bridge” from to .

In the field of fractional PDEs, there exists another definition of the fractional derivative, the Caputo derivative:

The two definitions are linked by the following relation [41], (pp. 62–77):

Thanks to this connection, we shall only consider the Riemann–Liouville derivatives, since the results can be directly generalized to problems with Caputo derivatives.

There exists another class of fractional PDEs, which are also widely studied in the literature [41,42], the fractional heat equations:

From [41] (pp. 65–75, Chapter 2), we know that (9a) can be rewritten as

Hence, the analysis for the model problem (5) can be also applied to (9a).

We now introduce the Schwarz waveform relaxation algorithms for (5). We divide the spatial domain into two subdomains, , with and . Here, denotes the size of overlap. Then, the SWR algorithms for (5) arewhere is a free parameter, is iteration index, and . It is clear that the proposed SWR algorithm does not directly depends on the orders of the differential equations. So, the algorithm is applicable to variable-order time-fractional PDEs. However, the convergence analysis will be different.

#### 3. Convergence Analysis at the Continuous Level

We are now in a position to analyze SWR algorithms (10). In this section, we perform an analysis at the continuous level, and in the next section, we analyze the algorithm at the fully discretized level. For any real-valued function , we denote by the Laplace transform of , defined by

Let . Then, it is easy to know that satisfies

Applying Laplace transform to the error equation (12), we get

The solutions can be expressed as . Then, using the fact that do not increase exponentially at infinity, we get

Substituting (14) into the transmission conditions in (13), we getand this gives and , where

Hence, . Mathematically, we want , which leads to the following min-max problem:

Lemma 1. *Assume , and . Then, the argument defined by (16) is analytic in the right half of the complex plane.*

*Proof. *Let with and . Then, we haveThis implies that, for and the function defined in (13) is analytic in the right half of the complex plane, since the quantity under the square root avoids the negative real axis. For , the denominator does not vanish in the right half of the complex plane; hence, is an analytic function. Since the product of two analytic functions is still analytic, is analytic.

Note that, for any with , the limit of as is zero. Hence, by using the maximum principle for complex analytic functions, we know that the maximum of is obtained on the imaginary axis, i.e., . Moreover, in a numerical computation, the angular frequency cannot be arbitrarily large or small; there is a range of frequencies, which can be represented on a time grid with size and a length of time window , as . Hence, it is reasonable to define , which we call convergence factor of the SWR algorithm at the space-time continuous level, as

The subscript “” appearing in and denotes “continuous.” Clearly, an ideal is the minimizer of the convergence factor , i.e.,

For any , it holds that

Thus, for any , we can rewrite as , where

This gives

Clearly, for any , it holds . Thus, min-max problem (20) is equivalent towhere

##### 3.1. Towards a Good Choice of : The Equioscillation Principle

A well-understood formula for choosing the free parameter is based on the so-called equioscillation principle, established recently by Gander and Halpern [30,35]:where and denotes the local maximizer of . The solution of (25), namely, , can be obtained by solving the following two nonlinear equations:together with a constraint (equations (26) may have several solutions and only the one satisfying this constraint corresponds to the solution of (25)). Solving such a nonlinear problem is quite easy from a numerical point of view. For integer order problems, such as the advection-reaction-diffusion equations, formula (25) gives the best choice of ; see [30,35]. In Figure 1(a), we have given an illustration of the equioscillation principle for the min-max problem (24a)-(24b).

**(a)**

**(b)**

**(c)**

For time-fractional PDEs, problems exist for the equioscillation principle in twofold:(1)It is unclear for the moment whether equation (25) is well defined or not. Is there a unique positive solution? Is the solution located in the relevant interval ? For integer order time-dependent PDEs, the answers for these questions are positive, at least in the asymptotic sense, i.e., , , and small; see [30,35]. However, in the fractional situation, we know nothing about (25).(2)More importantly, the parameter predicted by the equioscillation principle is unsatisfactory. For example, with problem parameters listed in the title of Figure 1(a), the parameter results in a convergence factor around 0.842, which implies that the SWR algorithm using for the Robin parameter converges slowly.

Since the applicability of the equioscillation principle to time-fractional problems is unclear, we introduce in Lemma 4 a new formula for computing . In Figure 1(b), with , , , , and two groups of and , we plot the function for , where we can see that compared to the new parameter leads to much smaller convergence factor and that with the new parameter the function does not equioscillate.

##### 3.2. Towards a Good Choice of : A New Formula

Because of the high complexity of the arguments , it is difficult to present a closed formula for the solution of min-max problem (24a) and (24b). In what follows, we focus on deriving a reliable approximate solution. We need two auxiliary lemmas.

Lemma 2. *Let , and*

Then, it holds that , where .

*Proof. *Clearly, we only need to consider the case , since for the function is a constant and therefore .

Let and . We then have and . Suppose has a positive root located at , i.e., . Then, we haveFor any , it holds thatand this, together with and , givesHence, if corresponds to a local extremum of , it must be a minimizer.

Lemma 3. *Let , , and with being the function defined in Lemma 2. Then, for all , we have .*

*Proof. *From (22), . Since ,Hence, is a decreasing function of ; this givesBy noticing , we haveHence, from Lemma 2, we get and this completes the proof.

Now, for min-max problem (24a) and (24b), by using Lemma 3, we have

For any , this leads to the following inequality:

In (34), we have used and . Let and . Then, the solution of the following min-max problem can be regarded as a reliable approximate solution of (24a) and (24b):where

We present the solution of (36a) and (36b) in a very general setting: we only assume and . The concrete proof is presented in Appendix.

Lemma 4. *Parameter optimization at the continuous level: assume and . Let , , , and be the solution of min-max problem (36a) and (36b). Then,*(a)*If , we have* *where and is defined by* *The argument is the unique root of .*(b)*If , we have* *where , , , , and is the unique solution of*(c)*If , then* *where is the unique solution of the equation (38b).*

Theorem 1. *Let be the function defined in Lemma 2 and *

Then, for specified , , and , the parameter involved in the Robin TCs can be chosen as , where , given by Lemma 4, is the solution of min-max problem (36a) and (36b) and denotes the size of overlap. The quantities and in (36b) can be explicitly written down, as

With , the convergence factor defined by (19) can be bounded bywhere is given by (37a), or (38a), or (39b).

*Proof. *bound (41) can be derived from (34) and (36b):

(34), (36b).

Hence, it remains to prove that the arguments and defined by (36b) can be explicitly expressed as (40). To this end, we defineClearly, for any and , the function is increasing for .

*Case 1. *. In this case, by noticingwe know that the function is increasing for . Hence, we have

*Case 2. *. In this case, from (43), we get , with . We have and for , it holds that if and ; therefore, and satisfy (40).

We now consider the case . Suppose the function has a local extremum located at . Then, we have . Hence, it holds that . We claim . Indeed, from (42), it is easy to get , and this, together with , implies . For , a tedious but routine calculation yieldsThis implies that has no local maximum(s) and thus also satisfy (40).

We next comment in two remarks the parameter optimizations studied here for the time-fractional problems. First, in Remark 1, we show that the analysis of finding the solution of min-max problem (36a) and (36b) is different from the analysis existing in the literature. Then, in Remark 2, we show that the parameter given by Lemma 4 is close to the best one and that we can get through numerical optimizations.

*Remark 1. *About min-max problem (36a) and (36b): for the integer order reaction-diffusion equation (or more generally, the advection-reaction-diffusion equation ), the related min-max problem iswhere are constants related to the problem/discretization parameters (see [35] for more details). An important property of this min-max problem is that, as we already mentioned above, the solution is determined by the equioscillation principle [30,35]. Although the basic mathematic tools used for analyzing (46) and (36a) and (36b) are the same (such as monotonicity, continuity of real-valued functions, and roots of quadratic/cubic polynomials), the concrete procedures are different. This difference can be seen from two points: (1)The two functions, in (40) and in (36a) and (36b), are different (2)Because of this difference, the solution of min-max problem (36a) and (36b) is not always determined by the equioscillation principleFor min-max problem (36a) and (36b), there are two reasons that the equioscillation principle does not hold. First, the local maximizer does not exist in the interval (see Figure 2(a)). Second, more typically, even though the local maximizer exists, (see Figure 2(b)). What we observed in these two subfigures does not happen occasionally, and it seems an inherent feature for the time-fractional problems. Precisely, as will be shown in our forthcoming paper, in the asymptotic sense, i.e., , and small, the equioscillation property always does not hold for min-max problem (36a) and (36b) when and .

**(a)**

**(b)**

*Remark 2. *About the parameter : the solution of min-max problem (36a) and (36b) depends on and , which are related to (problem parameters) and (discretization parameters). Hence, we denote bywhich we call “parameter formula” and its concrete expression is given by Lemma 5 as we have seen. Note that is not the solution of original min-max problem (24a) and (24b); it is only an approximate solution. Hence, it is important to present evidence to show that is really a good approximate solution. Our evidence is based on numerical optimization:(1)Find the best parameter through repetitive numerical optimizations for original min-max problem (24a) and (24b):(i)For , compute by the global optimization toolbox in Matlab, where (since sign , it suffices to consider ) and(ii)Compare the values , and the minimal one corresponds to .(2)Compute and then compare the values of , i.e., the convergence factor of the SWR algorithm with (the approximate solution) and obtained through numerical optimizations.Let , , , and . Then, in Figure 3, we plot (Figure 3(a)) and (Figure 3(b)) for , where a small step-size, , is used in the above numerical optimization procedure to get . In Figure 3(c), we plot , i.e., the error between the convergence factors using the two parameters. This figure clearly shows that the quantity , given by Lemma 4, is really a good approximate solution of original min-max problem (24a) and (24b).

**(a)**

**(b)**

**(c)**

#### 4. Convergence Analysis at the Discrete Level

To make the SWR algorithms practical for real computation, we need an analysis at the fully discretized level. To this end, we first make a transform for (5), which avoids derivatives. Integrating (10) from 0 to , we get