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 getwhere . For any , denotes fractional integral
4.1. Derivation of the Discrete SWR Algorithms
To get discrete algorithm, we apply the centered finite difference formula to (49), which giveswhere , , , and . To simplify the convergence analysis, we assume , i.e., the size of overlap equals a single space mesh size. Then, the artificial boundary interfaces (resp. ) corresponds to (resp. ) and the Robin TCs are discretized as
Here, we discretized the spatial derivatives in the Robin TCs as
Note that, (51c) is a first-order discretization, and therefore, a question arises naturally: does this simple discretization affect the accuracy of the converged solutions? Following the proof given in [43], Lemma 1, this worry can be removed.
In (51a), we need to deal with convolutions, aswhere (or ) and is not known beforehand but are computed as the computing proceeds from time step to time step. A naive implementation would require operations and memory for per expansion coefficients over time steps. In [44], fast convolution algorithms are constructed, which implements, over time steps, with operations and active memory.
The starting point of the fast convolution algorithms is to replace the kernel function in (52) by its inverse Laplace transform, along some contour :where denotes the Laplace forward transform of the kernel function . The contour is a curve in the sector oriented with an increasing imaginary part and going to infinity with an acute angle to the negative real half-axis, so that decays rapidly for growing along . In (53), reversing the order of integration gives
The function is the solution of the following linear scale ODE at :
Since is only available at discrete grids , it is reasonable to apply some time-integrator, e.g., the linear multistep method, to (55) to get at :
Multiplying (56) by and summing over from 0 to , we obtainwith the generating (formal) power series and . Solving from (57), we getwhere is called the quotient of the generating polynomial of a linear multistep method. Clearly, the solution of (56) is the -th coefficient of the power series , that is, . By substituting this into the right-hand side of (54), we get
Now, we expand in the formal power series, as
Then, the -th coefficient of is given by
Suppose the time-integrator applied to (55) is an -stable method with order . Then, from [40], Theorem 2.1, we know that (61) is an approximation to (52) and
The last step is to determine the weights . For general case, we have to rely on numerical calculation to get (see [44] for more details). For some special cases, such as the linear -method () and the two-step backward difference formula (BDF2, ), these weights can be written down, as
With the quadrature weights in our hands, we denote the discrete convolution (61) by , that is,
Then, the scheme of the fully discrete SWR algorithms on the -th subdomain is
Together with boundary conditions (51b) at .
4.2. Convergence Analysis
Since the quadrature weights are the coefficients of the formal power series of , the discrete Laplace transform of the sequence defined by (64) is , where and . Here, with and . Define the error sequence as . Then, the discrete Laplace transform of , which we denote by , satisfiesi.e.,where
Define
Then, interface conditions (51b) are transformed to
From (67), the arguments can be expressed as
Lemma 5. Let , , , , and . Then, for with and , we have(1)if (linear -method):(i) and , if ;(ii) and , if ;(2)if (BDF2 method): , .
Proof. Let and . We have, , , and . The proof for the two time-integrators is given as follows:(1)For the linear method: . In this case, we have where and . For , the bounds for and can be obtained easily. For , Clearly, it holds that and this implies Moreover, since and , we have This gives the bounds of and for .(2)For the BDF2 method: . In this case, we have Hence, . For , by noticing it is easy to get .
Lemma 6. Let , , and . Then, for the linear -method with and the BDF2 method, the arguments with are analytic for and satisfy and .
Proof. For , from Lemma 5, we have . Then, for , it holds that for , where is defined by (68). This implies , and thus, the composite functions are also analytic. Since , we have ; this, together with , gives .
Now, by using the boundedness conditions at infinity for each subsystem, the constants and in (71) can be determined as and for the first and second subsystems, respectively. We then have and . By using the transmission conditions in (70), we get the recurrence relation and , i.e.,where the subscript “” denotes “discrete” and
Mathematically, we expect , which leads to the following min-max problem:
Let . Then, from Lemma 6, we know . For , we have . This implies that, it is sufficient to consider in (79).
For , the quantity approaches to the origin, as approaches to in the complex plane. Hence, approaches to 1 as . This implies that , independent of .
For , we know from Lemma 5 that the origin is kept away from , and therefore, it is possible to expect . However, the region mapped by is a domain with curved boundary, which makes the analysis of solving (79) extremely complicated. Here, we take up the idea proposed by Vandewalle and Gander [45]: finding a rectangular region,which lies in the right-hand side of the vertical line (i.e., ) and covers the region mapped by ; see Figure 4 for an illustration. Then, regarding as a single complex variable, solve (79) on this rectangular region.
For the linear -method with and the BDF2 method, if , by using Lemma 5 and Lemma 7 given below, the boundaries of the rectangular region , i.e., , , and , can be expressed explicitly, aswhere are defined by (87a). The quantities , , and are defined by
The arguments and , appearing in (81b), are given bywhere and .In general, we have to rely on numerical optimizations to obtain these three quantities, i.e., , , and .
Lemma 7. . Assume and . Let and , where and are given by (81a) and (81b). Then, we have for .
Proof. With the arguments , , and given by (81b), we define . Then, the proof is divided into two steps: if and if . The function is defined by (68).
Step 1. for . Let . Then, we can rewrite aswhere . For any , the functions are defined by where . Since , we have . For fixed , because is an increasing function of , is also an increasing function of . Hence, for , . Define . Then, = for . We haveHence, for , it holds that and this implies that the function attains its maximum at . We therefore get the upper bound of , i.e., , as shown in (81b).
It remains to consider the lower and upper bounds of . For and , we have . Hence, for , it holds that and this givesThe first result gives the lower bound of , i.e., . With the function , the argument can be represented as with . We haveHence, attains its global maximum at ; this gives , the upper bound of .
Step 2. for . Let with and . Then, , whereClearly, it holds and . Moreover, some simple algebra showsHence, and from the second inequality given in (87b), we know that attains its maximum at ; this gives the upper bound of , i.e., . Similarly, we haveCombining the analysis in Step 1 and Step 2, we complete the whole proof.
With the rectangular region covering , we haveand based on this relation, it is reasonable to consider the following equation:
To analyze the solution of this min-max problem, we define
Then, for , we defined two functions and as
For , these two functions are defined by
Theorem 2. Parameter optimization at the discrete level: let , and . Then, the parameter involved in the discrete SWR algorithms can be chosen as , where , the solution of min-max problem (90), is given bywhere , , and are the unique root of . With the parameter , the convergence factor at the discrete level satisfies , where is given by (78a).
Proof. Let . Then, it holds that . For , a routine calculation yieldsHence, the solution of (90) should satisfy , since otherwise, if (resp. ), we can further minimize by increasing (resp. decreasing) .
We have , and therefore, for fixed , it holds that . Hence, Another routine calculation yieldsThen, for fixed and , if we vary from 1 to , we first meet , and then , and then . Hence, for fixed , the function does not have local maximum(s) for . This implies that for given the global maximum of on the region is obtained at the four corners, i.e.,With the functions defined by (91b) and (91c), we haveIt remains to analyze the monotonic properties of and . From (93), we know that and are, respectively, increasing and decreasing for . Similarly, is decreasing for and increasing for ; is decreasing for and increasing for . Hence, the functions and , defined by (91b) and (91c), are, respectively, increasing and decreasing for . This, together with some simple algebra, gives the solution of (96), i.e., , as shown in (92).
Remark 3. The optimization procedure at the discrete level consists of two steps:(1)For specified quotient and problem/discretization parameters, compute the boundaries of the rectangular region , i.e., , , and (2)Compute according to formula (92) and then let The quotient associated with a specified time-integrator determines the boundaries of , which further determines the parameter . To make a convergent SWR algorithm, should lay in and should be away from the origin, .
Remark 4. Large overlap size: if we use an overlap size with , the following min-max problem occurs:For and , the function has local maximums in the interior of , and therefore, the analysis is much more complicated compared to the case . The proof of Theorem 2 cannot be straightforwardly generalized to (97), and we will investigate this general min-max problem in other places.
4.3. An Application in the Field of Fractional Circuits
In the last part of this section, we show that the analysis performed at the discrete level can be directly applied to fractional diffusive (RC) circuits. In recent years, there is an increasing interest in the fractional order circuits, thanks to the fact that they can more precisely describe the dynamic behavior of the supercapacitors (also called ultracapacitors or electrochemical double-layer capacitors); see the recent papers [46,47] for more details.
For the model circuit shown in Figure 5, consisting of infinite number of series connected cell-circuits, the state (or modified-nodal-analysis) circuit equation is
(a)
(b)
where denotes the vector of the nodal voltage values, with the value of the fractional capacitor and the value of the resistance, as shown in Figure 5(a). The fractional derivative in (98) is defined in the Caputo sense (see (8)), and the argument is determined by the branches of the circuit shown on the right. Mathematically, (98) can be equivalently rewritten as , with the fractional integral defined by (50) and some other function . Clearly, this is a special case of (), which can be seen by letting and and then by replacing by and by there. Therefore, the semidiscrete SWR algorithm (51a) and (51b) corresponds to the waveform relaxation method for (98).
For an analysis of the continuous WR algorithm without time discretization, i.e., (51a) and (51b), we will also arrive at (78a), but withinstead of as shown in (78b). Hence, the rectangular region covering is a semi-infinite horizontal stripe, i.e., in (80). In this limit case, Theorem 2 is still applicable, but we need to replace by there and redefine aswhere and are the same quantities as we have stated in Theorem 2.
5. Numerical Results
In this section, we test and compare the performance of the parameters analyzed at the continuous and discrete level. Here, we restrict ourselves to the case . Since the parameter appearing in the continuous Robin TCs is transformed to in the discrete circumstance, we shall only mention throughout this section.
At the moment, for given problem/discretization parameters, we have two choices of , from the continuous analysis (i.e., Lemma 4) and from the discrete analysis (i.e., Theorem 2). Both are approximate values of the best one that we can make through numerical optimizations (as described in Remark 2). To compare the performance of these two parameters, in Figure 6 we show the convergence factor based on three choices of , where and are given in (78a) and (78b). Here, the problem/discretization parameters are chosen as , , , and . To get from Theorem 2, the boundaries of the rectangular region , i.e., , , and , are obtained through numerical optimizations (for the case , if we use the explicit expressions (81a)–(81c) to get these three quantities, the corresponding convergence factor looks similar as what we observed in Figure 6).
(a)
(b)
From Figure 6, we see that for the backward Euler (or BDF2) method, the convergence factors under the three choices of are almost equal. For a specified choice of , the convergence factors associated with the two time-integrators are almost equal. These can be regarded as theoretical predictions, and we now provide numerical results for validation. We consider the following equation with zero initial and boundary values:which is equivalent to model problem (5) after an integral transform. We choose , , and for the problem parameters, and and for the discretization parameters. For each experiment, the initial iterate is chosen randomly.
In Figure 7(a), we plot the convergence rates of the SWR algorithms, using from the continuous analysis. (For the parameters obtained from the discrete analysis, this error plot looks similar.) The three lines shown in this figure correspond to three time-integrators, trapezoidal (dot line), backward Euler (dot-dash line), and BDF2 (solid line). We see that (1) the convergence rates associated with backward Euler and BDF2 are very close and (2) both are significantly better than that of the trapezoidal method. The first observation is consistent with Figure 6, and the second one is consistent with our analysis at the discrete level, since for the trapezoidal method, we always have , independent of . Then, in Figure 7, we verify whether the theoretically analyzed parameter is close to the best numerical one or not. Precisely, we show the error obtained by running the SWR algorithm after 4 iterations using various values for the parameter in the Robin TCs. The choices predicted by the continuous and discrete analysis are indicated by a circle and a star, respectively. Clearly, both are very close to the best numerical one.
(a)
(b)
(c)
At the end of this section, we show in Figure 8 the effects of the discretization parameters and the length of time-interval on the convergence rate of the SWR algorithms using the backward Euler method as the time-integrator in the case of two and eight subdomains. For each and , we use two choices for the parameter involved in the Robin TCs, the one from the continuous analysis and the one from the numerical optimization. (For the parameter from the discrete analysis, the plots look similar.) For each and , the iteration number is measured when the error between the iterate and the reference solution is less than . The reference solution is obtained by applying the same discretization to (101) on the whole space-time domain. We see that the SWR method has very similar asymptotic convergence rate, depending on the discretization parameters and the length of time-interval, in the case of two and many subdomains. We can also see that the convergence rates of the SWR method using the theoretically analyzed parameter and the one from the numerical optimization are very similar in the two subdomain case, while for the eight-subdomain case, these two Robin parameters lead to slightly different outcomes. This is natural since the Robin parameter is analyzed in the two-subdomain case, and in the many subdomain case, an optimal choice of should depend on the number of subdomains (see, e.g., [48]).
(a)
(b)
6. Conclusions
We studied the Schwarz waveform relaxation (SWR) algorithms with Robin transmission conditions (TCs) for a class of representative time-fractional PDEs. The problem of determining the parameter for the Robin TCs was investigated at the continuous and discrete level. Particularly, at the continuous level, the parameter obtained by solving a transformed min-max problem (i.e., (36a) and (36b)) is found much better than the one predicted by the equioscillation principle [30,35]. The solution of the transformed min-max problem is not always determined by the equioscillation principle.
At the discrete level, we used the centered finite difference scheme to discretize the spatial variable and then we used the convolution quadratures [38–40], based on linear multistep methods, to get time discretizations. This framework permits us to treat the time discretizations as a whole; that is, the discrete Laplace transform of the sequences generated by the quadratures has a uniform expression, which is mainly determined by the quotient associated with a specified time-integrator. Three time-integrators were carefully studied (i.e., backward Euler, trapezoidal, and BDF2), and it was shown that the trapezoidal method leads to a convergence factor equaling 1 (independent of the Robin parameter) and that for the backward Euler and BDF2 methods, the convergence factors are (almost) equal and both are satisfactory.
Based on the work in this paper, our ongoing research is twofold. First, we attempt to generalize the two-subdomain analysis to the multisubdomain case; this is a necessary step to make the SWR algorithm practical for real parallel computation. Second, we try to solve the general min-max problem (97) (i.e., large overlap size with ) and perform an asymptotic analysis for the convergence factor with respect to the discretization parameters (and ) and , the length of the time-interval.
Appendix.
A. The Proof of Lemma 4
Proof. For the function defined by (36b), we haveHence, for , the solution of (36a) satisfies , since otherwise, if (resp. ), we can further uniformly reduce by increasing (resp. decreasing) the value of . Hence, it suffices to studywhere and .
A routine calculation shows that, for a given , the function has at most two local extremums, located atwhere is the discriminant. Since as , the larger one must be a local maximizer. This local maximizer does not always exist, and by letting , we know that it exists for , where . By using and , we know that for , i.e., , the function is strictly decreasing for .
Case A. If . In this case, we have for . Hence, , with . From (A.1), it is easy to know that the function (resp. is an increasing (resp. decreasing) function of when . Then, it can be shown that for the function attains its global minimum at , where is given by (37b); this implies . We havewhere for any real numbers and , the set (or ) is empty if .
Case B. . In this case, it is clear that . We claim for . Note that, for any , we have Hence, it is sufficient to prove . Since it is sufficient to prove . By noticing we finally get for . Now, by using and the continuity of with respect to , it is easy to getWe next claim for . From (A.1), we havewhere we used the fact that is the local maximizer of , that is . Then, it holds thatwhere the last inequality holds for any , becauseTherefore, and this implies that the function is decreasing for . Based on the above analysis, we next consider the following two cases:(1), , and . In this case, it must hold , since otherwise we have and this leads to contradiction. Hence, it holds that for . Then, by using , , and (A.8), we know that the equation has a unique root and that attains its global minimum at , because and ) are, respectively, decreasing and increasing functions of .(2)If any one of the conditions in 1 is not satisfied, we claimIndeed, this can be proved by using (A.8) and the fact that “both and are decreasing functions, but is an increasing function.” For example, if , by using (A.8), we know that (A.12) holds. For the case or , (A.12) can be proved similarly. Now, let and . Then, we haveThis, together with (A.12), implies that the function can reach the lower bound ; hence, the conclusion follows.
Case C. . In this case, we have and therefore the local maximizer exists for all . The first case in (39a) and (39b) has been analyzed in Case B. For , , and , the analysis is also similar to Case B. If , , , and , we haveThis implies that the equation has a unique solution , which is the global minimizer, because and are, respectively, decreasing and increasing functions for . For the last case, if , by using the fact that “both and are decreasing functions,” it is easy to get ; otherwise, it shall hold that for all and because of the decreasing property of , we also have .
Data Availability
The data used to support the findings of this study are included within the manuscript.
Conflicts of Interest
It is declared that there are no conflicts of interest among the authors regarding this manuscript.
Acknowledgments
This paper was supported financially by the National Science and Technology Major Project of China (no. 2016ZX05065 and no. 2016ZX05042-003), by the Sichuan Science and Technology Program (2019YJ0541), and by the Opening Project of Sichuan Province University Key Laboratory of Bridge Non-Destruction Detecting and Engineering Computing (2019QZJ03).