# Convergence Analysis of Schwarz Waveform Relaxation for Nonlocal Diffusion Problems

Authors: Ke Li; Dali Guo; Yunxiang Zhao
Journal: Mathematical Problems in Engineering
 (2021)

---

## 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.


---

## Body

## 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 as(1)x2t∼2Γ1+γtγ,t⟶+∞,where γ is the anomalous diffusion exponent. The process is usually referred to as subdiffusive when γ∈0,1. For the ordinary (or Brownian) diffusion (this corresponds to γ=1), from a continuous point of view the diffusion process is described by the well-known diffusion equation utx,t=uxxx,t, where ux,t represents the probability density of finding a “particle” at space x and at time t. It turns out that the probability density function ux,t that describes anomalous (sub) diffusive particles follows the fractional diffusion equation [1–4]:(2)∂ux,t∂t=Dt1−γ∂2ux,t∂x2+f,where Dt1−γ 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]:(3)Cm∂tux,t=Dtβd4RA∂x2ux,t−Dtαux,t+ERm,where α,β∈0,1. 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), d denotes the diameter of the cable, E is the equilibrium potential of the membrane, Rm is a membrane constant and Cm is a constant independent of the membrane area, and RA is an axial material constant, associated to the given cable. By taking τm=RmCm, ν=dRm/4RA and then by letting f=−DtαE/Rm,t˜=t/τm,x˜=xτm0.51−β/ν, cable equation (3) can be transformed to ∂t˜u=Dt˜β∂x˜2u−μDt˜αu+f with μ=τmα−1. Dropping the tilde superscript, we finally get(4)∂tu=Dtβ∂x2u−μDtαu+f.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 Section3, 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:(5)∂tu−Dtβ∂xxu+μDtαu=fx,t,x,t∈ℝ×0,T,ux,0=u0,x∈ℝ,where μ≥0 and α,β∈0,1. 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 γ∈0,1, we denote by Dtγ the Riemann–Liouville fractional derivative:(6)Dtγvt=1Γ1−γddt∫0t1t−τγvτdτ,where Γ is the gamma function Γτ=∫0+∞e−ttτ−1dt. From [41] (Chapter 2), for any v∈ℂ10,T it holds that limγ⟶1−Dtγv=v′t∀t∈0,T. Hence, the Riemann–Liouville fractional derivative builds a “bridge” from vt to v′t.In the field of fractional PDEs, there exists another definition of the fractional derivative, the Caputo derivative:(7)cDtγvt=1Γ1−γ∫0t1t−τγv′τdτ.The two definitions are linked by the following relation [41], (pp. 62–77):(8)Dtγvt=1Γ1−γ∫0tv′τt−τγdτ︸Caputo derivativecDtγvt+v0t−γΓ1−γ,∀γ∈0,1.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:(9a)Dtαu=∂xxu+f.From [41] (pp. 65–75, Chapter 2), we know that (9a) can be rewritten as(9b)∂tu=Dt1−α∂xxu+f˜,f˜depends onfandu0.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. Then, the SWR algorithms for (5) are(10)∂tujk−Dtβ∂xxujk+μDtαujk=fx,t,x∈Ωj,∂x+−1j+1pujk2−jL,t=∂x+−1j+1pu3−jk−12−jL,t,ujkx,0=u0,where p∈ℝ 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.


## 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 vt:ℝ+⟶ℝ, we denote by v^s the Laplace transform of v, defined by(11)v^s≔ℒv=∫0+∞e−stvtdt,s=σ+iω,σ>0.Letejkx,t=ujkx,t−ux,t|x∈Ωj. Then, it is easy to know that ejk satisfies(12)∂tejk−Dtβ∂xxejk+μDtαejk=0,∂x+−1j+1pejk2−jL,t=∂x+−1j+1pe3−jk−12−jL,t,ejkx,0=0.Applying Laplace transform to the error equation (12), we get(13)∂x2e^jk=λ2se^jk,λs=s1−β+μsα−βwith non−negative real part,∂x+−1j+1pe^jk2−jL,t=∂x+−1j+1pe^3−jk−12−jL,t.The solutions can be expressed ase^jk=Ajkeλsx−L+Bjke−λsx. Then, using the fact that e^1,2k do not increase exponentially at infinity, we get(14)e^1k=A1keλsx−L,e^2k=B2ke−λsx.Substituting (14) into the transmission conditions in (13), we get(15)p+λsA1k=p−λse−λsLB2k−1,p+λsB2k=p−λse−λsLA1k−1,and this gives A1k+2=ρ^cp,L,sA1k and B2k+2=ρ^cp,L,sB2k, where(16)ρ^cp,L,s=p−λsp+λs2e−2Lλs.Hence,e^jkx,s=ρ^cp,L,se^jk−2x,s,∀x∈Ωj. Mathematically, we want ρ^c≪1, which leads to the following min-max problem:(17)minp∈ℝmaxℜs≥0ρ^cp,L,s.Lemma 1.
Assumeα,β∈0,1, μ≥0 and p>0. Then, the argument ρ^c defined by (16) is analytic in the right half of the complex plane.Proof.
Lets=reiφ with r>0 and φ∈−π/2,π/2. Then, we have(18)ℜs1−β+μsα−β=r1−βcos1−βφ+μrα−βcosα−βφ.
This 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 anys=reiφ 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 ρ^cp,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 ρcp,L, which we call convergence factor of the SWR algorithm at the space-time continuous level, as(19)ρcp,L=maxω∈π/T,π/Δtρ^cp,L,iω2.The subscript “c” appearing in ρ^c and ρc denotes “continuous.” Clearly, an ideal p is the minimizer of the convergence factor ρcp,L, i.e.,(20)minp∈ℝmaxω∈π/T,π/Δtp−λiωp+λiω2e−2Lℜλiω.For anyω>0, it holds that(21)i±ω+μi±ωαi±ωβ=±i1−βω+μ±iα−βωαωβ=ωcos±1−β/2π+isin±1−β/2π+μωαcos±α−β/2π+isin±α−β/2πωβ=ω1−βcos1−β2π+μωα−βcosα−β2π±iω1−βsin1−β2π+μωα−βsinα−β2π.Thus, for anyω>0, we can rewrite λ±iω as λ±ω=λRω±iλIω, where(22)λRω=χR2ω+χI2ω+χRω2,λIω=χR2ω+χI2ω−χRω2,χRω=ω1−βcos1−β2π+μωα−βcosα−β2π,χIω=ω1−βsin1−β2π+μωα−βsinα−β2π.This gives(23)ρ^cp,L,±iω=λRω−p2+λI2ωλRω+p2+λI2ωe−2LλRω,∀ω∈πT,πΔt.Clearly, for anyp<0, it holds ρ^c>1. Thus, min-max problem (20) is equivalent to(24a)minq>0maxω∈π/T,π/Δtℛ˜ω,q,where(24b)q=2Lp,ℛ˜ω,q=λ˜Rω−q2+λ˜I2ωλ˜Rω+q2+λ˜I2ωe−λ˜Rω,λ˜R,Iω=2LλR,Iω.


### 3.1. Towards a Good Choice ofq=2Lp: The Equioscillation Principle

