Which transmission condition is the best in optimized Schwarz waveform relaxation for viscoelastic equation: Robin or characteristic?

Which transmission condition is the best in optimized Schwarz waveform relaxation for viscoelastic equation: Robin or characteristic? Fu Li1and Yingxiang Xu1* 1School of Mathematics and Statistics, Northeast Normal University, 5268 Renmin Street, Changchun, 130024, JiLin Province, People’s Republic of China. *Corresponding author(s). E-mail(s): [email protected]; Contributing authors: [email protected]; Abstract The viscoelastic equation introduces a viscoelastic damping term in the classical wave equation, thus exhibits properties of parabolic systems accompanied by hyperbolic properties. On the other hand, it is well known that for parabolic systems, when solved using the optimized Schwarz waveform relaxation (OSWR) method, Robin transmission condition is preferred, however characteristic transmission condition is good for hyperbolic system. In this paper, we investigate what transmission condition should be applied according to different strengths of damping. The results show that, in the OSWR algorithm, for weak damping the characteristic transmission condition is preferred, while for strong damping one should use Robin instead, as anticipated. It is also observed that the convergence depends on the relationship between the temporal and spatial mesh sizes applied, which shows overlap may not accelerate the convergence. In addition, we find numerically that a mix of both transmission conditions can only improve the convergence for strong damping, and it remains to be an interesting problem to look for transmission conditions that are more robust in the damping strength. Keywords: Domain decomposition, Optimized Schwarz waveform relaxation, Viscoelastic equation, Transmission condition MSC Classification: 65M55 , 65Y05 1 1 Introduction We focus on the optimized Schwarz waveform relaxation (OSWR) method for the one-dimensional viscoelastic equation ∂ttu(x, t)−ε∂txxu(x, t)−γ∂xxu(x, t) = f(x, t),(x, t)∈Ω×J, u(x, t) = ϕ(x, t),(x, t)∈∂Ω×J, u(x, 0) = ψ1(x), x ∈Ω, ∂tu(x, 0) = ψ2(x), x ∈Ω, (1) where Ω is the space domain, J= (0, T] is the time interval, ε > 0 is the viscoelastic damping coefficient and √γis the wave velocity with γ > 0. The source function f∈L2(0, T;L2(Ω)), the boundary function ϕ(x, t)∈C2(0, T;C2(∂Ω)) and the initial functions ψ1(x)∈H1 0(Ω) and ψ2(x)∈L2(Ω) are given [13]. Compared to the classical wave equation, the viscoelastic equation uses a damping term −ε∂txxuto provide a more accurate model in many applications, for example, the propagation of vibration waves through viscoelastic media [36]. We refer the reader to [48] for more applications in science and engineering. Recently, Gander et al. [23] numerically validated that the viscoelastic damping −ε∂txxuworks much better than the first order damping term −∂tu(which results in the telegrapher’s equation) for modeling the vibration of an elastic string. The well-posedness of the viscoelastic equation has been studied extensively [2,42, 44]. Very recently, there has been significant research focused on the decay properties and the asymptotic behavior of solutions [3,10,30]. However, analytical solutions for many practical PDEs are often elusive, due to complex data (e.g., intricate initial or boundary conditions, sophisticated sources, discontinuous coefficients) and complex domain geometries. Therefore, it is essential to study the numerical methods for solving the viscoelastic equation. In fact, all prevalent techniques for spatial discretization can be applied to the viscoelastic equation (1), see [46] for the finite element methods, [4] for the mixed finite element methods, [32] for the finite difference methods, [35] for the generalized finite difference methods (also known as finite volume element methods), [45] for the discontinuous Galerkin methods, and [47] for the weak Galerkin finite element methods, etc. Recently, meshless methods have also gained attention in solving the viscoelastic equation (1), as demonstrated in [40,41]. However, most numerical methods for solving equation (1) are based on time stepping. For example, the reader can refer to [38,50] for a Crank-Nicolson scheme, where an extrapolation approach and a proper orthogonal decomposition technique are used to reduce the degree of freedom. In fact, there is not much experience in parallel computation for solving the viscoelastic equation. Adey and Brebbia [1] proposed long ago to solve the viscoelastic equation in parallel in the Laplace transformed plane by the finite element method, and then a least square collocation method is applied to construct the whole solution in the time domain. And the authors [34] studied a diagonalization-based parallel-in-time algorithm for Crank-Nicolson’s discretization of the viscoelastic equation, and through rigorous analysis, found that the spectral radius of the iteration matrix is uniformly bounded by υ/(1 −υ), where υ∈(0,1 2) is a free parameter, independent of the model parameters and the discretization parameters. 2 However, except for the above-mentioned works, we do not find any published articles concerning the parallel computing of viscoelastic equations. To speedup the computation, we turn our attention to the waveform relaxation method based on domain decomposition, which is developed for evolution equations [26,22,27]. The domain decomposition method originated from Schwarz’s seminal work [43] in 1870, and was developed in 1990 by Lions [37] as a parallel solver, which finally led to the optimized Schwarz algorithms, cf.[20,14] and references cited therein. Nowadays the optimized Schwarz algorithm has received extensive attention from many scholars and has been successfully used in many stationary engineering models, such as [20,11,52]. The waveform relaxation algorithm was invented to solve the extremely large system of evolution differential equations that arise in circuit simulations [33]. Therefore, the waveform relaxation method has been applied to the solution of many different types of problems, for example, the fractional differential equation [31], the singular perturbation problem [53], etc. The Schwarz waveform relaxation (SWR) algorithm was proposed for space-time problems [26,22], and independently in [27], at the continuous level. The spatial domain is decomposed into subdomains, and time-dependent problems are solved iteratively on subdomains, exchanging information at the interfaces between subdomains. This approach [17] permits a different numerical treatment in both space and time of the subdomain problems, and subdomain information is exchanged less often. Similar to the Schwarz algorithm, the SWR method shows fast convergence when using optimized transmission conditions, as first shown in [18], and then treated in detail in [7,17,39,6,49,9] for the diffusive problems, [16,19] for the wave equation, see also the circuit problems in [15,21], the heterogeneous problems in [28,29], and the models with discontinuous coefficients [8]. In this paper, we study the OSWR algorithm for solving the viscoelastic equation (1), and three different transmission conditions are under consideration, aiming at analyzing the effects of the damping strength on the algorithm’s performance. We designed Robin transmission condition for addressing relatively strong damping strength, and a characteristic condition for relatively weak damping. Targeting on making the transmission condition more robust in the damping strength ε, we designed a hybrid transmission condition by mixing the above two. These three transmission conditions are rigorously optimized in an asymptotic sense, using the idea applied in [25,24]. We consider both cases with and without overlap and derive closed-form solutions for the optimized parameters in an asymptotic sense as well as the estimates of the corresponding convergence factors. Interesting findings include, first, the asymptotic behavior of the OSWR algorithm depends not only on the temporal mesh size τ, but also on δ, when τis related to the spatial mesh size hby τ∼hδ; second, an overlap does not always accelerate the convergence of the OSWR algorithm; lastly, the mixed transmission condition, as the combination of Robin and characteristic conditions, still encounter problems in dealing with problems with weak damp. This paper is organized as follows. In Sect. 2, we present the OSWR algorithm for the viscoelastic equation (1). To achieve fast convergence, we devise three transmission conditions for efficiently transferring information between subdomains. We also derive the convergence factor using Fourier analysis. In Sections 3,4, and 5we 3 analyze the asymptotic convergence results for the Robin, characteristic, and mixed transmission conditions, respectively. By solving the min-max problem for the convergence factor in the Fourier frequency domain in an asymptotic sense, we optimize the transmission parameters involved in the subdomain iterations. In Sect. 6, we illustrate our theoretical results using numerical experiments. We then conclude in Sect. 7. 2 The optimized Schwarz waveform relaxation algorithm To well illustrate the theoretical analysis, we consider only a one-dimensional problem and assume that the defining domain is an infinite domain Ω = Rfor ease of analysis. We note that the case of finite defining domains and high-dimensional problems can also be analyzed similarly, but requires a more intricate analysis. As established in [51], the discrepancy between the convergence factors of the infiniteand finite-domain analyses is bounded by a term that decays exponentially with the subdomain size. Thus, the infinite-domain analysis yields satisfactory results as long as the subdomains are not too small. Furthermore, we consider only the two-subdomain scenario here for the sake of analytical simplicity, and the analysis can be extended to the case of multiple subdomains. We thus suppose that the defining domain Ω is decomposed as Ω = Ω1∪Ω2with Ω1= (−∞, L) and Ω2= (0,+∞), where Lindicates the overlap between subdomains. The classical SWR algorithm for the model problem (1) is given by        ∂ttun 1(x, t)−ε∂txxun 1(x, t)−γ∂xxun 1(x, t) = f(x, t),(x, t)∈Ω1×(0, T], un 1(L, t) = un−1 2(L, t), t ∈(0, T], un 1(x, 0) = ψ1(x), x ∈Ω1, ∂tun 1(x, 0) = ψ2(x), x ∈Ω1,        ∂ttun 2(x, t)−ε∂txxun 2(x, t)−γ∂xxun 2(x, t) = f(x, t),(x, t)∈Ω2×(0, T], un 2(0, t) = un−1 1(0, t), t ∈(0, T], un 2(x, 0) = ψ1(x), x ∈Ω2, ∂tun 2(x, 0) = ψ2(x), x ∈Ω2. (2) By linearity, it suffices to consider only the case f= 0, which corresponds to analyzing the error equation directly. We consider the time Fourier transform ˆu(x, k) = Ft(u)(x, k) := Z∞ −∞ u(x, t)e−itkdt, 4 with all functions under consideration vanishing for t < 0. Taking this Fourier transform of the governing equations in (2) in the time direction gives1 −k2ˆun 1(x, k)−ikε∂xx ˆun 1(x, k)−γ∂xx ˆun 1(x, k) = 0,(x, k)∈Ω1×K, ˆun 1(L, k) = ˆun−1 2(L, k), k ∈K, −k2ˆun 2(x, k)−ikε∂xx ˆun 2(x, k)−γ∂xx ˆun 2(x, k) = 0,(x, k)∈Ω2×K, ˆun 2(0, k) = ˆun−1 1(0, k), k ∈K, (3) where krepresents the Fourier frequency and Kis the set of all admissible frequencies involved in a practical calculation. Due to the infinity assumption on the subdomains, the subdomain solutions to the governing equation in (3) hence read ˆun 1=A(k)eλ(k)x,ˆun 2=B(k)e−λ(k)x,(4) where ±λ(k) are the solutions to the characteristic equation −k2−(ikε +γ)λ2= 0 and have the following expressions λ(k) = s−k2 ikε +γ=a(k) + ib(k), with a(k) = |k|qpγ2+ε2k2−γ √2pγ2+ε2k2, b(k) = kqpγ2+ε2k2+γ √2pγ2+ε2k2 and b(k)2−a(k)2=γk2 γ2+ε2k2,2a(k)b(k) = εk3 γ2+ε2k2. To satisfy the interface conditions in (3), the subdomain solutions (4) are specified as ˆun 1(x, k) = ˆun−1 2(L, k)eλ(k)(x−L),ˆun 2(x, k) = ˆun−1 1(0, k)e−λ(k)x. Inserting these solutions into algorithm (3), we obtain by induction ˆu2n 1(0, k) = ρn cla ˆu0 1(0, k),ˆu2n 2(L, k) = ρn cla ˆu0 1(L, k), where the convergence factor ρcla(k, L) of the classical SWR algorithm is given by ρcla(k, L) := e−2λ(k)L=e−2a(k)L.(5) From this convergence factor, we know that an overlap L > 0 is necessary for the convergence of the overlapping SWR algorithm (2): if there is no overlap, algorithm (2) does not converge. It should be noted that the overlap accelerates convergence: the larger the overlap, the faster the algorithm converges. 1Strictly speaking, one should apply the Laplace transform here. However, a Fourier transform can be applied to provide the same convergence factors in a simpler way, as shown in [5]. 5 To achieve fast convergence, we modify the algorithm as follows by applying new transmission conditions        ∂ttun 1(x, t)−ε∂txxun 1(x, t)−γ∂xxun 1(x, t) = f(x, t),(x, t)∈Ω1×(0, T], ∂xun 1(L, t) + S1un 1(L, t) = ∂xun−1 2(L, t) + S1un−1 2(L, t), t ∈(0, T], un 1(x, 0) = ψ1(x), x ∈Ω1, ∂tun 1(x, 0) = ψ2(x), x ∈Ω1,        ∂ttun 2(x, t)−ε∂txxun 2(x, t)−γ∂xxun 2(x, t) = f(x, t),(x, t)∈Ω2×(0, T], −∂xun 2(0, t) + S2un 2(0, t) = −∂xun−1 1(0, t) + S2un−1 1(0, t), t ∈(0, T], un 2(x, 0) = ψ1(x), x ∈Ω2, ∂tun 2(x, 0) = ψ2(x), x ∈Ω2, (6) where Sj, j = 1,2, are linear operators in time that we will determine subsequently to get the best possible performance of the new Schwarz algorithm. As for the classical SWR method, taking a Fourier transform in time for f= 0, we obtain −k2ˆun 1(x, k)−ikε∂xx ˆun 1(x, k)−γ∂xx ˆun 1(x, k)=0,(x, k)∈Ω1×K, ∂xˆun 1(L, k) + σ1ˆun 1(L, k) = ∂xˆun−1 2(L, k) + σ1ˆun−1 2(L, k), k ∈K, −k2ˆun 2(x, k)−ikε∂xx ˆun 2(x, k)−γ∂xx ˆun 2(x, k)=0,(x, k)∈Ω2×K, −∂xˆun 2(0, k) + σ2ˆun 2(0, k) = −∂xˆun−1 1(0, k) + σ2ˆun−1 1(0, k), k ∈K, where σ1(k) and σ2(k) are Fourier symbols of the operators S1and S2, respectively. A discussion similar to the analysis for the classical SWR reveals that the subdomain solutions now satisfy ˆu2n 1(0, k) = ρn opt ˆu0 1(0, k),ˆu2n 2(L, k) = ρn opt ˆu0 2(L, k),(7) where the new convergence factor ρopt is given by ρopt(k, L, σ1, σ2) :=  −λ+σ1(k) λ+σ1(k)·−λ+σ2(k) λ+σ2(k)·e−2λL. The convergence factor ρopt differs from the classical version ρcla in the term in front of the exponential term, which can be controlled such that the SWR method (6) converges as fast as possible. Choosing the symbols σ1,2as σ1,2=λ(k) = s−k2 ikε +γ,(8) the new convergence factor vanishes identically, ρopt ≡0, and in this case the algorithm (6) will converge in two iterations [12]. We remark here that if the damping coefficient ε= 0, i.e., the model problem (1) degenerates to a wave equation, it holds λ(k) = ik √γ. 6 The optimal symbols σ1=σ2=λ(k) = ik √γthen correspond to local operators that can be implemented easily, and the OSWR algorithm for the wave equation can converge in two iterations, as shown in [19]. In fact, to obtain the transmission operators S1and S2, we need to back-transform the Fourier symbols σ1and σ2from the Fourier frequency domain into the operators in the physical domain S1(un 1) = F−1 t(σ1ˆun 1), S2(un 2) = F−1 t(σ2ˆun 2), where F−1 tdenotes the inverse Fourier transform. The use of the optimal symbols in (8) corresponds to nonlocal operators that incur a heavy computational load in each step of the algorithm, due to the convolution that must be evaluated for interface communication. Obviously, if the symbols σ1,2were polynomials in ik, the corresponding operators S1,2would be local differential operators in t, thereby simplifying the implementation of the information exchange. In what follows, we would like to approximate the optimal Fourier symbols σ1,2by linear functions of ik as follows σapp 1(k) = p1+q1·ik, σapp 2(k) = p2+q2·ik, which correspond to the following interface operators S1=p1+q1∂t, S2=p2+q2∂t. The convergence factor of the OSWR algorithm (6) now reads ρ(k, L, p1, p2, q1, q2) =  λ−p1−q1·ik λ+p1+q1·ik ·λ−p2−q2·ik λ+p2+q2·ik ·e−2λL. The OSWR methods then choose the free parameters pi≥0, qi≥0, i = 1,2, towards the best possible performance, by minimizing the convergence factor over all frequencies relevant to the problem min pi≥0,qi≥0max k∈Kρ(k, L, p1, p2, q1, q2),(9) where we comment that only positive Fourier frequencies should be considered, i.e. K= [kmin, kmax], since ρis symmetric in kdue to the presence of the damping term −ε∂txxu, which essentially differs from the case of typical wave equation where the convergence factor is asymmetric in k. However, the function λ(k) depends on kin a complicated fashion, which makes it very hard to optimize the convergence factor ρopt directly. As a help, we therefore propose, highlighted by the fact λ(k)∼qk 2ε(1 + i) as k→+∞, the following approximation for sufficiently large k > 0 ˜ρ(k, L, p1, p2, q1, q2) =  ˜ λ−p1−q1·ik ˜ λ+p1+q1·ik ·˜ λ−p2−q2·ik ˜ λ+p2+q2·ik ·e−2˜ λL, 7 where ˜ λ(k) = ˜a(1 + i) with ˜a=qk 2ε. The approximation error is shown below. Lemma 1 (Approximation error of the convergence factor).For sufficiently large k > 0, the difference between the convergence factor ρand its approximation ˜ρsatisfies the following estimates: (1) When there is a small overlap L > 0: |ρ(k, L, p1, p2, q1, q2)−˜ρ(k, L, p1, p2, q1, q2)| =γL √2ε3 2√k−γ(q1+q2)L ε2q1q2k+O(L2+L2 √k+L2 k). (2) When there is no overlap (L= 0): |ρ(k, 0, p1, p2, q1, q2)−˜ρ(k, 0, p1, p2, q1, q2)|=γ(q1+q2) √2ε3 2q1q2k3 2 +O(1 k2). Proof. A direct expansion of the difference for large kyields |ρ(k, L, p1, p2, q1, q2)−˜ρ(k, L, p1, p2, q1, q2)| =e−√2k εL√2γL 2ε3 2 1 k1 2−(γ(q1+q2)L ε2q1q2−γ2L2 4ε3)1 k−Ψ1 k3 2 +O(1 k2),(10) where Ψ = √2γ(q1+q2) 2ε3 2q1q2−√2γ(4q2 1+ 4q2 2+ 8q1q2+ 3γq2 1q2 2) 8ε5 2q2 1q2 2 L+√2γ2(q1+q2) 4ε7 2q1q2 L2−√2γ3 24ε9 2 L3. Taking the limit L→0 in (10) directly gives assertion (2). For small L, we substitute the Taylor expansion e−√2k εL= 1 −r2k εL+1 εkL2+O(L3) into (10) and collect terms with respect to kand L, which leads to assertion (1). When a uniform temporal mesh of size τis applied, the largest Fourier frequency kmax can be estimated as kmax =π τ, and the smallest Fourier frequency can be estimated as kmin =π T[17]. When an overlapping domain decomposition is applied, we assume the overlap Lis L=C1hwith hbeing the uniform spatial mesh size. Note that for a particular discretization, which should satisfy a stability or an accuracy constraint, we then assume τ=C2hδ, δ > 0. Lemma 1shows that ˜ρwell approximates the convergence factor ρat high frequencies. This enables a hybrid analytical strategy for the min-max problem (9): for high frequencies, we analyze the approximation ˜ρ directly, while for low frequencies, where the approximation may be less accurate, we analyze the monotonicity of ρitself. Integrating these two parts yields an asymptotic 8 solution to (9) as τ→0, without affecting the final asymptotic results. This methodology, developed in [24,25], has proven effective in tackling complex min-max problems arising in optimized Schwarz methods. In practice, the temporal mesh size τmust be chosen sufficiently small (depending on the problem) to ensure this asymptotic regime is valid. 3 Robin transmission conditions Although the model problem (1) is a wave type equation, the occurrence of the damping term −ε∂txxuleads to parabolic properties [3,30]. This sparks our interest in the performance of the Robin transmission conditions, even if they are not efficient for the hyperbolic equations. Assuming σ1=σ2=p, i.e., taking p1=p2=p, q1=q2= 0, the convergence factor of the OSWR algorithm (6) becomes ρR(k, L, p) =  λ(k)−p λ(k) + p 2 ·e−2λ(k)L=(a(k)−p)2+b(k)2 (a(k) + p)2+b(k)2·e−2a(k)L, and the approximate convergence factor now reads ˜ρR(k, L, p) =  ˜ λ(k)−p ˜ λ(k) + p 2 ·e−2˜ λ(k)L=(˜a(k)−p)2+ ˜a(k)2 (˜a(k) + p)2+ ˜a(k)2·e−2˜a(k)L.(11) The optimized parameter p∗is thus the solution of the min-max problem min p>0max k∈KρR(k, L, p).(12) However, directly solving the min-max problem (12) is highly challenging due to the complexity of the expression for ρRand the difficulty in computing its derivatives. To circumvent this, we adopt an alternative approach: first, we determine the optimized parameter p∗by solving the equi-oscillation equation ρR(kmin, L, p∗) = ρR(˜ k∗, L, p∗),(13) where ˜ k∗is the point at which ρR(k, L, p∗) attains its maximum. We then demonstrate that this p∗indeed solves the original min-max problem (12) by analyzing the behavior of ρRfor parameters different from p∗. Thus, the problem reduces to characterizing the maximum points of ρR. Given that Lemma 1establishes ˜ρRas an accurate high-frequency approximation of ρR, we can achieve this by examining the extremum points of ˜ρR. Lemma 2. For 0<p<√2−1 L, the approximate convergence factor ˜ρR(k, L, p)defined in (11), when regarded as a function in k∈K, reaches its unique local maximum at ¯ k=εp(1 + p1−2pL −p2L2) L, 9 zero, the model problem (1) degenerates to a wave equation, which should not be solved using Robin transmission conditions. We thus consider the parameter setting p1=p2= 0, q1=q2=q, which accounts for a characteristic transmission condition. The corresponding convergence factor can now be reformulated as ρC(k, L, q) = a(k)2+ (b(k)−qk)2 a(k)2+ (b(k) + qk)2e−2a(k)L, thus the free parameter qshould be determined towards fast convergence by minimizing the maximum of the convergence factor ρCover all admissible frequencies min q>0max k∈KρC(k, L, q).(25) Correspondingly, the interface transmission operator now reads q∂t±∂x, which we call the characteristic transmission condition. Note that Gander et al. [19] showed that this transmission condition has good convergence performance for the one-dimensional wave equation. Below for the viscoelastic equation, we find the optimized parameter involved in the characteristic transmission condition in the OSWR algorithm and establish the corresponding convergence results. Theorem 4. The equi-oscillation problem between kmin and ˜ k∗=Ch−θ,0< θ < 2 ρC(kmin, L, p∗) = ρC(˜ k∗, L, p∗) (26) is asymptotically solved for hsmall by q∗=       23 4(εC)−1 4ς−1 2 minhθ 4,if 0< θ < 4 3, (2εC)−1 2ς−1 min(qC2C2 1+ 23 2(εC)1 2ςmin +CC1)h1 3,if θ=4 3, 21 2ε−1 2C1 2C1ς−1 minh1−θ 2,if 4 3< θ < 2, (27) where amin =a(kmin),bmin =b(kmin),ςmin = 4kminbmin/(a2 min +b2 min)and L= C1h. Proof. The proof is similar to that of Theorem 1, and is provided in the appendix. For further analysis, we need the following approximate convergence factor ˜ρC(k, L, q) = ˜a(k)2+ (˜a(k)−qk)2 ˜a(k)2+ (˜a(k) + qk)2e−2˜a(k)L.(28) Lemma 3. For q > (1+√2)L ε, as a function of k∈Kthe approximate convergence factor ˜ρC(k, L, q)defined in (28)reaches its unique local maximum at ¯ k=εq +pε2q2−2εqL −L2 εq2L, 16 and reaches its unique local minimum at k=εq −pε2q2−2εqL −L2 εq2L. Proof. From the derivative of ˜ρC(k, h, p) in kwe get ¯ kand kby direct calculation. Theorem 5 (Characteristic transmission condition, overlapping case).Let L=C1h and τ=C2hδ. The min-max problem (25)is asymptotically solved for hsufficiently small by q∗=       24 3ε−1 3ς−2 3 minC 1 3 1h1 3,if δ > 4 3or δ=4 3and C1≥Cc, ˜ CqC 1 3 1h1 3,if δ=4 3and C1< Cc, 23 4(επ)−1 4C 1 4 2ς−1 2 minhδ 4,if δ < 4 3, (29) where Cc= 2−1 4(C2/π)3 4ε1 4ς 1 2 min, and ςmin = 4kminbmin/(a2 min +b2 min),˜ Cq= (2επC2)−1 2C−1 3 1ς−1 min(qπ2C2 1+ (2C2)3 2(επ)1 2ςmin +πC1). The maximum of the corresponding convergence factor ρChas the following asymptotic expansion for hsmall max k∈KρC(k, L, q∗) =       1−24 3ε−1 3ς 1 3 minC 1 3 1h1 3+O(h2 3),if δ > 4 3or δ=4 3and C1≥Cc, 1−ςmin ˜ CqC 1 3 1h1 3+O(h2 3),if δ=4 3and C1< Cc, 1−23 4(επ)−1 4C 1 4 2ς 1 2 minhδ 4+O(hδ 2),if δ < 4 3. Proof. The proof is very similar to that of Theorem 2and is relegated to the appendix. The following results on determining the optimized transmission parameters when the characteristic transmission conditions are applied to a non-overlapping OSWR method can be proven in a fashion similar to that of Theorem 3, so we also omit the details. Theorem 6 (Characteristic transmission conditions, non-overlapping case).For the non-overlapping case, L= 0, the parameter q∗given by q∗= 23 4ε−1 4ς−1 2 mink−1 4 max (30) solves the min-max problem (25)asymptotically for kmax large and the corresponding convergence factor satisfies the following estimate for kmax large enough max k∈KρC(k, 0, q∗)=1−23 4ε−1 4ς 1 2 mink−1 4 max +O(k−1 2 max) where ςmin = 4kminbmin/(a2 min +b2 min). 17 Proof. The proof of this theorem is deferred to the appendix. Remark 3. Assume the time step τis linked to the spatial mesh size hby the relation τ=C2hδ. Comparing Theorem 5with Theorem 6, we find that if δ < 4 3, an overlap cannot improve the performance of the OSWR algorithm with characteristic transmission conditions, which is similar to the result in Remark 2. 5 Mixed transmission conditions Previous analyses show that for relatively large εRobin condition is good, while for relatively small εcharacteristic condition is preferable. To incorporate the features of the two transmission conditions mentioned above, we hope to design a unified form of transmission conditions that can deal with arbitrary εwell. We thus choose p1=p2=p, q1=q2=q, the convergence factor simplifies to ρM(k, L, p, q) = (a(k)−p)2+ (b(k)−qk)2 (a(k) + p)2+ (b(k) + qk)2e−2a(k)L. The corresponding interface transmission operators are S1=S2=p+q∂t, which we call the mixed transmission condition since it is a mix of Robin and characteristic conditions. The free parameters pand qthen should be determined by solving the min-max problem min p>0,q>0max k∈KρM(k, L, p, q).(31) However, it is not an easy task, and we need help from the following approximate convergence factor ˜ρM(k, L, p, q) = (˜a(k)−p)2+ (˜a(k)−qk)2 (˜a(k) + p)2+ (˜a(k) + qk)2e−2˜a(k)L,(32) which has the properties below. Lemma 4. As a function of k, the approximate convergence factor ˜ρM(k, L, p, q) defined in (32)has four critical points for k > 0at which it attains local extrema: two local minima at k=k1and k=k2, and two local maxima at k=¯ k1and k=¯ k2. Their asymptotic behaviors for small Lare given by: k1∼εp2,¯ k1∼p q, k2∼1 εq2,¯ k2∼2 Lq .(33) Proof. The change of variables ˜ k=r2k εL, ˜p= 2Lp, ˜q=εq L 18 reduces the approximate convergence factor to ˜ρM(˜ k, ˜p, ˜q) = (˜ k−˜p)2+ (˜ k−˜q˜ k2)2 (˜ k+ ˜p)2+ (˜ k+ ˜q˜ k2)2e−˜ k. By computing the derivative with respect to ˜ kand setting it to zero, we find that the extreme points correspond to the roots of the polynomial P(˜ k) = ˜q4˜ k8−4˜q3˜ k6+(2˜p2˜q2+(−12˜q2−8˜q)˜p+ 8˜q+4)˜ k4+(12˜p2˜q−8˜p)˜ k2+ ˜p3(˜p+4). Following the approach in Theorem 4.12 in [7], we introduce X= ˜q˜ k2, which transforms P(˜ k) into P(X) = X4−4X3+ (2˜p2−12˜p−8˜p ˜q+8 ˜q+4 ˜q2)X2+ (12˜p2−8˜p ˜q)X+ ˜p3(˜p+ 4). For sufficiently large ˜q, small ˜p˜q, and small L, the dominant part of P(X) is P0(X) = X4−4X3+8 ˜qX2−8˜p ˜qX+ 4˜p3, which admits four distinct positive roots Xj(j= 1,...,4) with the asymptotic behavior: X1∼˜p2˜q 2, X2∼˜p, X3∼2 ˜q, X4∼4. A perturbation argument then implies that, for sufficiently small L, the polynomial P(˜ k) also possesses four distinct positive roots. Recalling that ˜ k=qX ˜q, we obtain their asymptotic behavior as ˜ k1∼˜p √2,˜ k2∼s˜p ˜q,˜ k3∼√2 ˜q,˜ k4∼2 √˜q, Reverting to the original variables yields the stated extreme points stated in (33). Similar to the Robin transmission conditions where the min-max problem (11) is solved by a two-point equi-oscillation, numerical investigation shows that the min-max problem (31) is solved by a three-point equi-oscillation, as illustrated in Fig. 1. We thus first determine the optimized parameters p∗and q∗by solving the three-point equi-oscillation problem ρM(kmin, L, p∗, q∗) = ρM(˜ k∗ 1, L, p∗, q∗) = ρM(˜ k∗ 2, L, p∗, q∗),(34) where ˜ k∗ 1and ˜ k∗ 2are the two local maximum points of ρM, then show the optimality of the resulting parameters. 19 100101102103104105 0.1 0.15 0.2 0.25 0.3 p q 1001021041061081010 0 0.1 0.2 0.3 0.4 0.5 p q Fig. 1: Illustration of the equi-oscillation property. Left: τ=h(kmax <˜ k∗ 2); Right: τ=h2(kmax >˜ k∗ 2). Parameters: ε= 2, γ = 1, L = 2h, h = 10−4. Theorem 7 (Asymptotic results for three-point equi-oscillation).The three-point equi-oscillation problem (34)between kmin,˜ k∗ 1=B1h−θ1and ˜ k∗ 2=B2h−θ2for 0< θ1< θ2≤8 5, is asymptotically solved for hsmall with L=C1hby            p∗= 21 4ε−1 4B 1 4 1a 1 2 minh−θ1 4, q∗= 2−1 4ε−3 4B 1 4 1B−1 2 2a−1 2 minhθ2 4−θ1 4if θ2≤8 5, θ1<θ2 2, p∗=√2aminB−1 4 1B 1 4 2h−θ2−θ1 4, q∗=ε−1 2(B1B2)−1 4hθ2+θ1 4if θ2≤8 5, θ1>θ2 2, p∗=B 1 4 1B 1 4 2aminζminh−θ1 4, q∗= 21 4ε−1 2B 1 4 1B−1 4 2ζminh3θ1 4,if θ2<8 5, θ1=θ2 2, p∗=˜ Cp(C1h)−1 5, q∗=˜ Cq(C1h)3 5,if θ2=8 5, θ1=4 5, where ζmin := (21 2B1+ 2ε1 2B 1 2 2amin)−1 2,˜ Cp:= B 1 2 1(2 + B2C 8 5 1˜ Cq−2ε(B1B2)1 2C 6 5 1˜ C2 q)/ (2εB 1 2 2C 2 5 1˜ Cq)and ˜ Cqis the unique positive root of the polynomial P(ξ) =2√2εC 14 5 1B 3 2 1B 3 2 2ξ3−4√2C 8 5 1B1B2ξ−4√2B1 +C 6 5 1(4√2εB 3 2 1B 1 2 2−√2C2 1B1B2 2+ 8ε3 2B 1 2 1B2amin)ξ2. (35) Proof. We are still seeking a solution in an asymptotic sense. We make the ansatz p∗=CpL−αand q∗=CqLβ, since numerical tests show that p∗tends to infinity and q∗tends to zero as L→0. Noting that L=C1h, we easily find asymptotic expansions of the three terms in (34), ρM(kmin, L, p∗, q∗) =1 −4amin Cp (C1h)α+O(hmin{2α,1}),(36) ρM(˜ k∗ 1, L, p∗, q∗) =˜ρM(˜ k∗ 1, L, p∗, q∗) + O(L1+ θ1 2−β) =1 −2√2εCp √B1Cα 1 hθ1 2−α−2p2εB1CqCβ 1hβ−θ1 2(37) +O(hmin{1−θ1 2,θ1−2α,2β−θ1}), 20 ρM(˜ k∗ 2, L, p∗, q∗) =˜ρM(˜ k∗ 2, L, p∗, q∗) + O(L1+ θ2 2−β) =1 −2√2 √εB2CqCβ 1 hθ2 2−β−√2B2C1 √εh1−θ2 2(38) +O(hmin{θ2−2β,1−β,2−θ2}). We need to consider the relationship between the exponents θ1 2−αand β−θ1 2, as well as θ2 2−βand 1 −θ2 2. Equalizing the exponents, i.e. solving α=θ1 2−α=β−θ1 2=θ2 2−β= 1 −θ2 2, one finds α∗=1 5, β∗=3 5, θ∗ 1=4 5, θ∗ 2=8 5. Then we have the following situations for 0< θ1< θ2≤8 5: 1. θ2≤8 5and θ1<θ2 2. Equalizing the dominant terms, that is solving 4amin Cp (C1h)α=2√2εCp √B1Cα 1 hθ1 2−α=2√2 √εB2CqCβ 1 hθ2 2−β, we get Cp= 21 4ε−1 4B 1 4 1a 1 2 minCα 1,Cq= 2−1 4ε−3 4B 1 4 1B−1 2 2a−1 2 minC−β 1, and α=θ1 4,β= θ2 2−θ1 4. 2. θ2≤8 5and θ1>θ2 2. Now, matching the exponents of the dominant terms 4amin Cp (C1h)α= 2p2εB1CqCβ 1hβ−θ1 2=2√2 √εB2CqCβ 1 hθ2 2−β, we arrive at Cp=√2aminB−1 4 1B 1 4 2Cα 1,Cq=ε−1 2(B1B2)−1 4C−β 1, and α=θ2−θ1 4, β=θ2+θ1 4. 3. θ2<8 5and θ1=θ2 2. Now, via solving 4amin Cp (C1h)α=2√2εCp √B1Cα 1 hθ1 2−α+ 2p2εB1CqCβ 1hβ−θ1 2=2√2 √εB2CqCβ 1 hθ2 2−β, we obtain α=θ1 4, β =3θ1 4,Cp=B 1 4 1B 1 4 2amin(21 2B1+ 2ε1 2B 1 2 2amin)−1 2Cα 1and Cq= 21 4ε−1 2B 1 4 1B−1 4 2(21 2B1+ 2ε1 2B 1 2 2amin)−1 2C−β 1 4. θ2=8 5, θ1=θ2 2=4 5. In this case, we need to solve 4amin Cp (C1h)α=2√2εCp √B1Cα 1 hθ1 2−α+ 2p2εB1CqCβ 1hβ−θ1 2 =2√2 √εB2CqCβ 1 hθ2 2−β+√2B2C1 √εh1−θ2 2, 21 to obtain α=1 5, β =3 5,Cp=B 1 2 1(2 + B2C 8 5 1Cq−2ε(B1B2)1 2C 6 5 1C2 q)/(2εB 1 2 2C 2 5 1Cq), Cqis the unique positive root of the polynomial (35). Theorem 8 (Mixed transmission condition, overlapping case).Let L=C1hand τ=C2hδ. The min-max problem (31)is asymptotically solved for hsufficiently small by (p∗= 2−1 5a 4 5 minC−1 5 1h−1 5, q∗= 2−2 5ε−1a−2 5 minC 3 5 1h3 5,if δ > 8 5,or δ=8 5and C1≥Cm,(39) p∗=CpqC−1 1h−1 5, q∗=Cpq ˜ C−2 pq C−1 1h3 5,if δ=8 5and C1< Cm,(40) (p∗= 2−1 8π1 8(εC2)−1 8a 3 4 minh−δ 8, q∗= 2−5 8π−3 8ε−5 8C 3 8 2a−1 4 minh3δ 8,if 0< δ < 8 5,(41) where Cm= 27 8π−5 8ε5 8C 5 8 2a 1 4 min and Cpq := 2((2πεC2)1 2amin −C2˜ C2 pq)π−1with ˜ Cpq being the unique positive root of polynomial P(ζ)=4√2ε1 2C2 2ζ4−16επ 1 2C 3 2 2aminζ2−π2C2 1aminζ+ 8√2ε3 2πC2a2 min. The maximum of the corresponding convergence factor ρMhas the following asymptotic expansion for hsmall max k∈KρM(k, L, p∗, q∗) =     1−4·21 5a 1 5 minC 1 5 1h1 5+O(h2 5),if δ > 8 5,or δ=8 5and C1≥Cm, 1−4aminC−1 pq C1h1 5+O(h2 5),if δ=8 5and C1< Cm, 1−4·21 8(εC2)1 8π−1 8a 1 4 minhδ 8+O(hδ 4),if 0< δ < 8 5. Proof. Similarly, we first, based on Theorem 7, determine the solutions p∗and q∗ to the equi-oscillation problem in the admissible frequency domain K. We make the ansatz: p∗=CpL−α, q∗=CqLβ. From Lemmas 1and 4, we know that the location of the interior maximum of ρM(k, L, p, q) in the asymptotic sense is ¯ k∗ 1∼p∗ q∗=CpC−1 q(C1h)−(α+β),¯ k∗ 2∼2 Lq∗= 2C−1 q(C1h)−(1+β). Note that, in a numerical calculation, an estimate for the maximum frequency is kmax =π C2hδ. We thus need to consider the relationship between the maximum frequency kmax and ¯ k∗ 2. If kmax ≥¯ k∗ 2, which requires that δ > 1+βor δ= 1+βand C1≥(2C2/(πCq)) 1 1+β, we should equi-oscillate between kmin,¯ k∗ 1and ¯ k∗ 2. Thus, let ˜ k∗ 1=¯ k∗ 1,˜ k∗ 2=¯ k∗ 2, i.e., 22 θ1=α+β, θ2= 1 + βin Theorem 7, equations (37), (38) now read ρM(¯ k∗ 1, L, p∗, q∗) = ˜ρM(¯ k∗ 1, L, p∗, q∗) + O(L1+ α−β 2) = 1 −4√2ε1 2(CpCq)1 2C β−α 2 1hβ−α 2+O(hβ−α),(42) ρM(¯ k∗ 2, L, p∗, q∗) = ˜ρM(¯ k∗ 2, L, p∗, q∗) + O(L3+β 2) = 1 −4ε−1 2C−1 2 qC 1−β 2 1h1−β 2+O(h1−β),(43) which lead to α∗=1 5, β∗=3 5, Cp= 2−1 5a 4 5 min, Cq= 2−2 5ε−1a−2 5 min.(44) It follows from these results that kmax ≥¯ k∗ 2is equivalent to δ > 8 5, or δ=8 5and C1≥27 8π−5 8ε5 8C 5 8 2a 1 4 min =: Cm. Then, when δ=8 5and C1< Cm, or 0 < δ < 8 5, we have kmax <¯ k∗ 2, and we should equi-oscillate between kmin,¯ k∗ 1and kmax. Choosing ˜ k∗ 2=kmax =π C2hδin Theorem 7, one finds (α∗=1 5, C∗ p=CpqC−4 5 1, β∗=3 5, C∗ q=Cpq ˜ C−2 pq C−8 5 1,if δ=8 5and C1< Cm,(45) (α∗=δ 8, C∗ p= 2−1 8π1 8(εC2)−1 8a 3 4 minC δ 8 1, β∗=3δ 8, C∗ q= 2−5 8π−3 8ε−5 8C 3 8 2a−1 4 minC−3δ 8 1,if 0 < δ < 8 5,(46) where ˜ Cpq and Cpq are defined below (41). Combining equations (44), (45) and (46), one finds the expression of the optimized parameters as shown in Theorem 8. Inserting the asymptotic values of p∗, q∗given in (39), (40), (41), and k=kmin into the convergence factor ρM(k, L, p, q) and expanding for hsmall, we obtain the asymptotic expansions in Theorem 8. We thus only need to show that p∗,q∗given in (39), (40) and (41) solve the minmax problem (31) asymptotically. Here we only show the detailed proof for the first case, kmax >¯ k∗ 2(c.f. the right plot in Fig. 1), and the other two cases (c.f. the left plot in Fig. 1) are similar. We first show that for kmin < k < k1the convergence factor ρM is a decreasing function in k, where k1=ε(p∗)2∼2−2 5εa 8 5 minL2 5is defined in Lemma 4. Again, we analyze using an asymptotic analysis as the following. Assume δkis a small increment of the frequency k, we then have by a direct calculation ρM(k+δk,0, p∗, q∗) = ρM(k, 0, p∗, q∗) 23 +                  −√2(√ε2C2 k+γ2(ε2C2 k+2γ2)−2γ3) C∗ p(ε2C2 k+γ2)3 2q√ε2C2 k+γ2−γL1 5δk+O(L2 5δk),if ϑ= 0, −√2 √εCkC∗ p L1 5+ϑ 2δk+O(L2 5δk),if 0 <ϑ<2 5, √2εC∗ p(Ck−ε(C∗ p)2) √Ck(√2εCkC∗ p+ε(C∗ p)2+Ck)2L2 5δk+O(L4 5δk) if ϑ=2 5, √2εC∗ pC−3 2 kL3ϑ 2−1 5δk+O(Lmin{2(ϑ−1 5),ϑ 2+3 5}δk) if 2 5<ϑ<4 5, (47) with k=CkL−ϑ. Obviously, equation (47) shows that in asymptotic sense ρM(k, 0, p∗, q∗) is monotonically decreasing in kfor k∈(kmin, k1). Thus, ρM(k, L, p∗, q∗), as a product of two decreasing positive functions ρM(k, 0, p∗, q∗) and ρcla(k, L), must decrease in kfor k < k1. As a consequence, kmin is a possible maximum point of ρM(k, L, p∗, q∗) for hsmall enough, as well as the two interior maxima points ¯ k∗ 1,2found in Lemma 4. We now need to show that (p, q) differs from (p∗, q∗) leading to a convergence factor ρMthat could be further improved. To this end, let p=CpL−α, q =CqLβand p∗=CpL−α∗, q∗=CqLβ∗, where α∗=1 5, β∗=3 5, as shown in (39). We now need to show that if (p, q)= (p∗, q∗), then there exists a k∗such that for h > 0 small enough ρM(k∗, L, p, q)>max kρM(k, L, p∗, q∗) = ρM(kmin, L, p∗, q∗), where the search for k∗can be steered by tracking how ρMevolves with pand q, as illustrated in Fig. 1. Here, we again employ the technique proposed by Gander and Xu [24,25]. We first show that if (α, β)= (α∗, β∗), then the asymptotic order of the resulting OSWR method will be enlarged. To do this, it is sufficient to treat the following cases: (a) α > α∗, β > −α. At k∗=kmin, we have ρM(k∗, L, p, q)=1−4amin Cp Cα 1hα−4kminbminCq C2 p Cβ+2α 1hβ+2α −2aminC1h+o(hmin{α,β+2α,1}). (b) α > α∗, β =−α. At k∗=kmin, we have ρM(k∗, L, p, q) = 1 −2aminC1h−4aminCp+ 4kminbminCq C2 p+k2 minC2 q Cα 1hα+o(hmin{α,1}). (c) α > α∗, β < −α. At k∗=kmin, we have ρM(k∗, L, p, q)=1−4bmin kminCq (C1h)−β−4aminCp k2 minC2 q (C1h)−(2β+α) −2aminC1h+o(hmin{−β,−(2β+α),1}). 24 (d) α=α∗, β > β∗. Here we have to treat two cases: If β < 1, for k∗=CkL−(α+β), we have ρM(k∗, L, p, q)=1−2√2ε(CqCk+Cp) √Ck (C1h)β−α 2+o(hβ−α 2). If β≥1, we take k∗=CkL−1to obtain ρM(k∗, L, p, q)=1−2√2εCp √Ck (C1h)3 10 +o(h3 10 ). (e) α≤α∗, β < β∗. At k∗=CkL−(1+β), we have ρM(k∗, L, p, q) = 1 −√2(2 + CkCq) √εCkCq (C1h)1−β 2+O(h1−β 2). (f) α < α∗, β > β∗. At k∗=CkL−4 5, it holds ρM(k∗, L, p, q)=1−2√2εCp √Ck (C1h)2 5−α−2p2εCkCq(C1h)β−2 5+o(hmin{2 5−α,β−2 5,3 5}). (g) α < α∗, β =β∗. At k∗=CkL−(α+β), we get ρM(k∗, L, p, q)=1−2√2ε √Ck (Cp+CkCq)(C1h)β−α 2+o(hβ−α 2). We therefore see that in each case above, at the given k∗, the convergence factor ρM(k∗, L, p, q) behaves asymptotically like 1 −Chδwith δ > 1 5. It remains thus to consider the case α=α∗and β=β∗but (Cp, Cq)= (C∗ p, C∗ q). From (36) we see ρM(kmin, L, p, q)> ρM(kmin, L, p∗, q∗) asymptotically if Cp> C∗ p. From (43) we see ρM(¯ k2, L, p, q)> ρM(¯ k∗ 2, L, p∗, q∗) if Cq> C∗ q, and from (42) we conclude that ρM(¯ k1, L, p, q)> ρM(¯ k∗ 1, L, p∗, q∗) if Cp< C∗ por Cq< C∗ q, which concludes the proof. For the non-overlapping case, the mixed transmission condition can also be optimized, the results are shown below. Theorem 9 (Mixed transmission condition, non-overlapping case).When there is no overlap, i.e. L= 0, the min-max problem (31)is solved by the parameters p∗= (2ε)−1 8a 3 4 mink 1 8 max, q∗= (2ε)−5 8a−1 4 mink−3 8 max in an asymptotic sense of kmax large, which leads the convergence factor to the following asymptotic estimate max k≥kmin ρM(k, 0, p∗, q∗) = 1 −4·(2ε)1 8a 1 4 mink−1 8 max +O(k−1 4 max). Proof. The proof is very similar to that for Theorems 3and 8. 25 Table 3: Number of iterations (Niter) required by the OSWR algorithm (6) with Robin, characteristic, and mixed transmission conditions compared with their predictions by the convergence factors (Npred). Parameters: τ=h,tol = 10−3. Robin Char. Mixed hmax ρRNpred Niter max ρCNpred Niter max ρMNpred Niter 1 40 0.3432 13 8 0.2996 12 8 0.1211 7 5 1 80 0.4087 16 10 0.3685 14 9 0.1505 8 5 1 160 0.4720 19 11 0.4349 17 11 0.1769 8 5 1 320 0.5317 22 13 0.4980 20 13 0.2016 9 6 1 640 0.5874 26 15 0.5571 24 15 0.2255 10 6 1 1280 0.6388 31 18 0.6117 29 17 0.2498 10 6 6.2 The effect of damping coefficient ε In this section, we investigate the robustness of various transmission conditions to the damping coefficient ε. Using the same exact solution as in Sect. 6.1, we vary the value of the damping coefficient εand see how the OSWR algorithm responds. We note that when εapproaches 0, the viscoelastic equation (1) will degenerate into a typical wave equation, on which Gander et al. [19] showed that the best transmission condition is the characteristic condition with parameter q= 1/√γ. In the experiment, we also examine the convergence effect of this parameter for comparison. To see the effect of the damping coefficient εmore clearly, we show the number of iterations as a function of εin the left plot of Fig. 7, with tol=10−3and τ=h= 0.01. Here, a fixed tolerance for varying εis still good, since the solutions do not vary with ε, and the maximum norm is applied. It is clear that for very small ε, the characteristic transmission condition outperforms the Robin and mixed, conforming well with the analysis from Remark 5, though it remains slightly inferior to the parameter selection q= 1/√γdetermined directly by the wave equation (we call it wave transmission condition in the rest of this paper). For not very small ε, say for example, 0.1<ε<10, the Robin and characteristic transmission conditions perform very similarly, but neither converges faster than the mixed transmission condition, which is consistent with our asymptotic analysis results. For even strong damping, say ε > 10, we find that for all transmission conditions considered in this paper the performance deteriorates. However, the mixed transmission condition is still the best, as indicated by Remark 5. Among the transmission conditions using only one parameter, Robin performs the best, which is again consistent with Remark 5. As a comparison, we repeat this experiment with a smaller mesh size, τ=h= 0.002, and show the results in the right plot of Fig. 7. The results show that each curve has a similar trend as εvaries, compared with the case τ=h= 0.01. However, for small ε, both the Robin and mixed 32 transmission conditions require significantly more iterations to converge. This observation is consistent with the asymptotic estimates in Table 1, which indicate that for small ε, reducing the mesh size hdrives the convergence factors closer to 1. However, the mixed transmission condition holds promise as it encompasses the characteristic condition–known for its excellent performance–as a special case. The key to realizing its potential under weak damping lies in developing a new strategy for determining the optimal parameter. 10-4 10-3 10-2 10-1 100101102 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 Iterations Robin Characteristic Mixed Wave 10-4 10-3 10-2 10-1 100101102 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 Iterations Robin Characteristic Mixed Wave Fig. 7: The number of iterations as a function of εwith L= 2h, τ =h, tol = 10−3. Left: h= 0.01; Right: h= 0.002. To evaluate the predictive capability of asymptotic analysis for determining optimal parameters in the case of weak damping, we choose ε= 0.01. Following the experimental methodology used for generating Fig. 4, we perform the numerical experiments with tol = 10−3, and two step sizes: τ=h= 0.01 and τ=h= 0.002. The results presented in Fig. 8show that, for both the Robin and characteristic conditions, a very fine mesh is required to achieve an accurate prediction based on asymptotic analysis. However, it appears that for the characteristic condition, the damping coefficient ε and the mesh size hmust be related in some way to achieve an accurate prediction. 6.3 Finite element discretization To evaluate the impact of spatio-temporal discretizations on the performance of the OSWR algorithm (6), we discretize space using linear finite elements and time using the backward Euler scheme. We then reproduce the two experiments from subsection 6.1 with this new discretization while keeping all parameter settings unchanged. The results, presented in Figures 9and 10, are largely consistent with those in Figures 4and 6, confirming that our analysis is only marginally affected by the choice of discretization. 33 0 100 200 300 400 p 20 30 40 50 60 70 Iterations 0 0.5 1 1.5 2 2.5 3 3.5 4 q 4 6 8 10 12 14 16 18 20 22 Iterations 6 7 7 8 8 8 8 8 9 9 9 9 9 9 10 10 10 10 10 10 11 11 11 11 11 11 12 12 12 12 12 12 13 13 13 13 13 13 14 14 14 14 14 14 15 15 15 15 15 15 16 16 16 16 16 16 1717 17 17 17 17 17 18 18 18 18 19 19 19 20 20 20 20 22 22 22 22 24 24 24 24 26 26 26 28 28 30 0.5 1 1.5 2 2.5 3 3.5 4 4.5 p 0.5 1 1.5 2 2.5 3 3.5 4 4.5 q (a) Parameters: ε= 0.01, τ =h, h = 0.01, tol = 10−3 0 100 200 300 400 p 20 40 60 80 100 120 140 160 180 Iterations 0 0.5 1 1.5 2 2.5 3 3.5 4 q 4 6 8 10 12 14 16 18 20 22 Iterations 7 7 8 8 8 8 9 9 9 9 9 9 10 10 10 10 10 10 11 11 11 11 11 11 12 12 12 12 12 12 13 13 13 13 13 13 14 14 14 14 14 14 15 15 15 15 15 15 15 16 16 16 16 16 16 16 17 17 17 17 17 17 17 1818 18 18 18 18 18 19 19 19 19 20 20 20 20 22 22 22 22 24 24 24 24 26 26 26 28 28 28 30 30 0.5 1 1.5 2 2.5 3 3.5 4 4.5 p 0.5 1 1.5 2 2.5 3 3.5 4 4.5 q (b) Parameters: ε= 0.01, τ =h, h = 0.002, tol = 10−3 Fig. 8: Optimized parameter (∗) found by asymptotic analysis, compared with the performance of other parameter values. Left: Robin; Middle: Characteristic; Right: Mixed. 0 2 4 6 8 10 12 p 5 10 15 20 25 30 35 40 45 Iterations 0 0.5 1 1.5 q 5 10 15 20 25 30 35 Iterations 6 7 7 7 8 8 8 8 9 9 9 9 10 10 10 10 10 11 11 11 11 11 12 12 12 12 13 13 13 13 14 14 14 14 14 15 15 15 15 15 16 16 16 16 16 17 17 17 17 17 18 18 18 18 18 19 19 19 19 19 20 20 22 0.5 1 1.5 2 2.5 3 3.5 4 4.5 p 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 q Fig. 9: Optimized parameter (∗) found by asymptotic analysis, compared with the performance of other parameter values when the viscoelastic equation is discretized by the finite element method. Left: Robin; Middle: Characteristic; Right: Mixed. Parameters: ε= 2, τ =h, h = 0.01, tol = 10−3. 7 Conclusion In this paper, the OSWR algorithm is developed for solving the viscoelastic equation. With the aid of an asymptotically accurate approximation of the convergence factor, we derive closed-form asymptotic formulas for Robin, characteristic and mixed transmission conditions, and obtain the corresponding convergence rate estimate, for both overlapping and non-overlapping cases. The theoretical results reveal for the first time that the efficacy of subdomain overlap in accelerating the convergence of the 34 10-2 10-1 100 101 102 103 104 Iterations Robin Characteristic Mixed Classical h-1/4 h-1/8 h-1 10-2 10-1 100 101 102 103 104 Iterations Robin Characteristic Mixed Classical h-1/3 h-1/5 h-1 Fig. 10: When discretized by the finite element method, the asymptotic behavior as the spatial grid is refined with the overlap L= 2h, together with the predicted rates from analysis, both for the classical SWR and OSWR with Robin, characteristic, and mixed transmission conditions. Left: τ=h; Right: τ=h2. OSWR algorithm for the viscoelastic wave equation is fundamentally governed by the value of δin τ=C2hδ, i.e., the scaling relationship between temporal and spatial discretizations. Similar results have been reported for parabolic problems in [17,7]. More importantly, though Robin and characteristic transmission conditions behave asymptotically the same, they are affected by the viscoelastic coefficient εdifferently: the characteristic transmission condition performs quite well across all values of ε, while the Robin condition outperforms the characteristic one only for very large ε. Furthermore, although the mixed transmission condition is theoretically expected to perform well under all circumstances, numerical results reveal performance degradation in weak damping scenarios when the grid is not sufficiently refined. It is therefore promising to investigate the optimal parameters for the mixed transmission condition under weak damping for mesh sizes commonly encountered in practice. Statements and Declarations •Funding: This work was supported by the National Key R&D Program of China (No. 2020YFA0714102), and the Scientific Research Project of Education Department of Jilin Province (No. JJKH20250297BS). •Conflict of interest/Competing interests: The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. •Ethics approval and consent to participate: Not applicable. •Consent for publication: Not applicable. •Data availability: All data generated are included in the manuscript. •Materials availability: Not applicable. •Code availability: Codes can be provided per request. •Author contribution: 35 – Fu Li: Conceptualization, Methodology, Analysis, Writing-Original Draft, Numerical experiments. – Yingxiang Xu: Conceptualization, Methodology, Analysis, Writing-Review & Editing, Supervision. Appendix A: Proof of Theorem 4 Proof. Using the ansatz q∗=CqLαwith L=C1hand asymptotically expanding both sides of (26) for hsmall, we find 1−4kminbmin a2 min +b2 min CqCα 1hα+··· = 1 −2√2 √εCCq C−α 1hθ 2−α−r2C εC1h1−θ 2+··· .(48) Using the same technique as in the proof of Theorem 1, we obtain 4kminbmin a2 min +b2 min CqCα 1hα=         2√2 √εCCqC−α 1hθ 2−α,if 0 < θ < 4 3, 2√2 √εCCqC−α 1hθ 2−α+q2C εC1h1−θ 2,if θ=4 3, q2C εC1h1−θ 2if 4 3< θ < 2, which leads to the solution        α∗=θ 4, C∗ q= 23 4(εC)−1 4C−θ 4 1ς−1 2 min,if 0 < θ < 4 3, α∗=1 3, C∗ q= (2εC)−1 2C−1 3 1ς−1 min(qC2C2 1+ 23 2(εC)1 2ςmin +CC1),if θ=4 3, α∗= 1 −θ 2, C∗ q= 21 2ε−1 2C1 2C θ 2 1ς−1 min,if 4 3< θ < 2, where ςmin = 4kminbmin/(a2 min +b2 min). Thus, the parameter q∗in (27) is announced. B: Proof of Theorem 5 Proof. The proof can also be divided into four steps. As the first step, we, based on Theorem 4, determine the solution q∗to the equi-oscillation problem in the admissible frequency domain K. We assume again q∗=CqLα. From Lemmas 1and 3, we know that the location of the unique interior maximum of ρC(k, L, q) in the asymptotic sense is ¯ k∗:= ¯ k(q∗) = εq∗+pε2(q∗)2−2εq∗L−L2 ε(q∗)2L=2 Cq (C1h)−(1+α)+O(h−2α), where we have used the assumption L=C1h. Similarly to the proof of Theorem 2, we need to consider the relationship between the maximum frequency kmax =π C2hδand ¯ k∗. 36 If kmax ≥¯ k∗, which occurs when δ > 1 + α, or when δ= 1 + αand C1≥ (2C2/(πCq)) 1 1+α, we should equi-oscillate between kmin and ¯ k∗. Thus, let ˜ k∗=¯ k∗= 2 CqC1h−(1+α)+O(h−2α), i.e., θ= 1 + α, C =2 CqC1+α 1 in Theorem 4, equation (48) now reads 1−ςminCqCα 1hα+··· = 1 −4 pεCq C 1 2−α 2 1h1 2−α 2+··· , which leads to α∗=1 3, C∗ q= 24 3ε−1 3ς−2 3 min.(49) Using these results, one finds that kmax ≥¯ k∗is equivalent to δ > 4 3,or δ=4 3and C1≥2−1 4(C2/π)3 4ε1 4ς 1 2 min =: Cc. Then, when δ < 4 3, or δ=4 3and C1< Cc, which means kmax <¯ k∗, thus we need to equi-oscillate between kmin and kmax. Choosing ˜ k∗=kmax =π C2hδin Theorem 4, one finds (α∗=1 3, C∗ q=˜ Cq,if δ=4 3and C1< Cc, α∗=δ 4, C∗ q= 23 4(επ)−1 4C 1 4 2C−δ 4 1ς−1 2 min,if δ < 4 3,(50) where ˜ Cq= (2επC2)−1 2C−1 3 1ς−1 min(qπ2C2 1+ (2C2)3 2(επ)1 2ςmin +πC1). By combining equations (49) and (50), we obtain the expression for the optimized parameter in (29). In the second step, we show that ρC(k, h, q∗) attains its maximum asymptotically at either k=kmin or ¯ k∗. Using the same proof as Theorem 2, we know that kmin and ¯ kare the only two maximum points of ρC(k, h, q∗) for k∈K. In the third step, we show that q∗in (29) solves the min-max problem (25) in an asymptotic sense. To this end, we only need to show that any parameters asymptotically different from q∗would yield a maximum value of ρCgreater than ρC(kmin, L, q∗). For brevity, we show only the detailed proof for the first case, i.e., kmax >¯ k∗, and the other two cases can be obtained in a similar way. Assuming q=CqLα, we only need to show that if α=α∗=1 3, or Cq=C∗ q, it holds maxkρC(k, L, q)>maxkρC(k, L, q∗) for hsmall enough. When α < α∗, an expansion for hsmall gives ρC(¯ k, L, q)=1−4 pεCq C 1−α 2 1h1−α 2+O(h1−α). Since 1−α 2>1 3for α < 1 3, we conclude that for sufficiently small hit holds ρC(¯ k, L, q)> ρC(¯ k, L, q∗)=1−O(h1 3). When α > α∗, we have for hsmall ρC(kmin, L, q) = 1 −ςminCqCα 1hα+O(hmin{2α,1}), 37 which is definitely greater than ρC(kmin, L, q∗) for sufficiently small h. We now consider the case α=α∗but Cq=C∗ q. By the asymptotic expansions above we find that ρC(kmin, L, q)> ρC(kmin, L, q∗) if Cq> C∗ qand ρC(¯ k, L, q)> ρC(¯ k, L, q∗) if Cq< C∗ q. Finally, to show the convergence factor estimate, we only expand ρC(kmin, L, q∗) with q∗defined in (29) for hsmall. C: Proof of Theorem 6 Proof. The convergence factor simplifies to ρC(k, 0, q) = a(k)2+ (b(k)−qk)2 a(k)2+ (b(k) + qk)2, and the free parameter qcan then be determined by solving the min-max problem min q>0max k∈KρC(k, 0, q). If we make ansatz q∗=c∗ qk−1 4 max, and assume kmax is large enough, then we have ρC(kmin,0, q∗) = 1 −ςminc∗ qk−1 4 max +O(k−1 2 max),(51a) ρC(kmax,0, q∗) = 1 −2√2/(√εc∗ q)k−1 4 max +O(k−1 2 max).(51b) where ςmin = 4kminbmin/(a2 min +b2 min) Matching the coefficients of the term k−1 4 max in (51a) and (51b), i.e., setting ςminc∗ q= 2√2/(√εc∗ q), one finds c∗ q= 23 4ε−1 4ς−1 2 min, which leads to (30). We now need to show that ρC(k, 0, q) attains its maximum asymptotically at either k=kmin or k=kmax. Actually, this analysis is very similar to those for overlapping case, we thus omit the detail. We are in the position to show that q∗given by (30) solves the min-max problem (25) with L= 0 asymptotically. Similar to the previous analysis, assuming that q= cqk−α max, we only need to show that if α=α∗=1 4or cq=c∗ q, it holds maxkρC(k, 0, q)> maxkρC(k, 0, q∗) for sufficiently large kmax. When α < α∗, we have ρC(kmax,0, q)=1−2√2/(√εcq)k−1−2α 2 max +O(k−2α max). Since 1−2α 2>1 4for α < 1 4=α∗, we claim that for sufficiently large kmax it holds ρC(kmax,0, q)> ρC(kmax,0, q∗)=1−O(k−1 4 max). When α > α∗, we consider the convergence factor at kmin, ρC(kmin,0, q)=1−ςmincqk−α max +O(k−2α max), 38 which is bigger than ρC(kmin,0, q∗) if kmax is large enough. We then consider the case α=α∗but cq=c∗ q. By the asymptotic expansions above we find that if cq> c∗ q, we have ρC(kmin,0, q)> ρC(kmin,0, q∗), and if cq> c∗ qwe obtain ρC(kmax,0, q)> ρC(kmax,0, q∗), which concludes the proof of asymptotic optimality. To show the convergence factor estimate, we only insert c∗ qand α∗into (51a). References [1] Adey, R.A., Brebbia, C.A.: Efficient method for solution of viscoelastic problems. J. Eng. Mech. Div. 99(6), 1119–1127 (1973) [2] Andrews, G.: On the existence of solutions to the equation utt =uxxt +σ(ux)x. J. Differential Equations 35(2), 200–231 (1980) [3] Barrera, J., Volkmer, H.: Asymptotic expansion of the L2-norm of a solution of the strongly damped wave equation. J. Differential Equations 267(2), 902–937 (2019) [4] B´ecache, E., Ezziani, A., Joly, P.: A mixed finite element approach for viscoelastic wave propagation. Comput. Geosci. 8(3), 255–299 (2005) [5] Beerends, R.J.: Fourier and Laplace transforms. Cambridge University Press, Cambridge (2003) [6] Bennequin, D., Gander, M.J., Gouarin, L., Halpern, L.: Optimized Schwarz waveform relaxation for advection reaction diffusion equations in two dimensions. Numer. Math. 134(3), 513–567 (2016) [7] Bennequin, D., Gander, M.J., Halpern, L.: A homographic best approximation problem with application to optimized Schwarz waveform relaxation. Math. Comp. 78(265), 185–223 (2009) [8] Blayo, E., Halpern, L., Japhet, C.: Optimized Schwarz waveform relaxation algorithms with nonconforming time discretization for coupling convection-diffusion problems with discontinuous coefficients. In: Domain decomposition methods in science and engineering XVI, pp. 267–274. Springer, Berlin, Heidelberg (2007) [9] Califano, G., Conte, D.: Optimal Schwarz waveform relaxation for fractional diffusion-wave equations. Appl. Numer. Math. 127, 125–141 (2018) [10] Chen, W., Ikehata, R.: Decay properties and asymptotic behaviors for a wave equation with general strong damping. J. Math. Anal. Appl. 519(1), 126765 (2023) [11] Dolean, V., Gander, M.J., Veneros, E.: Asymptotic analysis of optimized Schwarz methods for Maxwell’s equations with discontinuous coefficients. ESAIM Math. Model. Numer. Anal. 52(6), 2457–2477 (2018) [12] Engquist, B., Zhao, H.K.: Absorbing boundary conditions for domain decomposition. Appl. Numer. Math. 27(4), 341–365 (1998) [13] Evans, L.C.: Partial Differential Equations, vol. 19. American Mathematical Society (2022) [14] Gander, M.J.: Optimized Schwarz methods. SIAM J. Numer. Anal. 44(2), 699– 731 (2006) 39 [15] Gander, M.J., Al-Khaleel, M., Ruchli, A.E.: Optimized waveform relaxation methods for longitudinal partitioning of transmission lines. IEEE Transactions on Circuits and Systems I: Regular Papers 56(8), 1732–1743 (2008) [16] Gander, M.J., Halpern, L.: Absorbing boundary conditions for the wave equation and parallel computing. Math. Comp. 74(249), 153–176 (2005) [17] Gander, M.J., Halpern, L.: Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems. SIAM J. Numer. Anal. 45(2), 666–697 (2007) [18] Gander, M.J., Halpern, L., Nataf, F.: Optimal convergence for overlapping and non-overlapping Schwarz waveform relaxation. In: Eleventh international Conference of Domain Decomposition Methods. Greenwich (Great Britain) (1999) [19] Gander, M.J., Halpern, L., Nataf, F.: Optimal Schwarz waveform relaxation for the one dimensional wave equation. SIAM J. Numer. Anal. 41(5), 1643–1681 (2003) [20] Gander, M.J., Magoul`es, F., Nataf, F.: Optimized Schwarz methods without overlap for the Helmholtz equation. SIAM J. Sci. Comput. 24(1), 38–60 (2002) [21] Gander, M.J., Ruehli, A.E.: Optimized waveform relaxation methods for RC type circuits. IEEE Transactions on Circuits and Systems I: Regular Papers 51(4), 755–768 (2004) [22] Gander, M.J., Stuart, A.M.: Space-time continuous analysis of waveform relaxation for the heat equation. SIAM J. Sci. Comput. 19(6), 2014–2031 (1998) [23] Gander, M.J., Thibaut, L., Ausra, P.: Convergence of parareal for a vibrating string with viscoelastic damping. In: Domain Decomposition Methods in Science and Engineering XXVI 8(1), 1–8 (2021) [24] Gander, M.J., Xu, Y.: Optimized Schwarz methods for circular domain decompositions with overlap. SIAM J. Numer. Anal. 52(4), 1981–2004 (2014) [25] Gander, M.J., Xu, Y.: Optimized Schwarz methods with nonoverlapping circular domain decomposition. Math. Comp. 86(304), 637–660 (2017) [26] Gander, M.J., Zhao, H.: Overlapping Schwarz waveform relaxation for parabolic problems in higher dimension. In: Proceedings of Algoritmy, vol. 14, pp. 42–51 (1997) [27] Giladi, E., Keller, H.: Space-time domain decomposition for parabolic problems. Appl. Math. 217, 50 (1997) [28] Halpern, L., Japhet, C.: Discontinuous Galerkin and nonconforming in time optimized Schwarz waveform relaxation for heterogeneous problems. In: Domain Decomposition Methods in Science and Engineering XVII, pp. 211–219. Springer, Berlin, Heidelberg (2008) [29] Halpern, L., Japhet, C., Szeftel, J.: Optimized Schwarz waveform relaxation and discontinuous Galerkin time stepping for heterogeneous problems. SIAM J. Numer. Anal. 50(5), 2588–2611 (2012) [30] Ikehata, R., Takeda, H.: Asymptotic profiles of solutions for structural damped wave equations. J. Dynam. Differential Equations 31, 537–571 (2019) [31] Jiang, Y., Ding, X.: Waveform relaxation methods for fractional differential equations with the Caputo derivatives. J. Comput. Appl. Math. 238, 51–67 (2013) [32] Kalyani, V.K., Pallavika, Chakraborty, S.K.: Finite-difference time-domain 40 method for modelling of seismic wave propagation in viscoelastic media. Appl. Math. Comput. 237, 133–145 (2014) [33] Lelarasmee, E., Ruehli, A.E., Sangiovanni-Vincentelli, A.L.: The waveform relaxation method for time-domain analysis of large scale integrated circuits. IEEE transactions on computer-aided design of integrated circuits and systems 1(3), 131–145 (1982) [34] Li, F., Xu, Y.: A diagonalization-based parallel-in-time algorithm for CrankNicolson’s discretization of the viscoelastic equation. East Asian J. Appl. Math. 14(1), 47–78 (2024) [35] Li, H.: Generalized difference methods for one-dimensional viscoelastic problems. J. Korean Soc. Ind. Appl. Math. 9(2), 55–64 (2005) [36] Lin, Y.: A mixed type boundary problem describing the propagation of disturbances in viscous media I, weak solution for quasi-linear equations. J. Math. Anal. Appl. 135, 644–653 (1988) [37] Lions, P.L.: On the Schwarz alternating method. III: a variant for nonoverlapping subdomains. In: Third international symposium on domain decomposition methods for partial differential equations, vol. 6, pp. 202–223. SIAM, Philadelphia, PA (1990) [38] Luo, Z., Jin, S.: A reduced-order extrapolated Crank-Nicolson collocation spectral method based on proper orthogonal decomposition for the two-dimensional viscoelastic wave equations. Numer. Methods Partial Differential Equations 36(1), 49–65 (2020) [39] Martin, V.: An optimized Schwarz waveform relaxation method for the unsteady convection diffusion equation in two dimensions. Appl. Numer. Math. 52(4), 401–428 (2005) [40] Nikan, O., Avazzadeh, Z.: Coupling of the Crank-Nicolson scheme and localized meshless technique for viscoelastic wave model in fluid flow. J. Comput. Appl. Math. 398, 113695 (2021) [41] Oru¸c, ¨ O.: Two meshless methods based on local radial basis function and barycentric rational interpolation for solving 2D viscoelastic wave equation. Comput. Math. Appl. 79(12), 3272–3288 (2020) [42] Pata, V., Squassina, M.: On the strongly damped wave equation. Comm. Math. Phys. 253(3), 511–533 (2005) [43] Schwarz, H.A.: Ueber einen Grenz¨ubergang durch alternirendes Verfahren. Vierteljahrsschrift der Naturforschenden Gesellschaft in Z¨urich 15, 272–286 (1870) [44] Shibata, Y.: On the rate of decay of solutions to linear viscoelastic equation. Math. Methods Appl. Sci. 23(3), 203–226 (2000) [45] Shukla, K., Chan, J., Maarten, V.: A high order discontinuous Galerkin method for the symmetric form of the anisotropic viscoelastic wave equation. Comput. Math. Appl. 99, 113–132 (2021) [46] Stucky, P., Lord, W.: Finite element modeling of ultrasonic waves in viscoelastic media. In: Review of Progress in Quantitative Nondestructive Evaluation, pp. 113–120. Springer, Boston, MA (1997) [47] Wang, X., Gao, F., Sun, Z.: Weak Galerkin finite element method for viscoelastic 41