Convergence Analysis of Schwarz Waveform Relaxation for Nonlocal Diffusion Problems

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.


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 as where c is the anomalous diffusion exponent. e process is usually referred to as subdiffusive when c ∈ (0, 1). For the ordinary (or Brownian) diffusion (this corresponds to c � 1), from a continuous point of view the diffusion process is described by the well-known diffusion equation u t (x, t) � u xx (x, t), where u(x, t) represents the probability density of finding a "particle" at space x and at time t. It turns out that the probability density function u(x, t) that describes anomalous (sub) diffusive particles follows the fractional diffusion equation [1][2][3][4]: where D 1− c t 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][9][10][11][12][13][14][15][16][17]. Very recently, the so-called fractional cable equation was introduced for modeling the anomalous diffusion in spiny neuronal dendrites [18][19][20]: where α, β ∈ (0, 1). is 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), d denotes the diameter of the cable, E is the equilibrium potential of the membrane, R m is a membrane constant and C m is a constant independent of the membrane area, and R A is an axial material constant, associated to the given cable. By taking τ m � R m C m , ] � (dR m /4R A ) and then by letting , cable equation (3) can be transformed to z t u � D β t z 2 x u − μD α t u + f with μ � τ α− 1 m . 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. ese 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. is 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). ere are also other efficient algorithms for handling the fractional PDEs, such as the interesting parallel-in-time algorithms studied in [25][26][27][28][29].
For integer order differential equations, the SWR algorithms have been investigated widely and deeply by Gander and his colleagues [30][31][32][33][34][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. e 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][36][37] and the main finding is that a good parameter can be computed by solving one or two simple nonlinear equations.
is 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 timefractional 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][39][40].
is 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.

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 μ ≥ 0 and α, β ∈ (0, 1). is 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 c ∈ (0, 1), we denote by D where Γ is the gamma function Γ(τ) � +∞ 0 e − t t τ− 1 dt. From [41] (Chapter 2), for any v ∈ C 1 (0, T) it holds that lim c⟶1 − D c t v � v ′ (t)(∀t ∈ (0, T)). Hence, the Riemann-Liouville fractional derivative builds a "bridge" from v(t) to v ′ (t).
In the field of fractional PDEs, there exists another definition of the fractional derivative, the Caputo derivative: e two definitions are linked by the following relation [41], (pp. 62-77): anks to this connection, we shall only consider the Riemann-Liouville derivatives, since the results can be directly generalized to problems with Caputo derivatives.
ere exists another class of fractional PDEs, which are also widely studied in the literature [41,42], the fractional heat equations:

Mathematical Problems in Engineering
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, Ω � Ω 1 ∪ Ω 2 , with Ω 1 � (− ∞, L] and Ω 2 � [0, +∞). Here, L > 0 denotes the size of overlap. en, the SWR algorithms for (5) are where p ∈ R is a free parameter, k ≥ 1 is iteration index, and j � 1, 2. 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.

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 v(t): Applying Laplace transform to the error equation (12), we get e solutions can be expressed as e k j � A k j e λ(s)(x− L) + B k j e − λ(s)x . en, using the fact that e k 1,2 do not increase exponentially at infinity, we get Substituting (14) into the transmission conditions in (13), we get and this gives A k+2 Hence, e k j (x, s) � ρ c (p, L, s)e k− 2 j (x, s), ∀x ∈ Ω j . Mathematically, we want |ρ c | ≪ 1, which leads to the following min-max problem: min Lemma 1. Assume α, β ∈ (0, 1), μ ≥ 0 and p > 0. en, the argument ρ c defined by (16) is analytic in the right half of the complex plane.
Proof. Let s � re iφ with r > 0 and φ ∈ (− (π/2), (π/2)). en, we have is implies that, for μ ≥ 0 and α, β ∈ (0, 1) the function λ(s) 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 p > 0, the denominator λ(s) + p does not vanish in the right half of the complex plane; hence, (p − λ(s)/p + λ(s)) is an analytic function. Since the product of two analytic functions is still analytic, ρ c is analytic. □ Note that, for any s � re iφ with φ ∈ (− (π/2), (π/2)), the limit of |ρ c | as r ⟶ + ∞ is zero. Hence, by using the maximum principle for complex analytic functions, we know that the maximum of |ρ c (p, L, s)| is obtained on the imaginary axis, i.e., s � iω. 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 Δt and a length of time window (0, T), as |ω| ∈ [(π/T), (π/Δt)]. Hence, it is reasonable to define ρ c (p, L), which we call convergence factor of the SWR algorithm at the space-time continuous level, as