A well-understood formula for choosing the free parameterp is based on the so-called equioscillation principle, established recently by Gander and Halpern [30,35]:(25)ℛ˜ωmin,q=ℛ˜ω¯q,q,where ωmin=π/T and ω¯q denotes the local maximizer of ℛ˜. The solution of (25), namely, qequi, can be obtained by solving the following two nonlinear equations:(26)ℛ˜ωmin,q=ℛ˜ω,q,∂ωℛ˜ω,q=0,together with a constraint ∂ω2ℛ˜ω,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 advection-reaction-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).Figure 1
(a) An illustration of equioscillation, i.e., the parameterq=qequi is obtained by solving (26). The plot of ℛ˜ω,q in the relevant frequency domain ω∈π/T,π/Δt for two choices of q, (b) q=qequi and (c) q=qopt given by Lemma 3.
(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 π/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Δxr, 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 Figure1(a), the parameter qequi results in a convergence factor around 0.842, which implies that the SWR algorithm using p=qequi/2L for the Robin parameter converges slowly.Since the applicability of the equioscillation principle to time-fractional problems is unclear, we introduce in Lemma4 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 ℛ˜ω,q for ω∈π/T,π/Δt, where we can see that compared to qequi the new parameter qopt leads to much smaller convergence factor and that with the new parameter the function ℛ˜ω,q does not equioscillate.


### 3.2. Towards a Good Choice ofq=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 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α,β∈0,1, μ≥0 and(27)Sx=x2+2μxcos1−α/2π+μ2−xcos1−β/2π+μcosα−β/2πx2+2μxcos1−α/2π+μ2+xcos1−β/2π+μcosα−β/2π.Then, it holds thatmaxx∈x0,x1Sx=maxSx0,Sx1, where 0<x0<x1.Proof.
Clearly, we only need to consider the caseμ>0, since for μ=0 the function Sx is a constant and therefore maxx∈x0,x1Sx=maxSx0,Sx1.
LetR0x=x2+2μxcos1−α/2π+μ2 and R1x=xcos1−β/2π+μcosα−β/2π. We then have Sx=R0x−R1x/R0x+R1x and S′x=R1xR0′x−2R0xR1′x/R0xR0x+R1x2. Suppose S′x has a positive root located at x=x∗, i.e., R1x∗R0′x−2R0x∗R1′x∗=0. Then, we have(28)S″x∗=R1xR0″x−2R0xR1″x−R1′xR0′xR0xR0x+R1x2|x=x∗.
For anyx, it holds that(29)R0′x=2x+2μcos1−α2π,R1′x=cos1−β2π,R0″x=2,R1″x=0,and this, together with μ>0 and α,β∈0,1, gives(30)S″x∗=2μcosα−β/2π−cos1−β/2πcos1−α/2πR0x∗R0x∗+R1x∗2=2μsin1−β/2πsin1−α/2πR0x∗R0x∗+R1x∗2>0.
Hence, ifx∗ corresponds to a local extremum of Sx, it must be a minimizer.Lemma 3.
Letμ≥0, α,β∈0,1, and ϕ=maxSπ/T1−α,Sπ/Δt1−α with S being the function defined in Lemma 2. Then, for all ω∈π/T,π/Δt, we have λI2ω≤ϕλR2ω.Proof.
From (22), χIω/χRω=ω1−αsin1−β/2π+μsinα−β/2π/ω1−αcos1−β/2π+μcosα−β/2π. Since α,β∈0,1,(31)cos1−β2πsinα−β2π≤sin1−β2πcosα−β2π.
Hence,ω1−αcos1−β/2π+μcosα−β/2π/ω1−αsin1−β/2π+μsinα−β/2π is a decreasing function of ω1−α; this gives(32)χIπ/ΔtχRπ/Δt≤χIωχRω≤χIπ/TχRπ/T,∀ω∈πT,πΔt.
By noticingχR2ω+χI2ω=ω21−β+μ2ω2α−β+2μω1+α−2βcos1−α/2π, we have(33)λI2ωλR2ω≤χR2ω+χI2ω−χRωχR2ω+χI2ω+χRω=Sω1−α.
Hence, from Lemma2, we get λI2ω/λR2ω≤maxSπ/T1−α,Sπ/Δt1−α=ϕ and this completes the proof.Now, for min-max problem (24a) and (24b), by using Lemma 3, we have(34)λRω−p2+λI2ωλRω+p2+λI2ωe−2LλRω≤λRω−p2+ϕλR2ωλRω+p2+ϕλR2ωe−2LλRω.For anyp>0, this leads to the following inequality:(35)maxω∈π/T,π/ΔtλRω−p2+λI2ωλRω+p2+λI2ωe−2LλRω≤maxω∈π/T,π/ΔtλRω−p2+ϕλR2ωλRω+p2+ϕλR2ωe−2LλRω.In (34), we have used p>0 and λRω>0. Let λRω=y/2L and p=q/2L. Then, the solution of the following min-max problem can be regarded as a reliable approximate solution of (24a) and (24b):(36a)minq>0maxy∈y0,y1ℛy,q,withℛy,q≔y−q2+ϕy2y+q2+ϕy2e−y,where(36b)y0≔2Lminω∈π/T,π/ΔtλRω,y1≔2Lmaxω∈π/T,π/ΔtλRω.We present the solution of (36a) and (36b) in a very general setting: we only assume ϕ∈0,1 and y1>y0>0. The concrete proof is presented in Appendix.Lemma 4.
Parameter optimization at the continuous level: assumeϕ∈0,1 and y1>y0>0. Let q0=y01+ϕ, q1=y11+ϕ, ℛmaxq=maxy∈y0,y1ℛy,q, and qopt be the solution of min-max problem (36a) and (36b). Then,(a)
Ify0≥1+ϕ1+1/ϕ−1, we have(37a)qopt=q†,ℛmaxqopt=ℛ†qopt,whereℛ†q=maxℛy0,q,ℛy1,q and q† is defined by(37b)q†=q0,ifℛy0,q0≥ℛy1,q0,q1,ifℛy1,q1≥ℛy0,q1,q˜†,otherwise.The argumentq˜†∈q0,q1 is the unique root of ℛy0,q=ℛy1,q.(b)
Ify0<1+ϕ1+1/ϕ−1<y1, we have(38a)qopt=q0∗,ifqd>q†,y¯+q†<y1andℛ¯q†>ℛ†q†,q†,otherwise,ℛmaxqopt=ℛ¯qopt,ifqd>q†,y¯+q†<y1,andℛ¯q†>ℛ†q†,ℛ†qopt,otherwise,whereqd=1+ϕ1+1/ϕ−1, y¯+q=1−ϕq2+21+ϕq+2qdq/1+ϕ, ℛ¯q=ℛy¯+q,q, dq=−ϕq2−2ϕ1+ϕq+1+ϕ2, and q0∗∈q†,qd is the unique solution of(38b)ℛ¯q=ℛy0,q.(c)
Ify1≤1+ϕ1+1/ϕ−1, then(39a)qopt=q†,if,y¯+q†≥y1,ory¯+q†<y1,andℛ¯q†≤ℛ†q†,q0∗,if,y¯+q†<y1,ℛ¯q†>ℛ†q†,andy¯+q1≥y1,or,y¯+q1<y1,andℛ¯q1<ℛ†q1,q1,if,maxy¯+q†,y¯+q1<y1,ℛ¯q†≥ℛ†q†,andℛ¯q1≥ℛ†q1,(39b)ℛmaxqopt=ℛ†qopt,if,y¯+q†≥y1,ory¯+q†<y1,andℛ¯q†≤ℛ†q†,ℛ¯qopt,otherwise,whereq0∗∈q†,q1 is the unique solution of the equation (38b).Theorem 1.
LetS be the function defined in Lemma 2 and ϕ=maxSπ/T1−α,Sπ/Δt1−α.Then, for specifiedα,β∈0,1, μ≥0, and T,Δt>0, the parameter p involved in the Robin TCs can be chosen as pc∗=qopt/2L, where qopt, given by Lemma 4, is the solution of min-max problem (36a) and (36b) and L>0 denotes the size of overlap. The quantities y0 and y1 in (36b) can be explicitly written down, as(40)y0=2LλRπT,if,α≥β,2LminλRπT,λRπΔt,if,α<β,y1=2LλRπΔt,if,α≥β,2LmaxλRπT,λRπΔt,if,α<β.Withpc∗, the convergence factor ρc defined by (19) can be bounded by(41)ρcpc∗,L≤ℛmaxqopt,where ℛmaxqopt is given by (37a), or (38a), or (39b).Proof.
bound (41) can be derived from (34) and (36b):
(34), (36b)⇒ρcp,L≤ℛmaxq≔maxy∈y0,y1ℛy,q,∀q=2Lp>0.
Hence, it remains to prove that the argumentsy0 and y1 defined by (36b) can be explicitly expressed as (40). To this end, we define(42)S0x=x2+2μxcos1−α2π+μ2+xcos1−β2π+μcosα−β2π.
Clearly, for anyα,β∈0,1 and μ≥0, the function S0x is increasing for x>0.Case 1.
α≥β. In this case, by noticing(43)2λR2ω=ωα−βS0ω1−α,we know that the function λRω is increasing for ω∈π/T,π/Δt. Hence, we have(44)y0=2Lminω∈π/T,π/ΔtλRω=L2πTα−βS0πT1−α,y1=2Lmaxω∈π/T,π/ΔtλRω=L2πΔtα−βS0πΔt1−α.Case 2.
α<β. In this case, from (43), we get 2λR2ω=S˜x≔S0xx−β−α/1−α, with x=ω1−α. We have S˜′x=xS0′x−β−α/1−αS0xx−1+β−2α/1−α and for μ=0, it holds that S˜′x=1−β/1−αS0xx−1+β−2α/1−α>0 if x>0 and α,β∈0,1; therefore, y0 and y1 satisfy (40).
We now consider the caseμ>0. Suppose the function S˜x has a local extremum located at x=x∗. Then, we have x∗S0′x∗=β−α/1−αS0x∗. Hence, it holds that S˜″x∗=x∗S0″x∗+1−β/1−αS0′x∗x∗−1+β−2α/1−α. We claim S˜″x∗>0. Indeed, from (42), it is easy to get S0′x∗>0, and this, together with α,β∈0,1, implies 1−β/1−αS0′x∗>0. For S0″x∗, a tedious but routine calculation yields(45)S0″x∗=x∗2+2μx∗cos1−α/2π+μ2−x∗+μcos1−α/2π2x∗2+2μx∗cos1−α/2π+μ23/2=μ21−cos21−α/2πx∗2+2μx∗cos1−α/2π+μ2−3/2>0.
This implies thatS˜x has no local maximum(s) and thus y0,1 also satisfy (40).We next comment in two remarks the parameter optimizations studied here for the time-fractional problems. First, in Remark1, 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 ∂tu=∂xxu−μu+f (or more generally, the advection-reaction-diffusion equation ∂tu=∂xxu+b∂xu−μu+f), the related min-max problem is(46)minq>0maxy∈y0,y1ℛInt−ordery,q,withℛInt−ordery,q=y−q2+y2−y02y+q2+y2−y02e−y,where y1>y0>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 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,ℛInt−order in (40) and ℛ in (36a) and (36b), are different(2)Because of this difference, the solutionqopt of min-max 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¯+qopt does not exist in the interval y0,y1 (see Figure 2(a)). Second, more typically, even though the local maximizer y¯+qopt exists, ℛy0,qopt≠ℛy¯+qopt,qopt (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., L=OΔx, Δt=OΔxr and Δx small, the equioscillation property always does not hold for min-max problem (36a) and (36b) when r<1/1−β and α,β∈0,1.Figure 2
Two typical behaviors ofℛy,qopt, where the equioscillation principle does not hold. Here, α, β, and μ are the problem parameters, L is the overlap size, and T is the length of time-interval.
(a)(b)Remark 2.
About the parameterq=qopt: the solution qopt of min-max problem (36a) and (36b) depends on y0,1 and ϕ, which are related to α,β,μ,T (problem parameters) and Δt,L (discretization parameters). Hence, we denote qopt by(47)qopt=Qnewα,β,μ,Δt,L,T,which we call “parameter formula” and its concrete expression is given by Lemma 5 as we have seen. Note that qopt 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 qopt is really a good approximate solution. Our evidence is based on numerical optimization:(1)
Find the best parameterqbest through repetitive numerical optimizations for original min-max problem (24a) and (24b):(i)
Forqj=qmin+jΔq, compute R˜j≔maxω∈π/T,π/Δtℛ˜w,qj by the global optimization toolbox in Matlab, where j=0,1,…,J=qmax−qmin/Δq (since sign ∂qℛ˜ω,q=signq−2LλR2ω+λI2ω, it suffices to consider q∈qmin,qmax) and(48)qmin=2Lminω∈π/T,π/ΔtλR2ω+λI2ω,qmax=2Lmaxω∈π/T,π/ΔtλR2ω+λI2ω.(ii)
Compare the valuesR˜jj=0J, and the minimal one corresponds to qbest.(2)
Compute and then compare the values ofℛ˜maxq≔maxω∈π/T,π/Δtℛ˜w,q, i.e., the convergence factor of the SWR algorithm with q=qopt (the approximate solution) and q=qbest obtained through numerical optimizations.
LetΔt=0.02, T=40, L=0.04, and μ=8. Then, in Figure 3, we plot ℛ˜maxqopt (Figure 3(a)) and ℛ˜maxqbest (Figure 3(b)) for α,β∈0,1, where a small step-size, Δq=10−3, is used in the above numerical optimization procedure to get qbest. In Figure 3(c), we plot ℛ˜maxqopt−ℛ˜maxqbest, i.e., the error between the convergence factors using the two parameters. This figure clearly shows that the quantity qopt, given by Lemma 4, is really a good approximate solution of original min-max problem (24a) and (24b).Figure 3
Comparisons ofℛ˜maxq, i.e., the convergence factor of the SWR algorithm with Robin TCs, under two choices of the parameter q, (a) q=qopt and (b) q=qbest. (c) We show the error ℛ˜maxqopt−ℛ˜maxqbest.
(a)(b)(c)


## 3.1. Towards a Good Choice ofq=2Lp: The Equioscillation Principle

A well-understood formula for choosing the free parameterp is based on the so-called equioscillation principle, established recently by Gander and Halpern [30,35]:(25)ℛ˜ωmin,q=ℛ˜ω¯q,q,where ωmin=π/T and ω¯q denotes the local maximizer of ℛ˜. The solution of (25), namely, qequi, can be obtained by solving the following two nonlinear equations:(26)ℛ˜ωmin,q=ℛ˜ω,q,∂ωℛ˜ω,q=0,together with a constraint ∂ω2ℛ˜ω,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 advection-reaction-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).Figure 1
(a) An illustration of equioscillation, i.e., the parameterq=qequi is obtained by solving (26). The plot of ℛ˜ω,q in the relevant frequency domain ω∈π/T,π/Δt for two choices of q, (b) q=qequi and (c) q=qopt given by Lemma 3.
(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 π/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Δxr, 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 Figure1(a), the parameter qequi results in a convergence factor around 0.842, which implies that the SWR algorithm using p=qequi/2L for the Robin parameter converges slowly.Since the applicability of the equioscillation principle to time-fractional problems is unclear, we introduce in Lemma4 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 ℛ˜ω,q for ω∈π/T,π/Δt, where we can see that compared to qequi the new parameter qopt leads to much smaller convergence factor and that with the new parameter the function ℛ˜ω,q does not equioscillate.


## 3.2. Towards a Good Choice ofq=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 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α,β∈0,1, μ≥0 and(27)Sx=x2+2μxcos1−α/2π+μ2−xcos1−β/2π+μcosα−β/2πx2+2μxcos1−α/2π+μ2+xcos1−β/2π+μcosα−β/2π.Then, it holds thatmaxx∈x0,x1Sx=maxSx0,Sx1, where 0<x0<x1.Proof.
Clearly, we only need to consider the caseμ>0, since for μ=0 the function Sx is a constant and therefore maxx∈x0,x1Sx=maxSx0,Sx1.
LetR0x=x2+2μxcos1−α/2π+μ2 and R1x=xcos1−β/2π+μcosα−β/2π. We then have Sx=R0x−R1x/R0x+R1x and S′x=R1xR0′x−2R0xR1′x/R0xR0x+R1x2. Suppose S′x has a positive root located at x=x∗, i.e., R1x∗R0′x−2R0x∗R1′x∗=0. Then, we have(28)S″x∗=R1xR0″x−2R0xR1″x−R1′xR0′xR0xR0x+R1x2|x=x∗.
For anyx, it holds that(29)R0′x=2x+2μcos1−α2π,R1′x=cos1−β2π,R0″x=2,R1″x=0,and this, together with μ>0 and α,β∈0,1, gives(30)S″x∗=2μcosα−β/2π−cos1−β/2πcos1−α/2πR0x∗R0x∗+R1x∗2=2μsin1−β/2πsin1−α/2πR0x∗R0x∗+R1x∗2>0.
Hence, ifx∗ corresponds to a local extremum of Sx, it must be a minimizer.Lemma 3.
Letμ≥0, α,β∈0,1, and ϕ=maxSπ/T1−α,Sπ/Δt1−α with S being the function defined in Lemma 2. Then, for all ω∈π/T,π/Δt, we have λI2ω≤ϕλR2ω.Proof.
From (22), χIω/χRω=ω1−αsin1−β/2π+μsinα−β/2π/ω1−αcos1−β/2π+μcosα−β/2π. Since α,β∈0,1,(31)cos1−β2πsinα−β2π≤sin1−β2πcosα−β2π.
Hence,ω1−αcos1−β/2π+μcosα−β/2π/ω1−αsin1−β/2π+μsinα−β/2π is a decreasing function of ω1−α; this gives(32)χIπ/ΔtχRπ/Δt≤χIωχRω≤χIπ/TχRπ/T,∀ω∈πT,πΔt.
By noticingχR2ω+χI2ω=ω21−β+μ2ω2α−β+2μω1+α−2βcos1−α/2π, we have(33)λI2ωλR2ω≤χR2ω+χI2ω−χRωχR2ω+χI2ω+χRω=Sω1−α.
Hence, from Lemma2, we get λI2ω/λR2ω≤maxSπ/T1−α,Sπ/Δt1−α=ϕ and this completes the proof.Now, for min-max problem (24a) and (24b), by using Lemma 3, we have(34)λRω−p2+λI2ωλRω+p2+λI2ωe−2LλRω≤λRω−p2+ϕλR2ωλRω+p2+ϕλR2ωe−2LλRω.For anyp>0, this leads to the following inequality:(35)maxω∈π/T,π/ΔtλRω−p2+λI2ωλRω+p2+λI2ωe−2LλRω≤maxω∈π/T,π/ΔtλRω−p2+ϕλR2ωλRω+p2+ϕλR2ωe−2LλRω.In (34), we have used p>0 and λRω>0. Let λRω=y/2L and p=q/2L. Then, the solution of the following min-max problem can be regarded as a reliable approximate solution of (24a) and (24b):(36a)minq>0maxy∈y0,y1ℛy,q,withℛy,q≔y−q2+ϕy2y+q2+ϕy2e−y,where(36b)y0≔2Lminω∈π/T,π/ΔtλRω,y1≔2Lmaxω∈π/T,π/ΔtλRω.We present the solution of (36a) and (36b) in a very general setting: we only assume ϕ∈0,1 and y1>y0>0. The concrete proof is presented in Appendix.Lemma 4.
Parameter optimization at the continuous level: assumeϕ∈0,1 and y1>y0>0. Let q0=y01+ϕ, q1=y11+ϕ, ℛmaxq=maxy∈y0,y1ℛy,q, and qopt be the solution of min-max problem (36a) and (36b). Then,(a)
Ify0≥1+ϕ1+1/ϕ−1, we have(37a)qopt=q†,ℛmaxqopt=ℛ†qopt,whereℛ†q=maxℛy0,q,ℛy1,q and q† is defined by(37b)q†=q0,ifℛy0,q0≥ℛy1,q0,q1,ifℛy1,q1≥ℛy0,q1,q˜†,otherwise.The argumentq˜†∈q0,q1 is the unique root of ℛy0,q=ℛy1,q.(b)
Ify0<1+ϕ1+1/ϕ−1<y1, we have(38a)qopt=q0∗,ifqd>q†,y¯+q†<y1andℛ¯q†>ℛ†q†,q†,otherwise,ℛmaxqopt=ℛ¯qopt,ifqd>q†,y¯+q†<y1,andℛ¯q†>ℛ†q†,ℛ†qopt,otherwise,whereqd=1+ϕ1+1/ϕ−1, y¯+q=1−ϕq2+21+ϕq+2qdq/1+ϕ, ℛ¯q=ℛy¯+q,q, dq=−ϕq2−2ϕ1+ϕq+1+ϕ2, and q0∗∈q†,qd is the unique solution of(38b)ℛ¯q=ℛy0,q.(c)
Ify1≤1+ϕ1+1/ϕ−1, then(39a)qopt=q†,if,y¯+q†≥y1,ory¯+q†<y1,andℛ¯q†≤ℛ†q†,q0∗,if,y¯+q†<y1,ℛ¯q†>ℛ†q†,andy¯+q1≥y1,or,y¯+q1<y1,andℛ¯q1<ℛ†q1,q1,if,maxy¯+q†,y¯+q1<y1,ℛ¯q†≥ℛ†q†,andℛ¯q1≥ℛ†q1,(39b)ℛmaxqopt=ℛ†qopt,if,y¯+q†≥y1,ory¯+q†<y1,andℛ¯q†≤ℛ†q†,ℛ¯qopt,otherwise,whereq0∗∈q†,q1 is the unique solution of the equation (38b).Theorem 1.
LetS be the function defined in Lemma 2 and ϕ=maxSπ/T1−α,Sπ/Δt1−α.Then, for specifiedα,β∈0,1, μ≥0, and T,Δt>0, the parameter p involved in the Robin TCs can be chosen as pc∗=qopt/2L, where qopt, given by Lemma 4, is the solution of min-max problem (36a) and (36b) and L>0 denotes the size of overlap. The quantities y0 and y1 in (36b) can be explicitly written down, as(40)y0=2LλRπT,if,α≥β,2LminλRπT,λRπΔt,if,α<β,y1=2LλRπΔt,if,α≥β,2LmaxλRπT,λRπΔt,if,α<β.Withpc∗, the convergence factor ρc defined by (19) can be bounded by(41)ρcpc∗,L≤ℛmaxqopt,where ℛmaxqopt is given by (37a), or (38a), or (39b).Proof.
bound (41) can be derived from (34) and (36b):
(34), (36b)⇒ρcp,L≤ℛmaxq≔maxy∈y0,y1ℛy,q,∀q=2Lp>0.
Hence, it remains to prove that the argumentsy0 and y1 defined by (36b) can be explicitly expressed as (40). To this end, we define(42)S0x=x2+2μxcos1−α2π+μ2+xcos1−β2π+μcosα−β2π.
Clearly, for anyα,β∈0,1 and μ≥0, the function S0x is increasing for x>0.Case 1.
α≥β. In this case, by noticing(43)2λR2ω=ωα−βS0ω1−α,we know that the function λRω is increasing for ω∈π/T,π/Δt. Hence, we have(44)y0=2Lminω∈π/T,π/ΔtλRω=L2πTα−βS0πT1−α,y1=2Lmaxω∈π/T,π/ΔtλRω=L2πΔtα−βS0πΔt1−α.Case 2.
α<β. In this case, from (43), we get 2λR2ω=S˜x≔S0xx−β−α/1−α, with x=ω1−α. We have S˜′x=xS0′x−β−α/1−αS0xx−1+β−2α/1−α and for μ=0, it holds that S˜′x=1−β/1−αS0xx−1+β−2α/1−α>0 if x>0 and α,β∈0,1; therefore, y0 and y1 satisfy (40).
We now consider the caseμ>0. Suppose the function S˜x has a local extremum located at x=x∗. Then, we have x∗S0′x∗=β−α/1−αS0x∗. Hence, it holds that S˜″x∗=x∗S0″x∗+1−β/1−αS0′x∗x∗−1+β−2α/1−α. We claim S˜″x∗>0. Indeed, from (42), it is easy to get S0′x∗>0, and this, together with α,β∈0,1, implies 1−β/1−αS0′x∗>0. For S0″x∗, a tedious but routine calculation yields(45)S0″x∗=x∗2+2μx∗cos1−α/2π+μ2−x∗+μcos1−α/2π2x∗2+2μx∗cos1−α/2π+μ23/2=μ21−cos21−α/2πx∗2+2μx∗cos1−α/2π+μ2−3/2>0.
This implies thatS˜x has no local maximum(s) and thus y0,1 also satisfy (40).We next comment in two remarks the parameter optimizations studied here for the time-fractional problems. First, in Remark1, 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 ∂tu=∂xxu−μu+f (or more generally, the advection-reaction-diffusion equation ∂tu=∂xxu+b∂xu−μu+f), the related min-max problem is(46)minq>0maxy∈y0,y1ℛInt−ordery,q,withℛInt−ordery,q=y−q2+y2−y02y+q2+y2−y02e−y,where y1>y0>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 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,ℛInt−order in (40) and ℛ in (36a) and (36b), are different(2)Because of this difference, the solutionqopt of min-max 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¯+qopt does not exist in the interval y0,y1 (see Figure 2(a)). Second, more typically, even though the local maximizer y¯+qopt exists, ℛy0,qopt≠ℛy¯+qopt,qopt (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., L=OΔx, Δt=OΔxr and Δx small, the equioscillation property always does not hold for min-max problem (36a) and (36b) when r<1/1−β and α,β∈0,1.Figure 2
Two typical behaviors ofℛy,qopt, where the equioscillation principle does not hold. Here, α, β, and μ are the problem parameters, L is the overlap size, and T is the length of time-interval.
(a)(b)Remark 2.
About the parameterq=qopt: the solution qopt of min-max problem (36a) and (36b) depends on y0,1 and ϕ, which are related to α,β,μ,T (problem parameters) and Δt,L (discretization parameters). Hence, we denote qopt by(47)qopt=Qnewα,β,μ,Δt,L,T,which we call “parameter formula” and its concrete expression is given by Lemma 5 as we have seen. Note that qopt 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 qopt is really a good approximate solution. Our evidence is based on numerical optimization:(1)
Find the best parameterqbest through repetitive numerical optimizations for original min-max problem (24a) and (24b):(i)
Forqj=qmin+jΔq, compute R˜j≔maxω∈π/T,π/Δtℛ˜w,qj by the global optimization toolbox in Matlab, where j=0,1,…,J=qmax−qmin/Δq (since sign ∂qℛ˜ω,q=signq−2LλR2ω+λI2ω, it suffices to consider q∈qmin,qmax) and(48)qmin=2Lminω∈π/T,π/ΔtλR2ω+λI2ω,qmax=2Lmaxω∈π/T,π/ΔtλR2ω+λI2ω.(ii)
Compare the valuesR˜jj=0J, and the minimal one corresponds to qbest.(2)
Compute and then compare the values ofℛ˜maxq≔maxω∈π/T,π/Δtℛ˜w,q, i.e., the convergence factor of the SWR algorithm with q=qopt (the approximate solution) and q=qbest obtained through numerical optimizations.
LetΔt=0.02, T=40, L=0.04, and μ=8. Then, in Figure 3, we plot ℛ˜maxqopt (Figure 3(a)) and ℛ˜maxqbest (Figure 3(b)) for α,β∈0,1, where a small step-size, Δq=10−3, is used in the above numerical optimization procedure to get qbest. In Figure 3(c), we plot ℛ˜maxqopt−ℛ˜maxqbest, i.e., the error between the convergence factors using the two parameters. This figure clearly shows that the quantity qopt, given by Lemma 4, is really a good approximate solution of original min-max problem (24a) and (24b).Figure 3
Comparisons ofℛ˜maxq, i.e., the convergence factor of the SWR algorithm with Robin TCs, under two choices of the parameter q, (a) q=qopt and (b) q=qbest. (c) We show the error ℛ˜maxqopt−ℛ˜maxqbest.
(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 t, we get(49)ujk−ℐt1−β∂xxujk+μℐt1−αujk=f˜x,t,x∈Ωj,∂x+−1j+1pujk2−jL,t=∂x+−1j+1pu3−jk−12−jL,t,ujkx,0=u0,where f˜=∫0tfx,τdτ+u0x. For any γ∈0,1, ℐtγ denotes fractional integral(50)ℐtγvt=vt,ifγ=0,1Γγ∫0t1t−τ1−γvτdτ,ifγ∈0,1.


### 4.1. Derivation of the Discrete SWR Algorithms

To get discrete algorithm, we apply the centered finite difference formula to (49), which gives(51a)Uj,mkt−ℐt1−βδx2Uj,mkt+μℐt1−αUj,mkt=f˜xm,t,xm∈Ωj,j=1,2,where Uj,mkt≈ujkxm,t, δx2Uj,mkt=Uj,m+1kt−2Uj,mkt+Uj,m−1kt/Δx2, xm=mΔx, and m∈ℤ. To simplify the convergence analysis, we assume L=Δx, i.e., the size of overlap equals a single space mesh size. Then, the artificial boundary interfaces x=L (resp. x=0) corresponds to x1 (resp. x0) and the Robin TCs are discretized as(51b)U1,1kt−U1,0ktΔx+pU1,1kt=U2,1k−1t−U2,0k−1tΔx+pU2,1k−1t,U2,1kt−U2,0ktΔx−pU2,0kt=U1,1k−1t−U1,0k−1tΔx−pU1,0k−1t.Here, we discretized the spatial derivatives in the Robin TCs as(51c)∂xu1kx1,t≈U1,1kt−U1,0ktΔx,∂xu2kx0,t≈U2,1kt−U2,0ktΔx.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(52)∫0tKγt−τVτdτ,with kernel functionKγt=tγ−1Γγ,where γ=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 ONt2 operations and ONt memory for per expansion coefficients over Nt time steps. In [44], fast convolution algorithms are constructed, which implements, over Nt time steps, with ONtlogNt operations and OlogNt active memory.The starting point of the fast convolution algorithms is to replace the kernel functionKγt−τ in (52) by its inverse Laplace transform, along some contour Υ:(53)∫0t12πi∫Υest−τK^γsds︸Kγt−τVτdτ,where K^γs=ℒKγ=s−γ denotes the Laplace forward transform of the kernel function Kγ. 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 ets decays rapidly for growing s along Υ. In (53), reversing the order of integration gives(54)12πi∫ΥK^γs∫0test−τVτdτ︸≔y0,t;sds=12πi∫ΥK^γsy0,t;sds.The functiony0,t;s is the solution of the following linear scale ODE at τ=t:(55)y′τ=syτ+Vτ,τ∈0,t,y0=0.SinceVt is only available at discrete grids tnn=0Nt, it is reasonable to apply some time-integrator, e.g., the linear multistep method, to (55) to get y at t=tn:(56)∑j=0Kajyn+j−K=Δt∑j=0Kbjsyn+j−K+Vn+j−K.Multiplying (56) by ζn and summing over n from 0 to +∞, we obtain(57)∑j=0KajζK−jyζ=∑j=0KbjζK−jΔtsyζ+ΔtVζ,with the generating (formal) power series yζ=∑n=0+∞ynζn and Vζ=∑n=0+∞Vnζn. Solving yζ from (57), we get(58)ΛζΔt−syζ=Vζ,with,Λζ=a0ζK+⋯+aK−1ζ+aKb0ζK+⋯+bK−1ζ+bK,where Λζ is called the quotient of the generating polynomial of a linear multistep method. Clearly, the solution yn of (56) is the n-th coefficient of the power series Λζ/Δt−s−1Vζ, that is, y0,tn;s≈n−th coefficient ofΛζ/Δt−s−1Vζ. By substituting this into the right-hand side of (54), we get(59)12πi∫ΥK^γs∫0tnest−τVτdτds≈n−th coefficient of12πi∫ΥK^γsΛζΔt−s−1ds︸using Cauchy'sintegral formula to this partVζ=n−th coefficient ofK^γΛζΔtVζ.Now, we expandK^γΛζ/Δt in the formal power series, as(60)K^γΛζΔt=∑j=0+∞cjγ,Δtζj.Then, then-th coefficient of K⌢γΛζ/ΔtVζ is given by(61)∑j=0ncn−jγ,ΔtVj.Suppose the time-integrator applied to (55) is an A-stable method with order p. Then, from [40], Theorem 2.1, we know that (61) is an approximation to (52) and(62)∫0tnKγtn−τVτdτ=∑j=0ncn−jγ,ΔtVj+Ohp,∀n≥1.The last step is to determine the weightscjγ,Δtj=0+∞. For general case, we have to rely on numerical calculation to get cjγ,Δtj=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, Λζ=∑j=121−ζj/j), these weights can be written down, as(63)θ−method:cjγ,Δt=θΔtγ1−θθj∑l=0jθθ−1l−γlγj−l,θ∈0.5,1,Δtγ−1j−γl,θ=1,BDF2:cjγ,Δt=2Δt3γ−1j∑l=0j3−l−γj−l−γl.With the quadrature weightscnγ,Δtn=0+∞ in our hands, we denote the discrete convolution (61) by Kγ,ΔtVn, that is,(64)Kγ,ΔtVn=∑l=0ncn−lγ,ΔtVl.Then, the scheme of the fully discrete SWR algorithms on thej-th subdomain is(65)Uj,mkn−K1−β,Δtδx2Uj,mkn+μK1−α,ΔtUj,mkn=f˜xm,tn,xm∈Ωj.Together with boundary conditions (51b) at t=tn.


### 4.2. Convergence Analysis

Since the quadrature weightscnγ,Δtn=0+∞ are the coefficients of the formal power series of K⌢γΛζ/Δt, the discrete Laplace transform of the sequence Kγ,ΔtVnn=0∞ defined by (64) is K^γΛz/ΔtVz, where Vz=∑n=0∞Vnzn and z=e−Δts. Here, s=σ+iω with σ>0 and ω∈π/T,π/Δt. Define the error sequence as Ej,mknn=0+∞=Uj,mkn−Uj,mnn=0+∞. Then, the discrete Laplace transform of Ej,mknn=0+∞, which we denote by Ej,mkz, satisfies(66)Ej,mkz−K^1−βΛzΔtδx2Ej,mkz+μK^1−αΛzΔtEj,mkz=0,xm∈Ωj,i.e.,(67)Ej,m−1kz−2as˜Ej,mkz+Ej,m+1kz=0,where(68)s˜=ΛzΔt,as˜=1+Δx22s˜1−β+μs˜α−β.Define(69)η=1+Δxp.Then, interface conditions (51b) are transformed to(70)ηE1,1kz−E1,0kz=ηE2,1k−1z−E2,0k−1z,ηE2,0kz−E2,1kz=ηE1,0k−1z−E1,1k−1z.From (67), the arguments Ej,mk can be expressed as(71)Ej,mk=Ajkξ+m+Bjkξ−m,withξ±s˜=as˜±a2s˜−1.Lemma 5.
Letσ, Δt, T>0, ϵ=e−σΔt, and cmax=cosπ/T. Then, for z=e−sΔt with s=σ+iω and ω∈π/T,π/Δt, we have(1)
ifΛz=1−z/θ+1−θz (linear θ-method):(i)
ℜΛz∈1−cmax2θ−1/1−θ2+θ2+2θ1−θcmax,2/2θ−1 and ℑΛz≤1/1−2θ2, if θ∈1/2,1;(ii)
ℜΛz∈21−ϵ/1+ϵ,21+ϵ/1−ϵ and ℑΛz≤41+ϵ2/1−ϵ22sinarccos2ϵ/1+ϵ2, if θ=1/2;(2)
ifΛz=∑j=121−zj/j (BDF2 method): ℜΛz∈1−cmax2,3/2, ℑΛz≤2.202.Proof.
Letc=cosωΔt and c˜=sinωΔt. We have, ϵ∈0,1, c∈−1,cmax, and c˜∈−1,1. The proof for the two time-integrators is given as follows:(1)
For the linearθ−method: Λz=1−z/θ+1−θz. In this case, we have (72)Λz=1−ϵc+ic˜θ+ϵ1−θc+ic˜=ΛRϵ,c−iΛIϵ,c,c˜=Nϵ,cDϵ,c−iϵc˜Dϵ,c,whereNϵ,c=−ϵ21−θ+ϵ1−2θc+θ and Dϵ,c=ϵ21−θ2+2ϵθ1−θc+θ2. For θ=0.5, the bounds for ℜΛ and ℑΛ can be obtained easily. For θ∈0.5,1,(73)sign∂ϵΛR=sign−ϵ21−θ2c−2ϵ1−θθ−θ2c,sign∂ϵΛI=signc˜θ1−θ−ϵ=signc˜,sign∂cΛR=signϵ−θ1−θ,sign∂cΛI=sign−c˜.Clearly, it holds that∂cΛR<0 and this implies(74)maxc∈−1,cmaxΛRϵ,c=Nϵ,−1Dϵ,−1,minc∈−1,cmaxΛRϵ,c=Nϵ,cmaxDϵ,cmax.Moreover, since∂ϵΛRϵ,−1>0 and ∂ϵΛRϵ,cmax<0, we have(75)maxϵ∈0,1maxc∈−1,cmaxΛRϵ,c=N1,−1D1,−1=22θ−1,minϵ∈0,1minc∈−1,cmaxΛRϵ,c=N1,cmaxD1,cmax>N1,1D1,1=0.This gives the bounds ofℜΛ and ℑΛ for θ∈0.5,1.(2)
For the BDF2 method:Λz=∑j=121−zj/j. In this case, we have(76)Λz=123−4εc+ε22c2−1−iεc˜2−εc.Hence,ℜΛz∈1−cmax2,3/2. For ℑΛz, by noticing(77)ℑΛz=1−c22ε−ε2c≤1−c22−c,forε∈0,1,it is easy to getℑΛz≤2.202.Lemma 6.
Letα,β∈0,1, μ≥0, and ℍ≔s∈ℂ+:ℑs∈π/T,π/Δt. Then, for the linear θ-method with θ∈1/2,1 and the BDF2 method, the arguments ξ±s˜ with s˜=e−sΔt are analytic for s∈ℍ and satisfy ξ+>1 and ξ−<1.Proof.
Fors∈ℍ, from Lemma 5, we have s˜=Λe−sΔt/Δt∈ℂ+. Then, for μ≥0, it holds that ℜas˜>1 for s˜∈ℂ+, where as˜ is defined by (68). This implies ℜ1+a2s˜>1, and thus, the composite functions ξ±s˜ are also analytic. Since ℜ1+a2s˜>1, we have ℜξ+s˜>1⇒ξ+s˜>1; this, together with ξ+s˜ξ−s˜≡1, gives ξ−s˜<1.Now, by using the boundedness conditions at infinity for each subsystem, the constantsAk and Bk in (71) can be determined as Bk=0 and Ak=0 for the first and second subsystems, respectively. We then have E1,mk=Akξ+m and E2,mk=Bkξ−m. By using the transmission conditions in (70), we get the recurrence relation Akηξ+−1=Bk−1ηξ−−1 and Bkη−ξ−=Ak−1η−ξ+, i.e.,(78a)Ak=ρ^dsAk−2,Bk=ρ^dsBk−2,withρ^ds=η−ξ+s˜ηξ+s˜−12,where the subscript “d” denotes “discrete” and(78b)s˜=Λe−sΔtΔt,s∈ℍ≔s∈ℂ+:ℑs∈π/T,π/Δt.Mathematically, we expectρd≪1, which leads to the following min-max problem:(79)minη∈ℝmaxs˜=Λe−sΔt/Δt,s∈ℍη−ξ+s˜ηξ+s˜−12.Letξ+s˜=ξRs˜+iξIs˜. Then, from Lemma 6, we know ξRs˜>1. For η≤1, we have ξ+s˜2−11−η2≥0⇒η−ξRs˜2+ξI2s˜/ηξRs˜−12+ηξI2s˜≥1⇒η−ξ+s˜/ηξ+s˜−12≥1. This implies that, it is sufficient to consider η>1 in (79).Forθ=0.5, the quantity s˜=Λe−sΔt/Δt approaches to the origin, as s approaches to 0,±π in the complex plane. Hence, ξ+ approaches to 1 as s⟶0,±π. This implies that maxs∈ℍρ^ds≡1, independent of η.Forθ∈0.5,1, we know from Lemma 5 that the origin is kept away from s˜, and therefore, it is possible to expect maxs∈ℍρ^ds<1. 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,(80)Dσ¯min,σ¯max,ω¯max≔s∈ℂ+:σ¯min≤ℜs≤σ¯max,0≤ℑs≤ω¯max,which lies in the right-hand side of the vertical line s=1 (i.e., σ¯min>1) 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.Figure 4
An illustration of the region mapped byξ+s˜=as˜+a2s˜−1 (with s˜=Λe−sΔt/Δt and s∈ℍ) and the rectangular region Dσ¯min,σ¯max,ω¯max covering the former. Here, we consider the backward Euler method; for the BDF2 method, the plot looks similar.For the linearθ-method with θ∈0.5,1 and the BDF2 method, if α≥β, by using Lemma 5 and Lemma 7 given below, the boundaries of the rectangular region D, i.e., σ¯min, σ¯max, and ω¯max, can be expressed explicitly, as(81a)σ¯min=σ˜0+σ˜02−1,σ¯max=ξRσ˜1,ω˜max,ω¯max=ξIσ˜0,ω˜max,where ξR,I are defined by (87a). The quantities σ˜0, σ˜1, and ω˜max are defined by(81b)σ˜0=1+Δx22σ01−β+μσ0α−β,σ˜1=1+Δx22σ12+ωmax21−β/2cos1−βψ1+μσ12+ωmax2α−β/2cosα−βψ1,ω˜max=Δx22σ02+ωmax21−β/2cos1−βψ0+μσ02+ωmax2α−β/2cosα−βψ0,withψ0=arccosσ0σ02+ωmax2,ψ1=arccosσ1σ12+ωmax2.The argumentsσ0,1 and ωmax, appearing in (81b), are given by(81c)σ0=1−cmax2θ−1Δt1−θ2+θ2+2θ1−θcmax,σ1=22θ−1Δt,ωmax=12θ−12Δt,θ−method,σ0=1−cmax2,σ1=32,ωmax=2.202,BDF2,where θ∈1/2,1 and cmax=cosπ/T.In general, we have to rely on numerical optimizations to obtain these three quantities, i.e., σ¯min=mins˜=Λe−sΔt/Δt,s∈ℍℜξ+s˜, σ¯max=maxs˜=Λe−sΔt/Δt,s∈ℍℜξ+s˜, and ω¯max=maxs˜=Λe−sΔt/Δt,s∈ℍℜξ+s˜.Lemma 7.
α≥β. Assume 0<β≤α<1 and μ≥0. Let D0=s∈ℂ+:0<σ0≤ℜs≤σ1,0≤ℑs≤ωmax and D=s∈ℂ+:σ¯min≤ℜs≤σ¯max,0≤ℑs≤ω¯max, where σ¯min,max and ω¯max are given by (81a) and (81b). Then, we have ξ+s∈D for s∈D0.Proof.
With the argumentsσ˜0, σ˜1, and ω˜max given by (81b), we define D˜=s∈ℂ+:0<σ˜0≤ℜs≤σ˜1,0≤ℑs≤ω˜max. Then, the proof is divided into two steps: as∈D˜ if s∈D0 and s+s2−1∈D˜ if s∈D˜. The function as is defined by (68).Step 1.
as∈D˜ for s∈D0. Let s=σ+iω∈D0. Then, we can rewrite as as(82)as=1+Δx22eGR1−β,r,σ+μeGRα−β,r,σ+iΔx22eGI1−β,r,σ+μeGIα−β,r,σ,where r=σ2+ω2. For any γ∈0,1, the functions GR,I are defined by (83)GRγ,r,σ=γlogr+logcosγψr,σ,GIγ,r,σ=γlogr+logsinγψr,σ,where ψr,σ=arccosσ/r. Since σ>0, we have ψr∈0,π/2. For fixed σ∈σ0,σ1, because ψr,σ is an increasing function of r, GI is also an increasing function of r. Hence, for γ∈0,1, maxr∈σ,σ2+ωmax2GIγ,r,σ=GIγ,σ2+ωmax2,σ. Define ϑr=arcsinωmax/r. Then, GIγ,σ2+ωmax2,σ = γlogr+logsinγϑr for r=ωmax2+σ2. We have(84)sign∂rγlogr+logsinγϑr=signrsinγϑr−ωmax.
Hence, forγ∈0,1, it holds that ∂rγlogr+logsinγϑr≤0 and this implies that the function GIγ,σ2+ωmax2,σ attains its maximum at σ=σ0. We therefore get the upper bound of ℑas, i.e., ω˜max, as shown in (81b).
It remains to consider the lower and upper bounds ofℜas. For r=σ2+ω2 and σ∈σ0,σ1, we have sign∂rGRγ,r,σ=signrcosγψr−σ. Hence, for γ∈0,1, it holds that ∂rGRγ,r,σ≥0 and this gives(85)minσ≤r≤ωmax2+σ2GRγ,r,σ=γlogσ,maxσ≤r≤ωmax2+σ2GRγ,r,σ=GRγ,ωmax2+σ2,σ.
The first result gives the lower bound ofℜas, i.e., σ˜0. With the function ϑr=arcsinωmax/r, the argument GRγ,ωmax2+σ2,σ can be represented as γlogr+logcosγϑr with r=ωmax2+σ2. We have(86)sign∂rγlogr+logcosγϑr=sign1+ωmaxr2−ωmax2sinγϑrcosγϑr≥0.
Hence,GR attains its global maximum at σ=σ1; this gives σ˜1, the upper bound of ℜas.Step 2.
ξ+s∈D˜ for s∈D˜. Let s=σ+iω with σ∈σ˜0,σ˜1 and ω∈0,ω˜max. Then, ξ+s=ξRσ,ω+iξIσ,ω, where(87a)ξRσ,ω=σ+σ2−ω2−12+4σ2ω2+σ2−ω2−12,ξIσ,ω=ω+σ2−ω2−12+4σ2ω2−σ2−ω2−12.
Clearly, it holds∂ωξIσ,ω>0 and ∂σξRσ,ω>0. Moreover, some simple algebra shows(87b)sign∂ω2ξRσ,ω=signσ2+ω2+1−σ2−ω2−12+4σ2ω2≥0,sign∂σ2ξIσ,ω=signσ2+ω2−1−σ2−ω2−12+4σ2ω2≥0.
Hence,maxω∈0,ω˜maxξIσ,ω=ξIσ,ω˜max and from the second inequality given in (87b), we know that ξIσ,ω˜max attains its maximum at σ=σ˜0; this gives the upper bound of ℑξ+s, i.e., ξIσ˜0,ω˜max. Similarly, we have(88)minσ∈σ˜0,σ˜1minω≤ω˜maxξRσ,ω=ξRσ˜0,0=σ˜0+σ˜02−1,maxσ∈σ˜0,σ˜1maxω≤ω˜maxξRσ,ω=ξRσ˜1,ω˜max.
Combining the analysis in Step1 and Step 2, we complete the whole proof.With the rectangular regionD covering ξ+, we have(89)maxs˜=Λe−sΔt/Δt,s∈ℍη−ξ+s˜ηξ+s˜−12≤maxs∈Dσ¯min,σ¯max,ω¯maxη−sηs−12,and based on this relation, it is reasonable to consider the following equation:(90)minη>1maxs∈Dσ¯min,σ¯max,ω¯maxη−sηs−12,withσ¯min>1.To analyze the solution of this min-max problem, we define(91a)κx,y=Z+Z2−42,withZ=1+x2+y2x,η−=κσ¯max,0,η+=κσ¯min,ω¯max,Tσ,ω,η=σ−η2+ω2ση−12+η2ω2.Then, forηmin<η−≤η+<ηmax, we defined two functions T−η and T+η as(91b)T+η=Tσ¯min,0,η,η∈ηmin,η−,maxTσ¯min,0,η,Tσ¯max,0,η,η∈η−,η+,maxTσ¯min,0,η,Tσ¯max,0,η,Tσ¯min,ω¯max,η,η∈η+,ηmax,T−η=maxTσ¯max,ω¯max,η,Tσ¯max,0,η,Tσ¯min,ω¯max,η,η∈ηmin,η−,maxTσ¯max,ω¯max,η,Tσ¯min,ω¯max,η,η∈η−,η+,Tσ¯max,ω¯max,η,η∈η+,ηmax.Forηmin<η+<η−<ηmax, these two functions are defined by(91c)T+η=Tσ¯min,0,η,η∈ηmin,η+,maxTσ¯min,0,η,Tσ¯min,ω¯max,η,η∈η+,η−,maxTσ¯min,0,η,Tσ¯max,0,η,Tσ¯min,ω¯max,η,η∈η−,ηmax,T−η=maxTσ¯max,ω¯max,η,Tσ¯max,0,η,Tσ¯min,ω¯max,η,η∈ηmin,η+,maxTσ¯max,ω¯max,η,Tσ¯max,0,η,η∈η+,η−,Tσ¯max,ω¯max,η,η∈η−,ηmax.Theorem 2.
Parameter optimization at the discrete level: letL=Δx, α,β∈0,1 and μ≥0. Then, the parameter p involved in the discrete SWR algorithms can be chosen as pd∗=ηopt−1/Δx, where ηopt, the solution of min-max problem (90), is given by(92)ηopt=ηmin,if,T+ηmin≥T−ηmin,ηmax,if,T−ηmax≥T+ηmax,η∗,otherwise,where ηmin=κσ¯min,0, ηmax=κσ¯max,ω¯max, and η∗∈ηmin,ηmax are the unique root of T+η=T−η. With the parameter pd∗, the convergence factor ρd≔maxs∈ℍρ^ds at the discrete level satisfies ρd≤maxT+ηopt,T−ηopt, where ρ^d is given by (78a).Proof.
Lets=σ+iω∈Dσ¯min,σ¯max,ω¯max. Then, it holds that η−s/ηs−12=Tσ,ω,η. For η>1, a routine calculation yields(93)sign∂ηTσ,ω,η=signη−κσ,ω.
Hence, the solution of (90) should satisfy ηopt∈ηmin,ηmax, since otherwise, if ηopt<ηmin (resp. ηopt>ηmax), we can further minimize T by increasing (resp. decreasing) ηopt.
We havesign∂ωTσ,ω,η=sign2ση−1−η2, and therefore, for fixed σ, it holds that maxω∈0,ω¯maxTσ,ω,η=maxTσ,0,η,Tσ,ω¯max,η. Hence, maxs∈Dσ¯min,σ¯max,ω¯maxTσ,ω,η=maxσ∈σ¯min,σ¯maxmaxTσ,0,η,Tσ,ω¯max,η. Another routine calculation yields(94)∂σTσ,ω,η=2η2−1ηση−12+η2ω22σ2−1+η2ησ+1−ω2.
Then, for fixedη>1 and ω≥0, if we vary σ from 1 to +∞, we first meet ∂σTσ,ω,η<0, and then ∂σTσ,ω,η=0, and then ∂σTσ,ω,η>0. Hence, for fixed ω, the function Tσ,ω,η does not have local maximum(s) for σ∈σ¯min,σ¯max. This implies that for given η>1 the global maximum of T on the region D is obtained at the four corners, i.e.,(95)maxs∈DTσ,ω,η=maxTσ¯min,0,η,Tσ¯max,0,η,Tσ¯min,ω¯max,η,Tσ¯max,ω¯max,η.
With the functionsT±η defined by (91b) and (91c), we have(96)minη∈ηmin,ηmaxmaxs∈Dσ¯min,σ¯max,ω¯maxTσ,ω,η=minη∈ηmin,ηmaxmaxT+η,T−η.
It remains to analyze the monotonic properties ofT+η and T−η. From (93), we know that Tσ¯min,0,η and Tσ¯max,ω¯max,η are, respectively, increasing and decreasing for η∈ηmin,ηmax. Similarly, Tσ¯max,0,η is decreasing for η∈ηmin,η− and increasing for η∈η−,ηmax; Tσ¯min,ω¯max,η is decreasing for η∈ηmin,η+ and increasing for η∈η+,ηmax. Hence, the functions T+η and T−η, defined by (91b) and (91c), are, respectively, increasing and decreasing for η∈ηmin,ηmax. This, together with some simple algebra, gives the solution of (96), i.e., ηopt, 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 D, i.e., σ¯min, σ¯max, and ω¯max(2)
Computeηopt according to formula (92) and then let pd∗=ηopt−1/Δx
The quotientΛζ associated with a specified time-integrator determines the boundaries of D, which further determines the parameter ηopt. To make a convergent SWR algorithm, Λe−sΔt should lay in ℂ+ and should be away from the origin, ∀s∈s∈ℂ+:ℑs∈π/T,π/Δt.Remark 4.
Large overlap size: if we use an overlap sizeL=MΔx with M≥2, the following min-max problem occurs:(97)minη>1maxs∈Dσ¯min,σ¯max,ω¯maxTs,η,withTs,η=η−ssM−1ηs−12.
ForM≥2 and η>1, the function Ts,η has local maximums in the interior of D, and therefore, the analysis is much more complicated compared to the case M=1. 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 Figure5, consisting of infinite number of series connected cell-circuits, the state (or modified-nodal-analysis) circuit equation is(98)dbxtdtb=Axt+ft,A=c2⋱⋱⋱1−211−21⋱⋱⋱,Figure 5
(a) The fractional RC circuit in infinite sizeN=∞; (b) a realization of the floating fractional capacitor, the main electronic component of the circuit shown in (a).
(a)(b)wherext=⋯,x−1t,x0t,x1t,⋯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). The 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 xt=ℐtbAxt+f˜t, with ℐtb 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 Uj,mkt by xmt and δx2 by A 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 with(99)s~=s∈H≔s∈C+:Is∈πT,πΔt,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, Theorem 2 is still applicable, but we need to replace ηmax by +∞ there and redefine T±η as(100)T+η=Tσ¯min,0,η,η∈ηmin,η+,maxTσ¯min,0,η,Tσ¯min,ω¯max,η,η>η+,T−η=maxη−2,Tσ¯min,ω¯max,η,η∈ηmin,η+,η−2,η>η+,where η+ and ηmin are the same quantities as we have stated in Theorem 2.


## 4.1. Derivation of the Discrete SWR Algorithms

To get discrete algorithm, we apply the centered finite difference formula to (49), which gives(51a)Uj,mkt−ℐt1−βδx2Uj,mkt+μℐt1−αUj,mkt=f˜xm,t,xm∈Ωj,j=1,2,where Uj,mkt≈ujkxm,t, δx2Uj,mkt=Uj,m+1kt−2Uj,mkt+Uj,m−1kt/Δx2, xm=mΔx, and m∈ℤ. To simplify the convergence analysis, we assume L=Δx, i.e., the size of overlap equals a single space mesh size. Then, the artificial boundary interfaces x=L (resp. x=0) corresponds to x1 (resp. x0) and the Robin TCs are discretized as(51b)U1,1kt−U1,0ktΔx+pU1,1kt=U2,1k−1t−U2,0k−1tΔx+pU2,1k−1t,U2,1kt−U2,0ktΔx−pU2,0kt=U1,1k−1t−U1,0k−1tΔx−pU1,0k−1t.Here, we discretized the spatial derivatives in the Robin TCs as(51c)∂xu1kx1,t≈U1,1kt−U1,0ktΔx,∂xu2kx0,t≈U2,1kt−U2,0ktΔx.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(52)∫0tKγt−τVτdτ,with kernel functionKγt=tγ−1Γγ,where γ=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 ONt2 operations and ONt memory for per expansion coefficients over Nt time steps. In [44], fast convolution algorithms are constructed, which implements, over Nt time steps, with ONtlogNt operations and OlogNt active memory.The starting point of the fast convolution algorithms is to replace the kernel functionKγt−τ in (52) by its inverse Laplace transform, along some contour Υ:(53)∫0t12πi∫Υest−τK^γsds︸Kγt−τVτdτ,where K^γs=ℒKγ=s−γ denotes the Laplace forward transform of the kernel function Kγ. 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 ets decays rapidly for growing s along Υ. In (53), reversing the order of integration gives(54)12πi∫ΥK^γs∫0test−τVτdτ︸≔y0,t;sds=12πi∫ΥK^γsy0,t;sds.The functiony0,t;s is the solution of the following linear scale ODE at τ=t:(55)y′τ=syτ+Vτ,τ∈0,t,y0=0.SinceVt is only available at discrete grids tnn=0Nt, it is reasonable to apply some time-integrator, e.g., the linear multistep method, to (55) to get y at t=tn:(56)∑j=0Kajyn+j−K=Δt∑j=0Kbjsyn+j−K+Vn+j−K.Multiplying (56) by ζn and summing over n from 0 to +∞, we obtain(57)∑j=0KajζK−jyζ=∑j=0KbjζK−jΔtsyζ+ΔtVζ,with the generating (formal) power series yζ=∑n=0+∞ynζn and Vζ=∑n=0+∞Vnζn. Solving yζ from (57), we get(58)ΛζΔt−syζ=Vζ,with,Λζ=a0ζK+⋯+aK−1ζ+aKb0ζK+⋯+bK−1ζ+bK,where Λζ is called the quotient of the generating polynomial of a linear multistep method. Clearly, the solution yn of (56) is the n-th coefficient of the power series Λζ/Δt−s−1Vζ, that is, y0,tn;s≈n−th coefficient ofΛζ/Δt−s−1Vζ. By substituting this into the right-hand side of (54), we get(59)12πi∫ΥK^γs∫0tnest−τVτdτds≈n−th coefficient of12πi∫ΥK^γsΛζΔt−s−1ds︸using Cauchy'sintegral formula to this partVζ=n−th coefficient ofK^γΛζΔtVζ.Now, we expandK^γΛζ/Δt in the formal power series, as(60)K^γΛζΔt=∑j=0+∞cjγ,Δtζj.Then, then-th coefficient of K⌢γΛζ/ΔtVζ is given by(61)∑j=0ncn−jγ,ΔtVj.Suppose the time-integrator applied to (55) is an A-stable method with order p. Then, from [40], Theorem 2.1, we know that (61) is an approximation to (52) and(62)∫0tnKγtn−τVτdτ=∑j=0ncn−jγ,ΔtVj+Ohp,∀n≥1.The last step is to determine the weightscjγ,Δtj=0+∞. For general case, we have to rely on numerical calculation to get cjγ,Δtj=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, Λζ=∑j=121−ζj/j), these weights can be written down, as(63)θ−method:cjγ,Δt=θΔtγ1−θθj∑l=0jθθ−1l−γlγj−l,θ∈0.5,1,Δtγ−1j−γl,θ=1,BDF2:cjγ,Δt=2Δt3γ−1j∑l=0j3−l−γj−l−γl.With the quadrature weightscnγ,Δtn=0+∞ in our hands, we denote the discrete convolution (61) by Kγ,ΔtVn, that is,(64)Kγ,ΔtVn=∑l=0ncn−lγ,ΔtVl.Then, the scheme of the fully discrete SWR algorithms on thej-th subdomain is(65)Uj,mkn−K1−β,Δtδx2Uj,mkn+μK1−α,ΔtUj,mkn=f˜xm,tn,xm∈Ωj.Together with boundary conditions (51b) at t=tn.


## 4.2. Convergence Analysis

Since the quadrature weightscnγ,Δtn=0+∞ are the coefficients of the formal power series of K⌢γΛζ/Δt, the discrete Laplace transform of the sequence Kγ,ΔtVnn=0∞ defined by (64) is K^γΛz/ΔtVz, where Vz=∑n=0∞Vnzn and z=e−Δts. Here, s=σ+iω with σ>0 and ω∈π/T,π/Δt. Define the error sequence as Ej,mknn=0+∞=Uj,mkn−Uj,mnn=0+∞. Then, the discrete Laplace transform of Ej,mknn=0+∞, which we denote by Ej,mkz, satisfies(66)Ej,mkz−K^1−βΛzΔtδx2Ej,mkz+μK^1−αΛzΔtEj,mkz=0,xm∈Ωj,i.e.,(67)Ej,m−1kz−2as˜Ej,mkz+Ej,m+1kz=0,where(68)s˜=ΛzΔt,as˜=1+Δx22s˜1−β+μs˜α−β.Define(69)η=1+Δxp.Then, interface conditions (51b) are transformed to(70)ηE1,1kz−E1,0kz=ηE2,1k−1z−E2,0k−1z,ηE2,0kz−E2,1kz=ηE1,0k−1z−E1,1k−1z.From (67), the arguments Ej,mk can be expressed as(71)Ej,mk=Ajkξ+m+Bjkξ−m,withξ±s˜=as˜±a2s˜−1.Lemma 5.
Letσ, Δt, T>0, ϵ=e−σΔt, and cmax=cosπ/T. Then, for z=e−sΔt with s=σ+iω and ω∈π/T,π/Δt, we have(1)
ifΛz=1−z/θ+1−θz (linear θ-method):(i)
ℜΛz∈1−cmax2θ−1/1−θ2+θ2+2θ1−θcmax,2/2θ−1 and ℑΛz≤1/1−2θ2, if θ∈1/2,1;(ii)
ℜΛz∈21−ϵ/1+ϵ,21+ϵ/1−ϵ and ℑΛz≤41+ϵ2/1−ϵ22sinarccos2ϵ/1+ϵ2, if θ=1/2;(2)
ifΛz=∑j=121−zj/j (BDF2 method): ℜΛz∈1−cmax2,3/2, ℑΛz≤2.202.Proof.
Letc=cosωΔt and c˜=sinωΔt. We have, ϵ∈0,1, c∈−1,cmax, and c˜∈−1,1. The proof for the two time-integrators is given as follows:(1)
For the linearθ−method: Λz=1−z/θ+1−θz. In this case, we have (72)Λz=1−ϵc+ic˜θ+ϵ1−θc+ic˜=ΛRϵ,c−iΛIϵ,c,c˜=Nϵ,cDϵ,c−iϵc˜Dϵ,c,whereNϵ,c=−ϵ21−θ+ϵ1−2θc+θ and Dϵ,c=ϵ21−θ2+2ϵθ1−θc+θ2. For θ=0.5, the bounds for ℜΛ and ℑΛ can be obtained easily. For θ∈0.5,1,(73)sign∂ϵΛR=sign−ϵ21−θ2c−2ϵ1−θθ−θ2c,sign∂ϵΛI=signc˜θ1−θ−ϵ=signc˜,sign∂cΛR=signϵ−θ1−θ,sign∂cΛI=sign−c˜.Clearly, it holds that∂cΛR<0 and this implies(74)maxc∈−1,cmaxΛRϵ,c=Nϵ,−1Dϵ,−1,minc∈−1,cmaxΛRϵ,c=Nϵ,cmaxDϵ,cmax.Moreover, since∂ϵΛRϵ,−1>0 and ∂ϵΛRϵ,cmax<0, we have(75)maxϵ∈0,1maxc∈−1,cmaxΛRϵ,c=N1,−1D1,−1=22θ−1,minϵ∈0,1minc∈−1,cmaxΛRϵ,c=N1,cmaxD1,cmax>N1,1D1,1=0.This gives the bounds ofℜΛ and ℑΛ for θ∈0.5,1.(2)
For the BDF2 method:Λz=∑j=121−zj/j. In this case, we have(76)Λz=123−4εc+ε22c2−1−iεc˜2−εc.Hence,ℜΛz∈1−cmax2,3/2. For ℑΛz, by noticing(77)ℑΛz=1−c22ε−ε2c≤1−c22−c,forε∈0,1,it is easy to getℑΛz≤2.202.Lemma 6.
Letα,β∈0,1, μ≥0, and ℍ≔s∈ℂ+:ℑs∈π/T,π/Δt. Then, for the linear θ-method with θ∈1/2,1 and the BDF2 method, the arguments ξ±s˜ with s˜=e−sΔt are analytic for s∈ℍ and satisfy ξ+>1 and ξ−<1.Proof.
Fors∈ℍ, from Lemma 5, we have s˜=Λe−sΔt/Δt∈ℂ+. Then, for μ≥0, it holds that ℜas˜>1 for s˜∈ℂ+, where as˜ is defined by (68). This implies ℜ1+a2s˜>1, and thus, the composite functions ξ±s˜ are also analytic. Since ℜ1+a2s˜>1, we have ℜξ+s˜>1⇒ξ+s˜>1; this, together with ξ+s˜ξ−s˜≡1, gives ξ−s˜<1.Now, by using the boundedness conditions at infinity for each subsystem, the constantsAk and Bk in (71) can be determined as Bk=0 and Ak=0 for the first and second subsystems, respectively. We then have E1,mk=Akξ+m and E2,mk=Bkξ−m. By using the transmission conditions in (70), we get the recurrence relation Akηξ+−1=Bk−1ηξ−−1 and Bkη−ξ−=Ak−1η−ξ+, i.e.,(78a)Ak=ρ^dsAk−2,Bk=ρ^dsBk−2,withρ^ds=η−ξ+s˜ηξ+s˜−12,where the subscript “d” denotes “discrete” and(78b)s˜=Λe−sΔtΔt,s∈ℍ≔s∈ℂ+:ℑs∈π/T,π/Δt.Mathematically, we expectρd≪1, which leads to the following min-max problem:(79)minη∈ℝmaxs˜=Λe−sΔt/Δt,s∈ℍη−ξ+s˜ηξ+s˜−12.Letξ+s˜=ξRs˜+iξIs˜. Then, from Lemma 6, we know ξRs˜>1. For η≤1, we have ξ+s˜2−11−η2≥0⇒η−ξRs˜2+ξI2s˜/ηξRs˜−12+ηξI2s˜≥1⇒η−ξ+s˜/ηξ+s˜−12≥1. This implies that, it is sufficient to consider η>1 in (79).Forθ=0.5, the quantity s˜=Λe−sΔt/Δt approaches to the origin, as s approaches to 0,±π in the complex plane. Hence, ξ+ approaches to 1 as s⟶0,±π. This implies that maxs∈ℍρ^ds≡1, independent of η.Forθ∈0.5,1, we know from Lemma 5 that the origin is kept away from s˜, and therefore, it is possible to expect maxs∈ℍρ^ds<1. 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,(80)Dσ¯min,σ¯max,ω¯max≔s∈ℂ+:σ¯min≤ℜs≤σ¯max,0≤ℑs≤ω¯max,which lies in the right-hand side of the vertical line s=1 (i.e., σ¯min>1) 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.Figure 4
An illustration of the region mapped byξ+s˜=as˜+a2s˜−1 (with s˜=Λe−sΔt/Δt and s∈ℍ) and the rectangular region Dσ¯min,σ¯max,ω¯max covering the former. Here, we consider the backward Euler method; for the BDF2 method, the plot looks similar.For the linearθ-method with θ∈0.5,1 and the BDF2 method, if α≥β, by using Lemma 5 and Lemma 7 given below, the boundaries of the rectangular region D, i.e., σ¯min, σ¯max, and ω¯max, can be expressed explicitly, as(81a)σ¯min=σ˜0+σ˜02−1,σ¯max=ξRσ˜1,ω˜max,ω¯max=ξIσ˜0,ω˜max,where ξR,I are defined by (87a). The quantities σ˜0, σ˜1, and ω˜max are defined by(81b)σ˜0=1+Δx22σ01−β+μσ0α−β,σ˜1=1+Δx22σ12+ωmax21−β/2cos1−βψ1+μσ12+ωmax2α−β/2cosα−βψ1,ω˜max=Δx22σ02+ωmax21−β/2cos1−βψ0+μσ02+ωmax2α−β/2cosα−βψ0,withψ0=arccosσ0σ02+ωmax2,ψ1=arccosσ1σ12+ωmax2.The argumentsσ0,1 and ωmax, appearing in (81b), are given by(81c)σ0=1−cmax2θ−1Δt1−θ2+θ2+2θ1−θcmax,σ1=22θ−1Δt,ωmax=12θ−12Δt,θ−method,σ0=1−cmax2,σ1=32,ωmax=2.202,BDF2,where θ∈1/2,1 and cmax=cosπ/T.In general, we have to rely on numerical optimizations to obtain these three quantities, i.e., σ¯min=mins˜=Λe−sΔt/Δt,s∈ℍℜξ+s˜, σ¯max=maxs˜=Λe−sΔt/Δt,s∈ℍℜξ+s˜, and ω¯max=maxs˜=Λe−sΔt/Δt,s∈ℍℜξ+s˜.Lemma 7.
α≥β. Assume 0<β≤α<1 and μ≥0. Let D0=s∈ℂ+:0<σ0≤ℜs≤σ1,0≤ℑs≤ωmax and D=s∈ℂ+:σ¯min≤ℜs≤σ¯max,0≤ℑs≤ω¯max, where σ¯min,max and ω¯max are given by (81a) and (81b). Then, we have ξ+s∈D for s∈D0.Proof.
With the argumentsσ˜0, σ˜1, and ω˜max given by (81b), we define D˜=s∈ℂ+:0<σ˜0≤ℜs≤σ˜1,0≤ℑs≤ω˜max. Then, the proof is divided into two steps: as∈D˜ if s∈D0 and s+s2−1∈D˜ if s∈D˜. The function as is defined by (68).Step 1.
as∈D˜ for s∈D0. Let s=σ+iω∈D0. Then, we can rewrite as as(82)as=1+Δx22eGR1−β,r,σ+μeGRα−β,r,σ+iΔx22eGI1−β,r,σ+μeGIα−β,r,σ,where r=σ2+ω2. For any γ∈0,1, the functions GR,I are defined by (83)GRγ,r,σ=γlogr+logcosγψr,σ,GIγ,r,σ=γlogr+logsinγψr,σ,where ψr,σ=arccosσ/r. Since σ>0, we have ψr∈0,π/2. For fixed σ∈σ0,σ1, because ψr,σ is an increasing function of r, GI is also an increasing function of r. Hence, for γ∈0,1, maxr∈σ,σ2+ωmax2GIγ,r,σ=GIγ,σ2+ωmax2,σ. Define ϑr=arcsinωmax/r. Then, GIγ,σ2+ωmax2,σ = γlogr+logsinγϑr for r=ωmax2+σ2. We have(84)sign∂rγlogr+logsinγϑr=signrsinγϑr−ωmax.
Hence, forγ∈0,1, it holds that ∂rγlogr+logsinγϑr≤0 and this implies that the function GIγ,σ2+ωmax2,σ attains its maximum at σ=σ0. We therefore get the upper bound of ℑas, i.e., ω˜max, as shown in (81b).
It remains to consider the lower and upper bounds ofℜas. For r=σ2+ω2 and σ∈σ0,σ1, we have sign∂rGRγ,r,σ=signrcosγψr−σ. Hence, for γ∈0,1, it holds that ∂rGRγ,r,σ≥0 and this gives(85)minσ≤r≤ωmax2+σ2GRγ,r,σ=γlogσ,maxσ≤r≤ωmax2+σ2GRγ,r,σ=GRγ,ωmax2+σ2,σ.
The first result gives the lower bound ofℜas, i.e., σ˜0. With the function ϑr=arcsinωmax/r, the argument GRγ,ωmax2+σ2,σ can be represented as γlogr+logcosγϑr with r=ωmax2+σ2. We have(86)sign∂rγlogr+logcosγϑr=sign1+ωmaxr2−ωmax2sinγϑrcosγϑr≥0.
Hence,GR attains its global maximum at σ=σ1; this gives σ˜1, the upper bound of ℜas.Step 2.
ξ+s∈D˜ for s∈D˜. Let s=σ+iω with σ∈σ˜0,σ˜1 and ω∈0,ω˜max. Then, ξ+s=ξRσ,ω+iξIσ,ω, where(87a)ξRσ,ω=σ+σ2−ω2−12+4σ2ω2+σ2−ω2−12,ξIσ,ω=ω+σ2−ω2−12+4σ2ω2−σ2−ω2−12.
Clearly, it holds∂ωξIσ,ω>0 and ∂σξRσ,ω>0. Moreover, some simple algebra shows(87b)sign∂ω2ξRσ,ω=signσ2+ω2+1−σ2−ω2−12+4σ2ω2≥0,sign∂σ2ξIσ,ω=signσ2+ω2−1−σ2−ω2−12+4σ2ω2≥0.
Hence,maxω∈0,ω˜maxξIσ,ω=ξIσ,ω˜max and from the second inequality given in (87b), we know that ξIσ,ω˜max attains its maximum at σ=σ˜0; this gives the upper bound of ℑξ+s, i.e., ξIσ˜0,ω˜max. Similarly, we have(88)minσ∈σ˜0,σ˜1minω≤ω˜maxξRσ,ω=ξRσ˜0,0=σ˜0+σ˜02−1,maxσ∈σ˜0,σ˜1maxω≤ω˜maxξRσ,ω=ξRσ˜1,ω˜max.
Combining the analysis in Step1 and Step 2, we complete the whole proof.With the rectangular regionD covering ξ+, we have(89)maxs˜=Λe−sΔt/Δt,s∈ℍη−ξ+s˜ηξ+s˜−12≤maxs∈Dσ¯min,σ¯max,ω¯maxη−sηs−12,and based on this relation, it is reasonable to consider the following equation:(90)minη>1maxs∈Dσ¯min,σ¯max,ω¯maxη−sηs−12,withσ¯min>1.To analyze the solution of this min-max problem, we define(91a)κx,y=Z+Z2−42,withZ=1+x2+y2x,η−=κσ¯max,0,η+=κσ¯min,ω¯max,Tσ,ω,η=σ−η2+ω2ση−12+η2ω2.Then, forηmin<η−≤η+<ηmax, we defined two functions T−η and T+η as(91b)T+η=Tσ¯min,0,η,η∈ηmin,η−,maxTσ¯min,0,η,Tσ¯max,0,η,η∈η−,η+,maxTσ¯min,0,η,Tσ¯max,0,η,Tσ¯min,ω¯max,η,η∈η+,ηmax,T−η=maxTσ¯max,ω¯max,η,Tσ¯max,0,η,Tσ¯min,ω¯max,η,η∈ηmin,η−,maxTσ¯max,ω¯max,η,Tσ¯min,ω¯max,η,η∈η−,η+,Tσ¯max,ω¯max,η,η∈η+,ηmax.Forηmin<η+<η−<ηmax, these two functions are defined by(91c)T+η=Tσ¯min,0,η,η∈ηmin,η+,maxTσ¯min,0,η,Tσ¯min,ω¯max,η,η∈η+,η−,maxTσ¯min,0,η,Tσ¯max,0,η,Tσ¯min,ω¯max,η,η∈η−,ηmax,T−η=maxTσ¯max,ω¯max,η,Tσ¯max,0,η,Tσ¯min,ω¯max,η,η∈ηmin,η+,maxTσ¯max,ω¯max,η,Tσ¯max,0,η,η∈η+,η−,Tσ¯max,ω¯max,η,η∈η−,ηmax.Theorem 2.
Parameter optimization at the discrete level: letL=Δx, α,β∈0,1 and μ≥0. Then, the parameter p involved in the discrete SWR algorithms can be chosen as pd∗=ηopt−1/Δx, where ηopt, the solution of min-max problem (90), is given by(92)ηopt=ηmin,if,T+ηmin≥T−ηmin,ηmax,if,T−ηmax≥T+ηmax,η∗,otherwise,where ηmin=κσ¯min,0, ηmax=κσ¯max,ω¯max, and η∗∈ηmin,ηmax are the unique root of T+η=T−η. With the parameter pd∗, the convergence factor ρd≔maxs∈ℍρ^ds at the discrete level satisfies ρd≤maxT+ηopt,T−ηopt, where ρ^d is given by (78a).Proof.
Lets=σ+iω∈Dσ¯min,σ¯max,ω¯max. Then, it holds that η−s/ηs−12=Tσ,ω,η. For η>1, a routine calculation yields(93)sign∂ηTσ,ω,η=signη−κσ,ω.
Hence, the solution of (90) should satisfy ηopt∈ηmin,ηmax, since otherwise, if ηopt<ηmin (resp. ηopt>ηmax), we can further minimize T by increasing (resp. decreasing) ηopt.
We havesign∂ωTσ,ω,η=sign2ση−1−η2, and therefore, for fixed σ, it holds that maxω∈0,ω¯maxTσ,ω,η=maxTσ,0,η,Tσ,ω¯max,η. Hence, maxs∈Dσ¯min,σ¯max,ω¯maxTσ,ω,η=maxσ∈σ¯min,σ¯maxmaxTσ,0,η,Tσ,ω¯max,η. Another routine calculation yields(94)∂σTσ,ω,η=2η2−1ηση−12+η2ω22σ2−1+η2ησ+1−ω2.
Then, for fixedη>1 and ω≥0, if we vary σ from 1 to +∞, we first meet ∂σTσ,ω,η<0, and then ∂σTσ,ω,η=0, and then ∂σTσ,ω,η>0. Hence, for fixed ω, the function Tσ,ω,η does not have local maximum(s) for σ∈σ¯min,σ¯max. This implies that for given η>1 the global maximum of T on the region D is obtained at the four corners, i.e.,(95)maxs∈DTσ,ω,η=maxTσ¯min,0,η,Tσ¯max,0,η,Tσ¯min,ω¯max,η,Tσ¯max,ω¯max,η.
With the functionsT±η defined by (91b) and (91c), we have(96)minη∈ηmin,ηmaxmaxs∈Dσ¯min,σ¯max,ω¯maxTσ,ω,η=minη∈ηmin,ηmaxmaxT+η,T−η.
It remains to analyze the monotonic properties ofT+η and T−η. From (93), we know that Tσ¯min,0,η and Tσ¯max,ω¯max,η are, respectively, increasing and decreasing for η∈ηmin,ηmax. Similarly, Tσ¯max,0,η is decreasing for η∈ηmin,η− and increasing for η∈η−,ηmax; Tσ¯min,ω¯max,η is decreasing for η∈ηmin,η+ and increasing for η∈η+,ηmax. Hence, the functions T+η and T−η, defined by (91b) and (91c), are, respectively, increasing and decreasing for η∈ηmin,ηmax. This, together with some simple algebra, gives the solution of (96), i.e., ηopt, 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 D, i.e., σ¯min, σ¯max, and ω¯max(2)
Computeηopt according to formula (92) and then let pd∗=ηopt−1/Δx
The quotientΛζ associated with a specified time-integrator determines the boundaries of D, which further determines the parameter ηopt. To make a convergent SWR algorithm, Λe−sΔt should lay in ℂ+ and should be away from the origin, ∀s∈s∈ℂ+:ℑs∈π/T,π/Δt.Remark 4.
Large overlap size: if we use an overlap sizeL=MΔx with M≥2, the following min-max problem occurs:(97)minη>1maxs∈Dσ¯min,σ¯max,ω¯maxTs,η,withTs,η=η−ssM−1ηs−12.
ForM≥2 and η>1, the function Ts,η has local maximums in the interior of D, and therefore, the analysis is much more complicated compared to the case M=1. 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 Figure5, consisting of infinite number of series connected cell-circuits, the state (or modified-nodal-analysis) circuit equation is(98)dbxtdtb=Axt+ft,A=c2⋱⋱⋱1−211−21⋱⋱⋱,Figure 5
(a) The fractional RC circuit in infinite sizeN=∞; (b) a realization of the floating fractional capacitor, the main electronic component of the circuit shown in (a).
(a)(b)wherext=⋯,x−1t,x0t,x1t,⋯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). The 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 xt=ℐtbAxt+f˜t, with ℐtb 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 Uj,mkt by xmt and δx2 by A 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 with(99)s~=s∈H≔s∈C+:Is∈πT,πΔt,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, Theorem 2 is still applicable, but we need to replace ηmax by +∞ there and redefine T±η as(100)T+η=Tσ¯min,0,η,η∈ηmin,η+,maxTσ¯min,0,η,Tσ¯min,ω¯max,η,η>η+,T−η=maxη−2,Tσ¯min,ω¯max,η,η∈ηmin,η+,η−2,η>η+,where η+ and ηmin 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 caseL=Δ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+Δxqopt/2L=1+qopt/2 from the continuous analysis (i.e., Lemma 4) and η=ηopt 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 ρd≔maxs∈ℍρ^d based on three choices of η, where ρ^d and ℍ 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 Theorem 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).Figure 6
Convergence factorρd=maxs∈ℍρ^d based on different parameter η: η=1+Δxqopt with qopt analyzed at the continuous level (see Lemma 4), η=ηopt analyzed at the discrete level (see Theorem 2), and η=ηbest obtained from numerical optimizations (see Remark 2). (a) Λζ=1−ζ (backward Euler); (b) Λζ=∑j=121−ζj/j (BDF2, i.e., the two-step backward difference formula).
(a)(b)From Figure6, 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:(101)ux,t−ℐt1−β∂x2ux,t+μℐt1−αux,t=sinxt,x,t∈0,4×0,18,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 Figure7(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.) 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 maxs∈ℍρ^d=1, 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.Figure 7
(a) Measured convergence rates of the SWR algorithms for three time-integrators, using the parameterη=1.12031 obtained from the continuous analysis. (b) For BDF2 and (c) backward Euler, the error after 4 iterations of the SWR algorithms using various choices of η for the Robin TCs; the parameters analyzed at the continuous and discrete level are denoted by a circle and a star.
(a)(b)(c)At the end of this section, we show in Figure8 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 Δt and T, we use two choices for the parameter η involved in the Robin TCs, the one η=1+Δxqopt 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. 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 η=1+Δxqopt 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]).Figure 8
(a) The effects of the discretization parameters and (b) 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 the BDF2 method, these plots look similar.
(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 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 time-interval.



---

*Source: 5574403-2021-07-22.xml*