Mathematical Problems in Engineering
e subscript "c" appearing in ρ c and ρ c denotes "continuous." Clearly, an ideal p is the minimizer of the convergence factor ρ c (p, L), i.e., min p∈R max For any ω > 0, it holds that us, for any ω > 0, we can rewrite λ( ± iω) as is gives Clearly, for any p < 0, it holds |ρ c | > 1. us, min-max problem (20) is equivalent to where

Towards a Good Choice of q(� 2Lp): e Equioscillation
Principle. A well-understood formula for choosing the free parameter p is based on the so-called equioscillation principle, established recently by Gander and Halpern [30,35]: where ω min � (π/T) and ω(q) denotes the local maximizer of R. e solution of (25), namely, q equi , can be obtained by solving the following two nonlinear equations: R ω min , q � R(ω, q), together with a constraint z 2 ω R(ω, q) < 0 (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 advectionreaction-diffusion equations, formula (25) gives the best choice of q; see [30,35]. In Figure 1(a), we have given an illustration of the equioscillation principle for the min-max problem (24a)-(24b).
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 [(π/T), (π/Δt)]? For integer order time-dependent PDEs, the answers for these questions are positive, at least in the asymptotic sense, i.e., L � O(Δx), Δt � O(Δx r ), and Δx 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 4 Mathematical Problems in Engineering Figure 1(a), the parameter q equi results in a convergence factor around 0.842, which implies that the SWR algorithm using p � (q equi /2L) 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 q. In Figure 1(b), with Δt � 0.001, T � 50, L � 0.02, μ � 1.3, and two groups of α and β, we plot the function R(ω, q) for ω ∈ [(π/T), (π/Δt)], where we can see that compared to q equi the new parameter q opt leads to much smaller convergence factor and that with the new parameter the function R(ω, q) does not equioscillate.

Towards a Good Choice of q(� 2Lp): A New Formula.
Because of the high complexity of the arguments λ I,R (ω), it is difficult to present a closed formula for the solution of minmax problem (24a) and (24b). In what follows, we focus on deriving a reliable approximate solution. We need two auxiliary lemmas. Lemma 2. Let α, β ∈ (0, 1), μ ≥ 0 and ω min = π/T ω max = π/∆t (c) Figure 1: (a) An illustration of equioscillation, i.e., the parameter q � q equi is obtained by solving (26). e plot of R(ω, q) in the relevant frequency domain ω ∈ [(π/T), (π/Δt)] for two choices of q, (b) q � q equi and (c) q � q opt given by Lemma 3. en, it holds that max x∈[x 0 , Proof. Clearly, we only need to consider the case μ > 0, since for μ � 0 the function S(x) is a constant and therefore For any x, it holds that and this, together with μ > 0 and α, β ∈ (0, 1), gives Hence, if x * corresponds to a local extremum of S(x), it must be a minimizer. Hence, By noticing χ 2 Hence, from Lemma 2, we get □ Now, for min-max problem (24a) and (24b), by using Lemma 3, we have For any p > 0, this leads to the following inequality: In (34), we have used p > 0 and λ R (ω) > 0. Let λ R (ω) � y/(2L) and p � q/(2L). en, 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 ϕ ∈ (0, 1) and y 1 > y 0 > 0. e concrete proof is presented in Appendix.

Theorem 1. Let S be the function defined in Lemma 2 and
en, for specified α, β ∈ (0, 1), μ ≥ 0, and T, Δt > 0, the parameter p involved in the Robin TCs can be chosen as p * c � q opt /(2L), where q opt , given by Lemma 4, is the solution of min-max problem (36a) and (36b) and L > 0 denotes the size of overlap. e quantities y 0 and y 1 in (36b) can be explicitly written down, as With p * c , the convergence factor ρ c defined by (19) can be bounded by
We now consider the case μ > 0. Suppose the function S(x) has a local extremum located at x � x * . en, we have (42), it is easy to get S 0 ′ (x * ) > 0, and this, together with α, β ∈ (0, 1), implies is implies that S(x) has no local maximum(s) and thus y 0,1 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. en, 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 z t u � z xx u − μu + f (or more generally, the advection-reaction-diffusion where y 1 > y 0 > 0 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 realvalued functions, and roots of quadratic/cubic polynomials), the concrete procedures are different. is difference can be seen from two points:

Mathematical Problems in Engineering
(1) e two functions, R Int− order in (40) and R in (36a) and (36b), are different (2)Because of this difference, the solution q opt of minmax problem (36a) and (36b) is not always determined by the equioscillation principle For min-max problem (36a) and (36b), there are two reasons that the equioscillation principle does not hold. First, the local maximizer y + (q opt ) does not exist in the interval [y 0 , y 1 ] (see Figure 2(a)). Second, more typically, even though the local maximizer y + (q opt ) exists, R(y 0 , q opt ) ≠ R(y + (q opt ), q opt ) (see Figure 2(b)). What we observed in these two subfigures does not happen occasionally, and it seems an inherent feature for the timefractional problems. Precisely, as will be shown in our forthcoming paper, in the asymptotic sense, i.e., L � O(Δx), Δt � O(Δx r ) and Δx small, the equioscillation property always does not hold for min-max problem (36a) and (36b) when r < (1/1 − β) and α, β ∈ (0, 1).

Remark 2.
About the parameter q � q opt : the solution q opt of min-max problem (36a) and (36b) depends on y 0,1 and ϕ, which are related to α, β, μ, T (problem parameters) and Δt, L (discretization parameters). Hence, we denote q opt by which we call "parameter formula" and its concrete expression is given by Lemma 5 as we have seen. Note that q opt 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 q opt is really a good approximate solution. Our evidence is based on numerical optimization: (1) Find the best parameter q best through repetitive numerical optimizations for original min-max problem (24a) and (24b): (ii) Compare the values R j J j�0 , and the minimal one corresponds to q best .
(2) Compute and then compare the values of R max (q) ≔ max ω∈[(π/T),(π/Δt)] R(w, q), i.e., the convergence factor of the SWR algorithm with q � q opt (the approximate solution) and q � q best obtained through numerical optimizations.
Let Δt � 0.02, T � 40, L � 0.04, and μ � 8. en, in Figure 3, we plot R max (q opt ) (Figure 3(a)) and R max (q best ) (Figure 3(b)) for α, β ∈ [0, 1], where a small step-size, Δq � 10 − 3 , is used in the above numerical optimization procedure to get q best . In Figure 3(c), we plot R max (q opt ) − R max (q best ), i.e., the error between the convergence factors using the two parameters. is figure clearly shows that the quantity q opt , given by Lemma 4, is really a good approximate solution of original min-max problem (24a) and (24b).

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 t, we get

Derivation of the Discrete SWR Algorithms.
To get discrete algorithm, we apply the centered finite difference formula to (49), which gives To simplify the convergence analysis, we assume L � Δx, i.e., the size of overlap equals a single space mesh size. en, the artificial boundary interfaces x � L (resp. x � 0) corresponds to x 1 (resp. x 0 ) and the Robin TCs are discretized as Local maximizer ȳ + (q opt ) exists, but R (y 0 , q opt ) ≠ R (ȳ + (q opt , y opt ) q opt = 0.80286 by new formula Here, we discretized the spatial derivatives in the Robin TCs as , (51c) 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, as where c � 1 − α (or 1 − β) and V is not known beforehand but are computed as the computing proceeds from time step to time step. A naive implementation would require O(N 2 t ) operations and O(N t ) memory for per expansion coefficients over N t time steps. In [44], fast convolution algorithms are constructed, which implements, over N t time steps, with O(N t log(N t )) operations and O(log(N t )) active memory.
e starting point of the fast convolution algorithms is to replace the kernel function K c (t − τ) in (52) by its inverse Laplace transform, along some contour Υ: where K c (s) � L(K c ) � s − c denotes the Laplace forward transform of the kernel function K c . e 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 e ts decays rapidly for growing |s| along Υ. In (53), reversing the order of integration gives e function y(0, t; s) is the solution of the following linear scale ODE at τ � t: Since V(t) is only available at discrete grids t n N t n�0 , it is reasonable to apply some time-integrator, e.g., the linear multistep method, to (55) to get y at t � t n : Multiplying (56) by ζ n and summing over n from 0 to +∞, we obtain with the generating (formal) power series y(ζ) � +∞ n�0 y n ζ n and V(ζ) � +∞ n�0 V n ζ n . Solving y(ζ) from (57), we get where Λ(ζ) is called the quotient of the generating polynomial of a linear multistep method. Clearly, the solution y n of (56) is the n-th coefficient of the power series . By substituting this into the right-hand side of (54), we get Now, we expand K c (Λ(ζ)/Δt) in the formal power series, as Suppose the time-integrator applied to (55) is an A-stable method with order p. en, from [40], eorem 2.1, we know that (61) is an approximation to (52) and  For general case, we have to rely on numerical calculation to get c j (c, Δt) +∞ j�0 (see [44] for more details). For some special cases, such as the linear θ-method (Λ(ζ) � (1 − ζ)/(θ + (1 − θ)ζ)) and the two-step backward difference formula (BDF2, Λ(ζ) � 2 j�1 (1 − ζ) j /j), these weights can be written down, as With the quadrature weights c n (c, Δt) +∞ n�0 in our hands, we denote the discrete convolution (61) by (K c,Δt V)(n), that is, en, the scheme of the fully discrete SWR algorithms on the j-th subdomain is Together with boundary conditions (51b) at t � t n .
Proof. With the arguments σ 0 , σ 1 , and ω max given by (81b), en, the proof is divided into two steps: a(s) ∈ D if s ∈ D 0 and s + ����� e function a(s) is defined by (68).

(97)
For M ≥ 2 and η > 1, the function T(s, η) has local maximums in the interior of D, and therefore, the analysis is much more complicated compared to the case M � 1. e proof of eorem 2 cannot be straightforwardly generalized to (97), and we will investigate this general min-max problem in other places.

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 where x(t) � (· · · , x − 1 (t), x 0 (t), x 1 (t), · · ·) T denotes the vector of the nodal voltage values, c � (1/ ��� RC √ ) with C the value of the fractional capacitor and R the value of the resistance, as shown in Figure 5(a). e fractional derivative in (98) is defined in the Caputo sense (see (8)), and the argument b ∈ (0, 1] is determined by the branches of the circuit shown on the right. Mathematically, (98) can be equivalently rewritten as the fractional integral defined by (50) and some other function f. Clearly, this is a special case of (), which can be seen by letting 1 − β � b and μ � 0 and then by replacing U k j,m (t) by x m (t) and δ 2 x by A there. erefore, 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 with instead of s � Λ(e − sΔt )/Δt as shown in (78b). Hence, the rectangular region D covering ξ + (s) is a semi-infinite horizontal stripe, i.e., σ max � +∞ in (80). In this limit case, eorem 2 is still applicable, but we need to replace η max by +∞ there and redefine T ± (η) as where η + and η min are the same quantities as we have stated in eorem 2.

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 L � Δx. Since the parameter p appearing in the continuous Robin TCs is transformed to η � 1 + Δxp in the discrete circumstance, we shall only mention η throughout this section.
At the moment, for given problem/discretization parameters, we have two choices of η, η � 1 + (Δxq opt /(2L)) � 1 + (q opt /2) from the continuous analysis (i.e., Lemma 4) and η � η opt from the discrete analysis (i.e., eorem 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 ρ d ≔ max s∈H |ρ d | based on three choices of η, where ρ d and H are given in (78a) and (78b). Here, the problem/discretization parameters are chosen as μ � 4.8, T � 18, Δt � 0.025, and Δx � 0.04. To get η opt from eorem 2, the boundaries of the rectangular region D, i.e., σ min , σ max , and ω max , 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).
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. ese 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 α � 0.84, β � 0.5, and μ � 4.8 for the problem parameters, and Δt � 0.025 and Δx � 0.04 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 η � 1.12031 from the continuous analysis. (For the parameters obtained from the discrete analysis, this error plot looks similar.) e 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. e 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 max s∈H |ρ d | � 1, independent of η.
en, 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. e 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.
At the end of this section, we show in Figure 8 the effects of the discretization parameters and the length of timeinterval 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 Δt and T, we use two choices for the parameter η involved in the Robin TCs, the one η � 1 + Δxq opt from the continuous analysis and the one η � η best from the numerical optimization. (For the parameter η � η opt from the discrete analysis, the plots look similar.) For each Δt and T, the iteration number is measured when the error between the iterate and the reference solution is less than 10 − 8 . e 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. is is natural since the Robin parameter η � 1 + Δxq opt 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]).

Conclusions
We studied the Schwarz waveform relaxation (SWR) algorithms with Robin transmission conditions (TCs) for a class of representative time-fractional PDEs. e 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]. e 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][39][40], based on linear multistep methods, to get time discretizations.
is 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. ree 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 L � MΔx with M ≥ 2) and perform an asymptotic analysis for the convergence factor with respect to the discretization parameters Δt (and Δx) and T, the length of the timeinterval.
(2) If any one of the conditions in 1 is not satisfied, we claim max y∈ y 0 ,y 1 Indeed, this can be proved by using (A.8) and the fact that "both R(q) and R(y 1 , q) are decreasing functions, but R(y 0 , q) is an increasing function." For example, if q † � q d , by using (A.8), we know that (A.12) holds. For the case q † > q d or y(q † ) ≥ y 1 , (A.12) can be proved similarly. Now, let Q � q: q ∈ [q 0 , q d ], and y + (q) < y 1 and Q 0 � [q 0 , q 1 ]/Q. en, we have is, together with (A.12), implies that the function R max (q)(� max y∈[y 0 ,y 1 ] R(y, q)) can reach the lower bound R † (q † ); hence, the conclusion q opt � q † follows.

Data Availability
e 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.