Numer. Math. (2002) 93: 279–313 Digital Object Identifier (DOI) 10.1007/s002110100345

Numerische Mathematik

Space-time domain decomposition for parabolic problems Eldar Giladi1 , Herbert B. Keller2 1 2

Incyte Genomics, 3160 Porter Drive, Palo Alto, CA 94304, USA; e-mail: [email protected] Applied Mathematics 217-50, Caltech, Pasadena CA 91125, USA; e-mail: [email protected]

Received April 7, 1999 / Revised version received November 12, 2000 / c Springer-Verlag 2002 Published online April 17, 2002 –


Summary. We analyze a space-time domain decomposition iteration, for a model advection diffusion equation in one and two dimensions. The discretization of this iteration is the block red-black variant of the waveform relaxation method, and our analysis provides new convergence results for this scheme. The asymptotic convergence rate is super-linear, and it is governed by the diffusion of the error across the overlap between subdomains. Hence, it depends on both the size of this overlap and the diffusion coefficient in the equation. However it is independent of the number of subdomains, provided the size of the overlap remains fixed. The convergence rate for the heat equation in a large time window is initially linear and it deteriorates as the number of subdomains increases. The duration of the transient linear regime is proportional to the length of the time window. For advection dominated problems, the convergence rate is initially linear and it improves as the the ratio of advection to diffusion increases. Moreover, it is independent of the size of the time window and of the number of subdomains. Numerical calculations illustrate our analysis. Mathematics Subject Classification (1991): 65M55 1 Introduction The numerical solution of advection diffusion problems usually involves a two stage process. First the spatial operator is dicretized to yield a stiff system of O.D.Es in time which is then integrated. The waveform relaxation method has gained renewed popularity as a practical scheme for integrating Correspondence to: E. Giladi

280

E. Giladi, H.B. Keller

large systems of O.D.Es on parallel computers [1]-[15]. Here, we consider the numerical solution of advection diffusion equations by the waveformrelaxation iteration applied to the semi-discretized equation. Specifically, we consider the block red-black variant of this iteration with overlapping splittings [16,13]. Our goal is to determine the convergence rate of the algorithm. In order to do so we view this scheme as the discretization of a Schwarz overlapping domain decomposition iteration, which we analyze prior to the discretization. In this fashion we obtain much sharper results then those obtained from the general theory of waveform relaxation. Our analysis focuses on the case where the time interval is finite. Jeltsch and Pohl [16] extended the theory of waveform relaxation to the case of overlapping splittings. They showed that the rate of convergence is super-linear, in a finite time interval. However, the rate of convergence rapidly deteriorates for the semi-discretized advection diffusion equations, as the spatial mesh parameter converges to zero. Our analysis of the continuous iteration is independent of discretization parameters. Moreover, we find that the asymptotic convergence rate is governed by the diffusion of the error, across the overlap between subdomains. This yields a much higher rate of convergence then the one implied by the general theory for the waveformrelaxation iteration. We also find that the convergence rate depends on the size of the overlap between subdomains, the size of the time window and the diffusion coefficient in the equation. It is independent of the number of subdomains, provided the size of the overlap between subdomains remains fixed. Gander and Stuart [17] simultaneously and independently analyzed the iteration for the semi-discretized one dimensional heat equation on an infinite time interval, by formulating the waveform relaxation using Schwarz overlapping domain decomposition. They found that on an infinite time interval convergence is linear and that the rate deteriorates as the number of subdomains increases. Our analysis for the heat equation shows that on large time intervals, convergence is initially linear with the the same convergence rate as for the infinite time interval. Moreover, the duration of the transient linear regime is proportional to the length of the time window. In advection dominated problems, we find that convergence is initially linear with a convergence factor that is governed by the ratio of advection to diffusion. Specifically, it is an exponentially decaying function of magnitude of this ratio, and it is independent of the size of the time window. We note that in a practice the equation is discretized and a numerical solution is computed in each subdomain. When the spatial operator is strongly convection dominated, construction of discrete approximation methods is a difficult task. Some recent developments in the steady state case can be found in [18,19].

Space-time domain decomposition for parabolic problems

281

The results obtained for the one dimensional model problem are also valid in two space dimensions in rectangles, albeit in a different norm. In Sect. 2 we present the algorithm and the 1-D model equations. In Sect. 3 we summarize the results of this paper for the 1-D model equations and we present numerical calculations which demonstrate the accuracy of our results. Then in Sect. 4 we present some preliminary error analysis. In Sect. 5 we derive the asymptotic convergence rate of the iteration, for the one dimensional problem. We also analyze the transient linear rate, for advection dominated equations. In Sect. 6 we analyze the transient linear rate for the heat equation, in a large time window. In Sect. 7 we extend the analysis to the two dimensional advection diffusion equation. We end with some conclusions and future directions in Sect. 8. 2 The model problem and the algorithm We consider the one-dimensional advection-diffusion equation in Ω×[0, ∞) (1)

∂u ∂2u ∂u = a2 2 + b − c2 u + f (x, t), t > 0, ∂t ∂x ∂x

(2)

u(x, 0) = h(x),

(3) (4)

u(xl , t) = gl (t),

u(xh , t) = gh (t),

Ω = [xl , xh ],

and it’s two dimensional counterpart, without the zeroth order term, in Ω × [0, ∞) ∂u = a2 ∆u + b[sin(θ), cos(θ)] · ∇u + f (x, y, t), (5) ∂t (6) u(x, y, 0) = h(x, y),

(7) (8)

u(xl , y, t) = gxl (y, t),

u(xh , y, t) = gxh (y, t),

u(x, yl , t) = gyl (x, t),

u(x, yh , t) = gyh (x, t),

Ω = [xl , xh ] × [yl , yh ].

In equations (1) and (5) the coefficients a, b and c are assumed constants. We solve the problems (1)–(3) and (5)–(7) by a domain decomposition (DD) iteration in which Ω × [0, ∞) is partitioned into space-time subdomains. Here Ω is defined in (4) and (8) for the one and two dimensional case respectively. Time is partitioned into windows of length T and we denote the resulting space-time windows by (9)

Wn = Ω × [nT, (n + 1)T ],

n = 0, 1, · · · .

282

E. Giladi, H.B. Keller

Fig. 1. Decomposition of the domain [x0 , x3 ] × [0, 2T ] into two time windows and six space-time subdomains. The light grey areas are the overlap between the black sub-domain and it’s red neighbors. The overlap of the bottom red sub-domain with it’s black neighbor is shaded in dark grey

Each space-time window is further decomposed into subdomains by partitioning space. In one dimension Ω = [xl , xh ] is decomposed into subintervals of length S, and each subdomain consists of a subinterval and an overlap of size ∆ over adjacent subdomains. In two dimensions Ω = [xl , xh ] × [yl , hh ] is partitioned into strips of width S, and each subdomain consists of a strip and an overlap of width ∆ over adjacent strips. Other possible decompositions are described in Sect. 7. Adjacent space-time subdomains overlap in space, and the width of the overlap is denoted by ∆. Figure 1 depicts the decomposition of [xl , xh ] × [0, 2T ] into two time windows which in turn are decomposed into three space-time subdomains. In the general case we decompose [xl , xh ] into N segments using the points x0 = xl , x1 , · · · , xN = xh where xi − xi−1 = S. The problem is solved successively in W0 , W1 , · · · to accuracy , by a sub-domain iteration with red black ordering. For example, Fig. 1 indicates the decomposition of W0 = [x0 , x3 ] × [0, T ] into three space-time subdomains. The red subdomains are [x0 , x1 +∆]×[0, T ] and [x2 −∆, x3 ]×[0, T ], and the black subdomain is [x1 − ∆, x2 + ∆] × [0, T ]. We denote the red and black iterates by urk (x, t) and ubk (x, t) and we assume that ub0 (x, t) is given. At iteration k ≥ 1, urk (x, t) is computed by solving (1) in [x0 , x1 +∆]×[0, T ] and [x2 − ∆, x3 ] × [0, T ], subject to the conditions urk (x, 0) = h(x),

(10) (11)

urk (x0 , t) = gl (t),

urk (x1 + ∆, t) = ubk−1 (x1 + ∆, t),

(12)

urk (x3 , t) = gh (t),

urk (x2 − ∆, t) = ubk−1 (x2 − ∆, t).

Space-time domain decomposition for parabolic problems

283

Then ubk (x, t) is computed by solving (1) in [x1 − ∆, x2 + ∆] × [0, T ], subject to ubk (x, 0) = h(x),

(13)

(14)

ubk (x2 + ∆, t) = urk (x2 + ∆, t), ubk (x1 − ∆, t) = urk (x1 − ∆, t).

We denote the k’th iterate by uk (x, t),
r uk (x, t) x ∈ [x0 , x1 ] ∪ [x2 , x3 ] (15) uk (x, t) = ubk (x, t) x ∈ [x1 , x2 ], for k ≥ 1. This definition generalizes in a straightforward way to the general case with N subdomains. We note that uk (x, t) is discontinuous in general at x = x1 and x = x2 in this example, and on the lines x = xk in the general case. The initial value in Wn , u(x, nT ), is obtained from Wn−1 for n ≥ 1 and from the initial condition of the problem, for n = 0. Successive time windows do not overlap in time. The error at iteration k is defined as ek = u − uk ,

(16) and we introduce the notation

|ek |t = max |ek (x, τ )|,

(17) (18)

0≤τ ≤t

ek t =

max

x∈Ω, 0≤τ ≤t

|ek (x, τ )|.

3 Summary of results We begin by summarizing the results of this paper for the one dimensional equation, and by presenting numerical calculations which validate these results. Analogous results are derived for the two dimensional heat equation in a rectangle in Sect. 7. 3.1 Results The first result shows that the asymptotic convergence rate of the iteration in a finite time window is super-linear. Theorem 1 Consider the space-time DD iteration for the one dimensional problem (1)–(3), and suppose S > 2∆. Then for any 0 ≤ t ≤ T   |b|(xN −x0 ) [2(k − 1)∆]2

ek t 2 (1 + o(1)), ≤ e 2a exp − a2 t

e0 t (19) k → ∞,

284

E. Giladi, H.B. Keller

The analysis of the semi-discrete iteration by the theory of waveform relaxation and multi-splitting [16], yields the following bound for the decay of the error in the iteration  

ek t (CT )k ≤O

e0 t k!    1 = O exp −k log k + k(log CT + 1) + k , (20) k → ∞. 2 Here, the constant C grows without bound as the mesh parameter decreases to 0, in view of the stiff nature of the system of ODEs that results from the discretization of a parabolic PDE by the method of lines. Moreover, the dependence of the rate of convergence on ∆ and the diffusion coefficient a2 is not apparent in (20). The quadratic exponent in expression (19), which results from the Gaussian nature of diffusion, predicts a substantially faster rate of convergence. We also note that the bound (19) for the continuous iteration is independent of a mesh parameter. The algorithm is particularly well suited to problems with small diffusion coefficients. Indeed, when the ratio ∆2 a2 T is large the convergence rate is high. Hence, a small diffusion coefficient allows the use of a small overlap and a large time window, while maintaining a high rate of convergence. Moreover, it allows the use of a large number of subdomains, without adversely affecting the convergence rate, in view of the inequality S > 2∆. Sharper estimates for the heat equations are presented in [20]. We next consider the case where the ratio of advection to diffusion, |b|/a2 , is large. We find in this case a transient linear regime with a rate of convergence which depends on this ratio. Theorem 2 Consider the space-time DD iteration for the one dimensional problem (1)–(3), and define γ2 =

b2 + c2 . 4a2

For all t ≥ 0 and k ≥ a/(4∆γπ) + 1 (21)

     |b|(xN −x0 ) Sγ 2k−2 2∆γ

ek t 2 2a + exp − exp − ≤ B(k)e .

e0 t a a

Space-time domain decomposition for parabolic problems

Here, B(k) is defined by B(k) = 1 +

1 4



285

a . 2π(k − 1)∆γ

The average convergence factor at iteration k, ρ¯k , is defined by  

ek t 1/k ρ¯k = (22) ,

e0 t and the asymptotic convergence factor, ρ¯, is then ρ¯ = limk→∞ ρ¯k .

(23)

Theorem 2 yields a bound for the asymptotic convergence factor, in the linear regime      Sγ 2 2∆γ + exp − (24) . ρ¯ ≤ exp − a a The improvement of this factor as the ratio a|b|2 increases, follows from the inequality |b| γ ≥ (25) . a 2a2 We also note that this factor is independent of the window size T . The next results characterize the convergence of the iteration for the heat equation, in large time windows. In this case the error initially decays at a linear rate before the transition to the super-linear regime. We derive a lower bound for the error that decay at a linear rate. The lower bound is valid for a finite number of iterations. Here we assume that the initial error e0 (x, t) ≥ 0 and that there exists a small positive number ν and a time Te such that (26)

e0 (x, t) ≥ (1 − ν) e0 t ,

for all x ∈ [x0 + S/2, xN − S/2] and t ≥ Te . Theorem 3 Consider the space-time DD iteration for the one dimensional heat equation, (1) with b = c = 0, in [x0 , xN ] × [0, T ] with N = 2, 3 subdomains. Suppose the initial error satisfies assumption (26). Then for any η > 0 there exists a time Tη such that if k ≤ (T − Te )/(2Tη ) (27)

(1 − η)k ρk−1 (1 − ν)C ≤

where 1 , (28) C = 1 + ∆/S

 ρ=

ek t ≤ ρk−1 C

e0 t

1 − ∆/S 1 + ∆/S

3−L

,

L = N − 1.

286

E. Giladi, H.B. Keller

Theorem 4 Consider the space-time DD iteration for the one dimensional heat equation with N ≥ 4 subdomains in [x0 , xN ] × [0, T ]. Suppose that the initial error satisfies assumption (26). Then for any η > 0 there exists a time Tη such that if k ≤ (T − Te )/(2Tη ) (29)

ek t ≥

e0 t



1 + 2∆/S 1− L

k

(1 − η)k (1 − ν),

L = N −1 for N odd. For N even inequality (29) holds with L−1 substituted for L. We see that convergence factor rapidly deteriorates as the number of subdomains increases. In a practical application, it is desirable to choose T sufficiently small so as to avoid the transient linear regime. Indeed, the super-linear convergence rate does not degrade with the number of subdomains, provided ∆ is fixed. Moreover, the super-linear rate improves as diffusion decreases while the transient linear regime is independent of the diffusion coefficient. 3.2 Numerical calculations We begin by comparing the error estimates of Theorem 1 for the one dimensional problem, with the decay of the error in an actual computation. We solve the problem in W0 = [0, 1] × [0, 0.25], with 3 subdomains of width S = 1/3. The diffusion coefficient is a2 = 1/2 and b = 0. In all our calculations the initial and boundary conditions are (30)

u(x, 0) = sin(πx/(xh − xl )),

(31)

u(xl , t) = u(xh , t) = 0.

In our computations we discretize the equation using the Crank-Nicholson method with mesh parameters ∆x = 1/210 and ∆t = 1/50. We first compute the solution in the whole domain and denote it by uh . Then we solve the problem using the space-time iteration. We denote the k’th iterate by uhk and define the error ehk by (32)

ehk = uh − uhk .

The initial iterate uh0 is chosen to be random in the interior of the domain and satisfies the conditions (30) and (31). Figure 2 indicates (33)

log10

ehk T ,

eh0 T

Space-time domain decomposition for parabolic problems

287

2

0

Log10(|max(error)|)

−2

−4

−6

−8

−10

−12

−14 0

10

20

30 Iteration number

40

50

60

eh  The log10 ekh T 0 T

Fig. 2. is indicated in solid and + for ∆ = 3S/10, ∆ = S/10, ∆ = S/20, from left to right respectively. The dashed line indicates the estimate of the bound 0

−2

−4

−6

−8

−10

−12

−14 0

20

40

60

80

100

120

eh 

Fig. 3. The log10 ekh T for T = 1, T = 2, T = 3 and T = 4, is indicated in + from left 0 T to right, respectively. The dashed line indicates the transient linear bound. The number of subdomains is 3

as a function of the iteration k, for three different values of the overlap ∆. The dashed lines indicate our estimate for the error. Our bound is in very good agreement with the computation and it’s asymptotic nature is apparent. We now compare the observed transient linear regime for the heat equation, in a large time interval, with the predictions of Theorems 3 and 4. In this computation ∆ = S/7 and T varies, T = 1, 2, 3, 4. Figure 3 indicates the error (33) as a function of the iteration number k when the number of

288

E. Giladi, H.B. Keller 2

0

−2

−4

−6

−8

−10

−12 0

20

40

60

80

100

120

eh  log10 ekh T 0 T

for T = 1, T = 2, T = 3 and T = 4, is indicated in + from left Fig. 4. The to right, respectively. The number of subdomains is 17

subdomains is three. The transient linear upper bound of Theorem 3 is also indicated in this figure. It’s validity for a finite number of iterations is apparent. Moreover, the duration of the linear regime is clearly proportional to T , as predicted by the analysis. In a second set of calculations we solve the problem in W0 = [0, 17/3]×T , with 17 subdomains of width S = 1/3. The results are presented in Fig. 4. A comparison with Fig. 3 shows the deterioration of the linear bound as the number of subdomains increases, as predicted by Theorem 4. In contrast, the asymptotic convergence rate, depends only on the size of the overlap as indicated in Theorem 1. We now compare the observed transient linear rate in advection dominated equations, with the analysis of Theorem 2. Here, the diffusion coefficient is a2 = 1/2 and the time window T = 4, is relatively large. The number of subdomains is three. The advection coefficient varies, b = 30, 20, 10. In Fig. 5 we observe rapid convergence, in contrast to the case b = 0, indicated in Fig. 3. The improvement in the convergence rate, as the ratio b/a2 increases, is apparent. In our computations, we measure the average convergence factor, defined in (22). Table 1 presents the measured convergence factor versus the bound for the asymptotic convergence factor of (24). There is good agreement between the two. 4 Error analysis for the 1-D model problem 4.1 The error equations We now derive equations for the error in the DD iteration for problem (1)– (3). We first assume that W0 = [x0 , x2 ] × [0, T ] is decomposed into two

Space-time domain decomposition for parabolic problems

289

2 0 −2 −4 −6 −8 −10 −12 −14 −16 0

Fig. 5. The

2

eh  log10 ehk T 0 T

4

6

8

10

12

14

16

18

20

for b = 30, b = 20 and b = 10, is indicated in solid and + from

left to right, respectively. The number of subdomains is three. Here a2 = 1/2 and T = 4

Fig. 6. The decomposition of W0 = [x0 , x2 ] × [0, T ] into two subdomains

space-time subdomains, as depicted in Fig. 6. The errors in the red and black iterates are defined by (34)

ebk = u − ubk ,

with uk defined by

erk = u − urk ,


uk (x, t) =

ek = u − uk ,

urk (x, t) x ∈ [x0 , x1 ] ubk (x, t) x ∈ [x1 , x2 ]

.

The domain of definition of ebk is [x1 − ∆, x2 ] × [0, T ], that of erk is [x0 , x1 + ∆], while ek is defined in the whole domain. In order to determine erk we subtract equations (1), (10) and (11) for urk from equations (1)–(3) for u to obtain (35)

∂er ∂ 2 er ∂erk = a2 2k + b k − c2 erk , ∂t ∂x ∂x

290

E. Giladi, H.B. Keller

Table 1. The bound for the asymptotic convergence factor versus the measured average convergence factor b

Bound for ρ¯

Measured ρ¯k

10 20 30

0.1777 0.0225 0.0033

0.1766 0.0288 0.0046

(36)

erk (x, 0) = 0,

(37)

erk (x0 , t) = 0, erk (x1 + ∆, t) = ebk−1 (x1 + ∆, t).

In a similar fashion we find that the error in [x1 − ∆, x2 ] × [0, T ] satisfies ∂eb ∂ 2 eb ∂ebk = a2 2k + b k − c2 ebk , ∂t ∂x ∂x

(38) and the conditions (39) (40)

ebk (x, 0) = 0,

ebk (x1 − ∆, t) = erk (x1 − ∆, t), ebk (x2 , t) = 0.

The recursive system of (35)–(40) determines the evolution of the error in iteration k as a function of the error at iteration k − 1. When the number of subdomains N ≥ 3 we derive a recursive system of N equations of the type (35) and (38), coupled through their boundary conditions, using the same method as for the two subdomain case. In order to obtain an explicit expression for the error from the recursive system (35)–(40) we first derive the solution to the generic equation (41)

∂e ∂2e ∂e = a2 2 + b − c2 e, ∂t ∂x ∂x

with e(x, 0) = 0 and Dirichlet conditions at x = xl and x = xh , in Appendix A. We find e(x, t) = eα(x−xl ) K(t, xh − x, a, b, c) ∗ e(xl , t) (42)

+eα(x−xh ) K(t, x − xl , a, b, c) ∗ e(xh , t),

where (43) K(t, z, a, b, c) =

∞ 

{F [t, λn (z), a, γ] − F [t, µn (z), a, γ]} ,

n=0

(44) (45)

b2 b + c2 , α = − 2 2 4a 2a λn (z) = (2n + 1)L − z and µn (z) = (2n + 1)L + z, γ2 =

Space-time domain decomposition for parabolic problems

291

  η η2 2 F (t, η, a, γ) = √ 3/2 exp − 2 − γ t , 4a t 2 πat

(46) and (47)

L = xh − xl .

4.2 Preliminary analysis We begin by deriving simple bounds for expression (42) which we use in the convergence analysis of the iteration, in the next sections. Lemma 1 Let F be defined in (46). Suppose that e(t) ≥ 0 is a nondecreasing function and that the scalars λ, µ satisfy 0 < λ < µ. Then for any γ  T (48) {F (τ, λ, a, γ) − F (τ, µ, a, γ)} e(T − τ )dτ ≥ 0. 0

Proof. The ratio λ F (t, λ, a, γ) = exp F (t, µ, a, γ) µ

(49)



µ2 − λ2 4a2 t

 ,

is a monotonically decreasing function of t that tends to infinity as t → 0 and to λ/µ as t → ∞. Hence, there exists a unique scalar t∗ such that (50)

F (t, λ, a, γ) ≥ F (t, µ, a, γ) for 0 < t ≤ t∗ ,

and (51)

F (t, λ, a, γ) ≤ F (t, µ, a, γ) for t∗ ≤ t.

We now rewrite the integral (48) as the sum 

min(t∗ ,T )

0

 (52)

2

2

[F (τ, λ, a, γ) − F (τ, µ, a, γ)]eγ τ e(T − τ )e−γ τ dτ +

T

min(t∗ ,T )

2

2

[F (τ, λ, a, γ) − F (τ, µ, a, γ)] eγ τ e(T − τ )e−γ τ dτ,

and we note that e(T − t∗ )e−γ (53)

2 t∗

= min e(T − t)e−γ

2t

t∈[0,t∗ ]

2

= max e(T − t)e−γ t , t∈[t∗ ,T ]

292

E. Giladi, H.B. Keller

in view of the monotonicity of e(t). Combining equations (50)–(53) yields the following lower bound for the left side of (48)  T {F (τ, λ, a, γ) − F (τ, µ, a, γ)} e(T − τ )dτ 0   T  λ λ2 ∗ −γ 2 t∗ ≥ e(T − t )e √ 3/2 exp − 2 4a t 0 2 πat    T µ2 µ − (54) dt. √ 3/2 exp − 2 4a t 0 2 πat The change of variables τ = 2aλ√t and τ = 2aµ√t in the first and second integral of expression (54), respectively, yields the following bound and proves the Lemma  T {F (τ, λ, a, γ) − F (τ, µ, a, γ)} e(T − τ )dτ 0

(55)

2 2 ∗ ≥ √ e(T − t∗ )e−γ t π



µ √ 2a T λ √ 2a T

2

e−τ dτ.  

We now introduce the following truncated sums for the kernel (43) KN (t, z, a, b, c) = F [t, λ0 (z), a, γ] +

N −1 

F [t, λn+1 (z), a, γ]

n=0

(56)

(57)

−F [t, µn (z), a, γ] ,

KN (t, z, a, b, c) =

N 

(F [t, λn (z), a, γ] − F [t, µn (z), a, γ]) ,

n=0

with γ defined in (44), and where the sum in (56) is absent when N = 0. We now prove Corollary 1 Suppose that e(t) ≥ 0 is a non-decreasing function and that 0
≤ KN (t, z, a, b, c) ∗ e(t).

Proof. We obtain both inequalities in (58) by determining the sign of the remainders K(t, z, a, b, c)∗e(t)−KN (t, z, a, b, c)∗e(t) and K(t, z, a, b, c)∗ e(t) − KN (t, z, a, b, c) ∗ e(t), using Lemma 1 and the fact that (59)

λn (z) < µn (z) < λn+1 (z).

Space-time domain decomposition for parabolic problems

Here, λn and µn are defined in (45).

293

 

The special case N = 0 of Corollary 1 is Corollary 2 Suppose that e(t) ≥ 0 is a non-decreasing function and that 0
K(t, z, a, b, c) ∗ e(t) ≤ F [t, λ0 (z), a, γ] ∗ e(t).

We note that in this inequality, the expression on the left corresponds to the solution of (41) in a bounded domain, with a homogeneous condition on one of the boundaries. The formula on the right corresponds to the solution in the quarter plane, with the requirement that |e| → 0 as |x| → ∞. It is obtained by neglecting all corrections introduced by the images in the sum (43), for the derivative of the quarter plane Green’s function. In view of the non-negativity of F and e in the right of inequality (60) we find that (61)

F [t, λ0 (z), a, γ] ∗ e(t) ≤ max e(τ ) I(t, λ0 (z), a, γ), τ ≤t



where I(t, η, a, γ) =

(62)

t

0

F (τ, η, a, γ)dτ.

In order to estimate I, we introduce the change of variables η √ σ = (63) 2a τ in (62) and we find using [21]

(64)



 √ η √ +γ t e erfc 2a t   √ η − γη a √ +e erfc . −γ t 2a t

1 I(t, η, a, γ) = 2



γη a

Here erfc(z) is the complementary error function  ∞ 2 2 √ (65) e−t dt, erfc(z) = π z which admits the expansion (see: [21]) 2

(66) and the bound

e−z erfc(z) = √ (1 + o(1)), z → ∞, πz 2

e−z erfc(z) ≤ √ , πz We now derive bounds for I. (67)

z > 0.

294

E. Giladi, H.B. Keller

Lemma 2 Let I be defined in (62) and suppose that η, a > 0. 1. If γ = 0 then for any t > 0 (68)

I(t, η, a, γ) ≤

2. If γ > 0 and η ≥ (69)

  η2 t 2a exp − 2 π η 4a t

a γπ

then for any t > 0      η2 a 1 exp − 2 I(t, η, a, γ) ≤ 1 + 2 2πηγ 4a t

3. If γ > 0 then for all t ≥ 0 (70)



I(t, η, a, γ) ≤

 1+

1 2



a 2πηγ



ηγ

exp − a

Proof. The bound (67) yields     √ γη η 1 η2 2 a

√ +γ t ≤ √ e erfc √ exp − 4a2 t − γ t . 2a t π 2aη√t + γ t (71) The exponential and the pre-exponential factors in this expression are maximized at η t = (72) , 2aγ and this yields the bound    ηγ

√ γη a η √ +γ t ≤ (73) exp − . e a erfc 2πηγ a 2a t We use this inequality and the fact that erfc(z) ≤ 2 in expression (64) to obtain (70). The substitution (72) for t in the pre-exponential factor of expression (71) yields      √ γη η2 η a a √ +γ t ≤ e erfc exp − 2 . (74) 2πηγ 4a t 2a t η The bound (67) implies that for t ≤ 4aγ      √ η2 a η − γη √ −γ t ≤ 2 (75) exp − 2 , e a erfc γηπ 4a t 2a t η while inequality erfc(z) ≤ 2 implies that for t ≥ 4aγ     √ η η2 − γη a √ (76) e erfc − γ t ≤ 2 exp − 2 . 4a t 2a t

Space-time domain decomposition for parabolic problems

295

We now use inequalities (74)–(76) in (64) to obtain (69). Inequality (68) follows from identity (64) and the bound (67) for erfc(z).   Lemma 3 Suppose that e(x, t) is the solution to (41) in [xl , xh ] × [0, T ] with e(x, 0) = 0 and Dirichlet conditions at x = xl and x = xh . Suppose that e¯(x, t) is also a solution to (41) in [xl , xh ] × [0, T ] with e¯(x, 0) = 0 and that e¯(xl , t) ≥ |e(xl , t)| , e¯(xh , t) ≥ |e(xh , t)| . (77) Then for all (x, t) ∈ [xl , xh ] × [0, T ] e¯(x, t) ≥ |e(x, t)|.

(78)

Proof. Both v = e¯−e and w = e¯+e satisfy (41) with non-negative boundary conditions. Hence v and w are non-negative, as shown in Appendix B for (41). This proves inequality (78).   Lemma 4 Suppose e(t) ≥ 0 is a non-decreasing function and that h(t) ≥ 0. The convolution (79) H(t) = h(t) ∗ e(t), is a non-negative and non-decreasing function of t. Proof. Suppose t2 ≥ t1 ≥ 0  H(t2 ) − H(t1 ) =

0

t1



+

(80)

h(τ ) [e(t2 − τ ) − e(t1 − τ )] dτ t2

t1

h(τ )e(t2 − τ )dτ.

In view of the monotonicity and non-negativity of e both integrands in (80) are non-negative. Hence, H(t2 ) ≥ H(t1 ) and this proves the lemma.   Lemma 5 Suppose θ > 0 and let  ak (j) =

(81)

2k − 2 j



[(4(k − 1)∆ + jθ]2 exp − 4a2 t

 .

Then (82)

2k−2  j=0



[(2(k − 1)∆]2 ak (j) = exp − a2 t

 (1 + o(1)), k → ∞.

296

E. Giladi, H.B. Keller

Proof. We derive the asymptotic estimate (82) by essentially using Laplace’s method for sums [22]. First we note that for all j = 0, . . . , 2k − 3   θ2(k − 1)∆ ak (j + 1) (83) ≤ 2(k − 1) exp − . ak (j) a2 t It follows that for sufficiently large k, say k > k0 ak (j + 1) < 1, ak (j)

(84)

j = 0, . . . , 2k − 3,

and ak (j) is a monotonically decreasing function of j. Hence for all k > k0 (85)

2k−2  j=0

   ak (1) [(2(k − 1)∆]2 1 + 2(k − 1) , ak (j) ≤ exp − a2 t ak (0)

and this proves the Lemma, in view of inequality (83).

 

5 Convergence analysis for the 1-D model problem We now consider the DD iteration of Sect. 2 for the one dimensional model problem in W0 = [x0 , xN ] × [0, T ]. We first derive a general bound, which we then use to obtain the asymptotic convergence rate of the iteration, and the transient linear regime for advection dominated equations. 5.1 General analysis We now prove Proposition 1 Consider the space-time DD iteration for the one dimensional problem (1)–(3), and suppose S > 2∆. Then for any 0 ≤ t ≤ T 2k−2   2k − 2 

ek t |α|(xN −x0 ) (86) I[t, 4(k − 1)∆ + jθ, a, γ], ≤ e j

e0 t j=0

(87)

θ = S − 2∆,

where α and γ are defined in (44). Proof. We derive bounds for |erk (x, t)| and |ebk (x, t)| which we denote by e¯rk (x, t) and e¯bk (x, t), respectively. We use an inductive construction which ensures that the bounds are non-decreasing functions of t. Initially we define (88)

e¯b0 (x, t) = Aeαx e0 t ,

Space-time domain decomposition for parabolic problems

 (89)

A =

297

e−αxN α < 0 , e−αx0 α ≥ 0

and we note that e¯b0 is a non-decreasing function of t. For k ≥ 1 we derive the bound for two subdomains the red and black subdomains [x0 , x1 ]×[0, T ] and [x1 , x2 ] × [0, T ], as indicated in Fig. 6. The general case can be derived in a similar fashion and is presented in detail in [20]. The error in the red sub-domain at iteration k ≥ 1, erk , is determined by b ek−1 , through equations (35)–(37). In order to define e¯rk we first substitute e¯bk−1 for ebk−1 in boundary condition (37), and solve the resulting problem using formula (42) to obtain (90) |erk (x, t)| ≤ eα(x−x1 −∆) K(t, x − x0 , a, b, c) ∗ e¯bk−1 (x1 + ∆, t). This inequality follows from Lemma 3 and the inequality |ebk−1 | ≤ e¯bk−1 . We now define e¯rk as the bound for K ∗ e¯bk−1 obtained from Corollary 2 (91) e¯rk (x, t) = eα(x−x1 −∆) F [t, λ0 (x − x0 ), a, γ] ∗ e¯bk−1 (x1 + ∆, t). The monotonicity of e¯rk follows from the monotonicity of e¯bk−1 and Lemma 4. Similar arguments yield (92) |ebk (x, t)| ≤ eα(x−x1 +∆) K(t, x2 − x, a, b, c) ∗ e¯rk (x1 − ∆, t), (93) e¯bk (x, t) = eα(x−x1 +∆) F [t, λ0 (x2 − x), a, γ] ∗ e¯rk (x1 − ∆, t). In equations (91) and (93) λ0 (z) = S + ∆ − z, in view of definition (45) for λ0 . The maximum principle for (35) (see: Appendix B) and Lemma 3 imply (94)

|erk (x, t)| ≤ max |ebk−1 (x1 + ∆, τ )|,

(95)

|ebk (x, t)|

0≤τ ≤t

≤ max |erk (x1 − ∆, τ )|. 0≤τ ≤t

It follows from these inequalities and the monotonicity of e¯bk−1 that (96)

ek t ≤ e¯bk−1 (x1 + ∆, t).

We evaluate this bound by combining equations (91) and (93) repeatedly, using the identity (97)

F (t, η1 , a, γ) ∗ F (t, η2 , a, γ) = F (t, η1 + η2 , a, γ),

which follows from identity (178) in appendix A and the convolution Theorem for the Laplace transform, to obtain (98) e¯bk (x1 + ∆, t) = Aeα(x1 +∆) F [t, 4(k − 1)∆, a, γ] ∗ e0 t .

298

E. Giladi, H.B. Keller

Equations (96), (98) and the monotonicity of e0 t yield  t α(x1 +∆) (99) ek t ≤ Ae

e0 t F [τ, 4(k − 1)∆, a, γ] dτ, 0

and it follows from definitions (62) for I and (89) for A that

ek t (100) ≤ e|α|(x2 −x0 ) I(t, 4(k − 1)∆, a, γ).

e0 t

 

5.2 Asymptotic convergence rate We now show that the rate of convergence of the space-time DD iteration is super-linear. Proof of Theorem 1: We use the bounds (68) and (69) in inequality (86) of Proposition 1 to obtain  2k−2  |b|(xN −x0 ) 

ek t 2k − 2 2 2a ≤ B(k)e j

e0 t j=0   [4(k − 1)∆ + jθ]2 (101) . × exp − 4a2 t This inequality is valid for all k when γ = 0 with √ a t B(k) = √ (102) , π2(k − 1)∆ and it is valid for k ≥ a/(4∆γπ) + 1 when γ > 0 with  a 1 B(k) = 1 + (103) . 4 2π(k − 1)∆γ Inequality (19) follows from the bound (101), in view of Lemma 5.

 

5.3 Advection dominated equations We now analyze the iteration when the ratio of advection to diffusion, |b|/a2 , is large. Proof of Theorem 2: Proposition 1 and inequality (70) of Lemma 2 yield the bound  2k−2  |b|(xN −x0 ) 

ek t 2k − 2 2 ≤ B(k)e 2a j

e0 t j=0

(104)

× exp (−[4(k − 1)∆ + jθ]γ/a) ,

which is valid for k ≥ a/(4∆γπ)+1. This bound is equivalent to inequality (21), in view of the definition of θ in (87).  

Space-time domain decomposition for parabolic problems

299

6 The heat equation in a large time window We now study the convergence of the DD iteration for the heat equation, in large time windows. The error initially decays at a linear rate before it transitions to the super-linear regime. 6.1 The large time limit We begin by proving a few preliminary Lemmas. Lemma 6 The function K(t, z, a, b, c) in (43) is continuous in t and nonnegative for t ≥ 0 and 0 ≤ x ≤ L. Proof. Continuity follows from the continuity of F in (46) and the uniform convergence of the sum (43), which we show by noting that (105)

 3/2 a2 1 6 √ 2 2. max F (t, λn (z), a, γ) ≤ 0≤t e 8 πL n

In order to show that K(t, z, a, b, c) ≥ 0 we consider expression (42), which is the solution to (41) with Dirichlet data. This solution is non-negative when the data is non-negative, as shown in Appendix B. Hence K ≥ 0, in view of it’s continuity.   We now introduce the function φ(z, L, a, T, t) =

(106)

∞  n=0

χ[ nL+β (n+1)L−β (t), √ , √ ] a T

a T

where

L−z , 2 and χ is the indicator function of an interval. In the following Lemmas we essentially show that for a continuous integrable function f   ∞ z ∞ f φdt = f dt. lim T →∞ 0 L 0 β=

(107)

Lemma 7 Let φ(z, L, a, T, t) be defined in (106) and suppose that 0 ≤ tj−1 < tj ≤ t and that 0 < z < L then  (108)

lim

T →∞ 0

t

χ[tj−1 ,tj ] φ(z, L, a, T, τ )dτ =

z L

 0

t

χ[tj−1 ,tj ] .

300

E. Giladi, H.B. Keller

Proof. For a given T , we define the integers n1 and n2 such that (109)

(n1 − 1)L n1 L √ ≤ tj−1 ≤ √ , a T a T

n2 L (n + 1)L √ ≤ tj ≤ 2 √ . a T a T

Hence,

(110)

 t   t    χ[t ,t ] φ(z, L, a, T, τ ) − z  χ [t ,t ] j−1 j j−1 j   L 0 0  n2 L   a√T  2L z   ≤ √ +  φ(z, L, a, T, τ )dτ − (tj − tj−1 ) .  n√1 L  L a T a T

Now, we note that  (111)

n2 L √ a T n1 L √ a T

φ(z, L, a, T, τ )dτ =

z (n2 − n1 )L √ , L a T

and we use this expression and (109) in inequality (110) to obtain  t   t   z  χ[t ,t ] φ(z, L, a, T, τ ) − χ[tj−1 ,tj ]  j−1 j  L 0 0 2 (112) ≤ √ (L + z). a T Taking the limit T → ∞ in (112) proves the Lemma.

 

Lemma 8 Let K(t, z, a, b, c) be defined in (43) and suppose that 0 < z < L. Let  T g(T ) = (113) K(τ, z, a, 0, 0)dτ, 0

then g(T ) is non-decreasing and (114)

lim g(T ) =

T →∞

z . L

Proof. The monotonicity of g follows from the non-negativity of the kernel K(t, z, a, 0, 0) which was noted in Lemma 6. In order to prove the limit (114), we use the definition of K in (43) and evaluate the integral (113) using equations (64) and (65) to obtain  ∞ 2 2 g(T ) = √ (115) e−t φ(z, L, a, T, t)dt, π 0

Space-time domain decomposition for parabolic problems

301

where φ is defined in (106). Now, we use the representation (115) of g(T ) and the fact that erfc(0) = 1 to note that for any  > 0 there exists a t such that     t  z   2  z t −τ 2   −τ 2 + √  e dτ − e φ(z, L, a, T, τ )dτ  .  − g(T ) ≤ L 3 π L 0 0 (116) Then, we define (117)

dt = t /n,

(118)

tj = jdt,

fn (τ ) =

n  j=1

j = 0, 1, . . . , n, 2

χ[tj−1 ,tj ] e−tj−1 ,

which we use in inequality (116) and show that for sufficiently large n  2 z 2    − g(T ) ≤  + √ L 3 π   t   t z   × (119) fn (τ )dτ − φ(z, L, a, T, τ )fn (τ )dτ  , L 0 0 in view of the integrability of fn and φfn . The right of inequality (119) is bounded by , for sufficiently large T , in view of Lemma 7. This proves the Lemma, since  was arbitrary.   Corollary 3 Let K be defined in (43) and 0 ≤ z ≤ L then 1. For any t > 0 (120)

z |e|t L 2. If in addition e ≥ 0 and e(t) ≥ e(Te ) for all t ≥ Te . Then for any η > 0 there exists a number Tη such that for all t ≥ Tη + Te z (121) e(t) ∗ K(t, z, a, 0, 0) ≥ (1 − η)e(Te ) L |e(t) ∗ K(t, z, a, 0, 0)| ≤

Proof. Inequality (120) follows from the mean value theorem for integrals and Lemma 8. For a given η there exists Tη such that  Tη z (122) K(τ, z, a, 0, 0)dτ ≥ (1 − η), L 0 in view of Lemma 8. It follows that for t ≥ Tη + Te  Tη (123) e(t) ∗ K(t, z, a, 0, 0) ≥ e(t − τ )K(τ, z, a, 0, 0)dτ, 0

and that e(t − τ ) ≥ e(Te ) for all τ in the range of integration. Combining this fact with inequalities (123) and (122) yields the bound (121).  

302

E. Giladi, H.B. Keller

Corollary 4 Suppose e(x, t) satisfies (41), with b = c = 0 in [xl , xh ] × [0, T ]. Let gl (t) = e(xl , t), gh (t) = e(xh , t) and suppose that e(x, 0) = 0. Then 1. For any t > 0 (124)

|e(x, τ )|t ≤

x − xl xh − x |gh |t + |gl |t xh − xl xh − xl

2. If in addition gl (t) ≥ gl (Te ) ≥ 0 and gh (t) ≥ gh (Te ) ≥ 0, for all t ≥ Te . Then for any η > 0 there exists a number Tη such that for all t ≥ Tη + Te   x − xl xh − x (125) e(x, t) ≥ (1 − η) gh (Te ) + gl (Te ) . xh − xl xh − xl Proof. The proof follows from formula (42) and Corollary 3.

 

6.2 Convergence analysis If T is sufficiently large, the error initially decays at a linear rate. Moreover, the number of iterations in the linear regime is a linear function of T . Proof of Theorem 3: We derive the upper bound for the error using an inductive construction analogous to that of Sect. 5. In this construction we use inequality (120) instead of the right inequality in (58). This yields the right inequality in (27). In order to derive a lower bound we define  ηˆ = 1 − 1 − η, (126) and Tη such that inequality (125) holds with ηˆ substituted for η. We first consider the case of two subdomains. We show by induction that for t ≥ Te + (2k − 1)Tη   1 − ∆/S 2(k−1) x − x0 (1 − ηˆ)2k−1 erk (x, t) ≥ S+∆ 1 + ∆/S (127) ×(1 − ν) e0 t , using Corollary 4 and assumption (26). Combining this result with definition (126) and the inequality (128)

erk (x1 , t) ≤ ek t

yields the left inequality in (27) for N = 2. The analysis for three subdomains is analogous.  

Space-time domain decomposition for parabolic problems

303

The upper bound agrees with the result of Gander and Stuart [17] for the infinite time window. We now analyze the case of an arbitrary number N of subdomains. We derive a linear lower bound which is valid for a finite number of iterations. Proof of Theorem 4: We define the hat function  2x−S S/2 ≤ x ≤ N S/2  LS (1 − ν) e0 t     φ0 (x, N ) = (2L+1)S−2x (1 − ν) e0 t N S/2 ≤ x ≤ N S − S/2 LS      0 elsewhere, (129) with L = N − 1. We also define ηˆ and Tη as in Theorem 3. We first assume that N is odd. In this case the sub-domain [xL/2 , xL/2+1 ] × [0, T ] and the sub-domains near the boundaries x = x0 and x = xN are red. Then we note that (130)

eb0 (x, t) ≥ φ0 (x − x0 , N ), t ≥ Te ,

in view of assumption (26). Moreover φ0 (x − x0 , N ) is convex in all subdomains, with the exception of the red sub-domain [xL/2 , xL/2+1 ] × [0, T ]. When combined with Corollary 4 the convexity of φ0 yields (131)

er1 (x, t) ≥ (1 − ηˆ)φ0 (x − x0 , N ), t ≥ Te + T η ,

for x ∈ [xL/2 , xL/2+1 ]. An explicit calculation with Corollary 4 yields   1 + 2∆/S φ0 (x − x0 , N ), er1 (x, t) ≥ (1 − ηˆ) 1 − L t ≥ Te + Tη , (132) for x ∈ [xL/2 , xL/2+1 ]. This inequality is valid in all red subdomains, in view of inequality (131). Then, we use Corollary 4 and the convexity of φ0 in all black subdomains to obtain   1 + 2∆/S φ0 (x − x0 , N ), eb1 (x, t) ≥ (1 − ηˆ)2 1 − L t ≥ Te + 2Tη . (133) Repeated application of the above arguments and the use of definitions (129) and (126) yields inequality (29). When N is even we use the same arguments as above with φ(x − x0 , N − 1) substituted for φ(x − x0 , N ), thereby reducing the problem to the previous case.  

304

E. Giladi, H.B. Keller

Fig. 7. The decomposition of space into three vertical strips. The dark areas represent the overlap of a red sub-domain over a black one. The light grey areas represent the overlap of the black sub-domain over a red one

Fig. 8. The decomposition of space into three vertical strips each of which is then decomposed into three cubes. The shaded areas represent the overlap between subdomains

Space-time domain decomposition for parabolic problems

305

We note that inequality (124) implies that the magnitude of the error in this iteration is bounded by e0 t . In conjunction with Theorem 4 this provides a linear lower and upper envelop with convergence factor that is arbitrarily close to 1 as the number of subdomains N → ∞. 7 The two dimensional advection diffusion equation We now analyze the space-time domain decomposition iteration for equations (5)–(7) in W0 = [xl , xh ] × [yl , yh ] × [0, T ]. We consider a decomposition of space into strips of width S, as depicted in Fig. 7, for three subdomains. We note that each strip may be further decomposed into yellow and green cubes, as depicted in Fig. 8. In this case, we solve the problem in green subdomains of red strips followed by the yellow ones. Then we perform a green yellow sweep in black strips.

7.1 Error analysis for a decomposition into strips We now derive equations for the error in the iteration, for the case of two subdomains. Hence, Ω = [x0 , x2 ] × [yl , yh ] and the red and black subdomains are [x0 , x1 ]×[yl , yh ]×[0, T ] and [x1 , x2 ]×[yl , yh ]×[0, T ], respectively. We proceed as in Sect. 4 to obtain an equation for the error in the red sub-domain and it’s overlap (134)

∂erk = a2 ∆erk + b[sin θ, cos θ] · ∇erk , ∂t

(135)

erk (x, y, 0) = 0,

(136) (137)

erk (x0 , y, t) = erk (x, yl , t) = erk (x, yh , t) = 0, erk (x1 + ∆, y, t) = ebk−1 (x1 + ∆, y, t).

In the black sub-domain and it’s overlap the error satisfies (138)

∂ebk = a2 ∆ebk + b[sin θ, cos θ] · ∇ebk , ∂t

(139)

ebk (x, y, 0) = 0,

(140) (141)

ebk (x2 , y, t) = ebk (x, yl , t) = ebk (x, yh , t) = 0, ebk (x1 − ∆, y, t) = erk (x1 − ∆, y, t).

306

E. Giladi, H.B. Keller

We solve equations (134)–(141) by separation of variables [23] to obtain erk (x, y, t)

=

ebk (x, y, t) =

(142) where

∞  n=1 ∞  n=1

r Mk,n (x, t)Nn (y), b Mk,n (x, t)Nn (y),



b cos θ 2 e− 2a2 y sin (ωn (y − yl )) , yh − yl πn . ωn = yh − yl

Nn (y) = (143)

r and M b satisfy The Fourier coefficients Mk,n k,n r ∂Mk,n

∂t (144)

= a2

r ∂ 2 Mk,n

∂x2 n = 1, 2, . . .

+ b sin θ

r ∂Mk,n

∂x

r − Cn2 Mk,n ,

r r Mk,n (x0 , t) = Mk,n (x, 0) = 0, r b (x1 + ∆, t) = Mk−1,n (x1 + ∆, t), Mk,n

(145)

b ∂Mk,n

∂t (146)

2

=a

b ∂ 2 Mk,n

∂x2 n = 1, 2, . . .

+ b sin θ

b ∂Mk,n

∂x

b , − Cn2 Mk,n

b b Mk,n (x2 , t) = Mk,n (x, 0) = 0,

(147)

b r (x1 − ∆, t) = Mk,n (x1 − ∆, t) Mk,n

(148)

b M0,n (x, t) = M0,n (x, t),

where (149)

 M0,n (x, t) =

yh

yl

x1 ≤ x ≤ x2 ,

p(y)e0 (x, y, t)Nn (y)dy,

e0 (x, y, t) is the initial error,     b cos θ 2 πna 2 2 Cn = (150) + , 2a yh − yl and (151)

p(y) = e

b cos θy a2

.

Space-time domain decomposition for parabolic problems

307

The original problem (134)–(141) has been transformed into a sequence of problems (144)–(147) that govern the evolution of each mode of the error in the iteration. We denote the coefficient of the n’th mode of the error at iteration k by Mk,n (x, t)  r Mk,n (x, t) x ∈ [x0 , x1 ] Mk,n (x, t) = (152) b (x, t) x ∈ [x , x ] , Mk,n 1 2 and the error at iteration k by ek (x, y, t) (153)

ek (x, y, t) =

∞ 

Mk,n (x, t)Nn (y).

n=1

When the number of subdomains N ≥ 3, we proceed in an analogous fashion to obtain for each mode of the error a recursive system of N differential equations, coupled through their boundary conditions. 7.2 Convergence analysis The convergence behavior of the iteration for the two dimensional model problem is essentially the same as in the one-dimensional case. Indeed, let

· 2,t denote the norm  xh  y h 1/2 (154) p(y)e2 (x, y, t)dxdy

e 2,t = max 0≤τ ≤t

xl

yl

with p defined in (151) then Theorem 5 Consider the space time DD iteration in strips for the advection diffusion equation, problem (5)–(7). Suppose that the initial iterate u0 ∈ C 1 (W0 ). 1. If S > 2∆ then for any 0 ≤ t ≤ T   4(k − 1)2 ∆2 (155) ek 2,t ≤ C exp − (1 + o(1)), k → ∞ a2 t 2. If |b| > 0 then for all t ≥ 0 and k ≥ 1 + [(2∆(b/a2 )π]−1      S|b| 2k−2 ∆|b| + exp − 2 (156) ek 2,t ≤ C D(k) exp − 2 a 2a In these equations (157)

1 D(k) = 1 +  4 π(k − 1)∆|b|/a2

and C is a constant which depends on the initial error.

308

E. Giladi, H.B. Keller

Proof. The recursive system of equations for Mk,n , derived in Sect. 7.1, is identical to the system of equations of Sect. 4. Hence Theorems 1 and 2 apply to each mode separately, and they yield   4(k − 1)2 ∆2

Mk,n t ≤ M0,n t C1 exp − a2 t (158) ×(1 + o(1)),

Mk,n t ≤ M0,n t C1 D(k)      ∆|b| S|b| 2k−2 × exp − 2 + exp − 2 (159) , a 2a where (160)

γ

2

 =

b 2a

2

 +

nπa yh − yl

2

,

C1 = e

|b sin θ|(xN −x0 ) 2a2

.

Inequality (159) follows from inequality (21) in view of the relation γ/a ≥ |b|/(2a2 ), which is valid for all n. Moreover, the o(1) term in inequality (158) is independent of n. This is shown by inspecting the proof of Theorem 1 and using the left bound for γ/a nπ γ π ≤ ≤ . yh − yl yh − yl a

(161)

We now use the representation (153) for ek , and the orthonormality of {Nn (y)}∞ n=1 with respect to the inner product  yh (162) p(y)f (y)g(y)dy, (f, g) = yl

(see: [24]) to obtain (163)

ek 2,t ≤

√

xh − xl

∞ 

n=1

1/2

Mk,n 2t

.

Then we note that

(164)

∞ 

1/2 2

( M0,n t )

n=1

where

√ 3/2

(165) C2 = (yh − yl )

π ≤ √ C2 , 6

 b cos θy 2 max e 2a2 π W0

   ∂e0 b cos θ     ∂y + 2a2 e0  .

Space-time domain decomposition for parabolic problems

309

We obtain this inequality by integrating (149) for M0,n by parts to obtain C2 n This manipulation is justified since uy has a bound uniform in y in the interior of W0 and we can choose u0 to have such a bound as well, then e0,y is uniformly bounded in the interior of W0 . The right side of inequality (166) is independent of both x and t and this yields the bound (164). Combining the bounds (158) and (159) for Mk,n t with inequalities (163) and (164) proves the Theorem.   |M0,n (x, t)| ≤

(166)

8 Conclusions and future directions Space-time DD performs best on equations with a small diffusion coefficient and advection dominated equations. The super-linear convergence rate, presented in Theorems 1 and 5 is strongly related to the Gaussian nature of diffusion. This leads us to believe that analogous results can be derived for equations with variable coefficients and possibly non-linear problems, in general domains. Acknowledgements. We wish to thank the anonymous referees for carefully reading the manuscript and for their comments. This work was supported in part by NSF, under cooperative agreement # CCR-9120008.

A Solution to error equation In this appendix we derive the solution to equation ∂2e ∂e ∂e = a2 2 + b − c2 e ∂t ∂x ∂x in [xl , xh ] × [0, ∞) subject to the conditions (167)

(168)

e(x, 0) = 0,

e(xl , t) = gl (t),

e(xh , t) = gh (t).

We use the Laplace transform method (see: [23]) to obtain the equivalent problem b s + c2 Φxx + 2 Φx − (169) Φ = 0, a a2 (170) Φ(xl , s) = Gl (s), Φ(xh , s) = Gh (s),  t  t −sτ Φ(x, s) = e e(τ )dτ, Gl (s) = e−sτ gl (τ )dτ, 0 0  t (171) e−sτ gh (τ )dτ. Gh (s) = 0

310

E. Giladi, H.B. Keller

The solution to equations (169)–(170) is sinh(β(xh − x)) sinh(βL) sinh(β(x − xl )) +Gh (s)eα(x−xh ) sinh(βL)

Φ(x, s) = Gl (s)eα(x−xl ) (172)

(173)

α=−



b , 2a2

β=

b2 s + c2 + 4a4 a2

1/2 .

We use the substitution

  e−βL eβu − e−βu sinh(βu) = , sinh(βL) 1 − e−2βL

(174)

in this expression and expand 1/(1 − exp(−2βL)) in Taylor series to obtain Φ(x, s) = Gl (s)eα(x−xl )

∞  

exp[−βλn (xh − x)]

n=0

 − exp[−βµn (xh − x)] + Gh (s)eα(x−xh ) (175)

×

∞ 

{exp[−βλn (x − xl )] − exp[−βµn (x − xl )]} ,

n=0

where the functions λn and µn are defined in (45). In order to determine e we now compute the inverse transform of each term in the sums (175) using the identities (see: [21])    √  k k2 −1 (176) L exp(−k s) = √ 3/2 exp − 2 , k > 0, 4t 2 πt and (177)

L−1 {G(s + h)} = e−ht L−1 {G(s)}

to obtain with η = λn or η = µn (178)

L−1 {exp (−βη)} = F (t, η, a, γ).

Here L−1 denotes the inverse Laplace transform and (179) (180)

b2 , 4a2   η2 2 exp − 2 − γ t . 4a t

γ 2 = c2 + η F (t, η, a, γ) = √ 3/2 2 πat

Space-time domain decomposition for parabolic problems

311

We now use identity (178) and the convolution Theorem for the Laplace transform in expression (175) to obtain e(x, t) = eα(x−xl ) K(t, xh − x, a, b, c) ∗ gl (t) +eα(x−xh ) K(t, x − xl , a, b, c) ∗ gh (t),

(181) where (182)

K(t, u, a, b, c) =

∞ 

{F [t, λn (u), a, γ] − F [t, µn (u), a, γ]}

n=0

and γ is defined in (179). B A qualified maximum principle for the error equation We now show that solutions to (183) below ∂2e ∂e ∂e = a2 2 + b − c2 e, ∂t ∂x ∂x

(183)

satisfy a maximum principle in [xl , xh ] × [0, T ] when the data is nonnegative. Specifically we assume that e ≥ 0 on (184)

∂ˆ = ({xl } × [0, T ]) ∪ ({xh } × [0, T ]) ∪ ([xl , xh ] × {0})

and we show that ˆ 1. The maximizer of e ∈ ∂. 2. The function e is non-negative throughout [xl , xh ] × [0, T ]. We use similar arguments as in [23] for the heat equation. Let (x0 , t0 ) be ˆ Then at the maximizer of e and suppose by contradiction that (x0 , t0 ) ∈ ∂. (x0 , t0 ) ∂e ∂e ∂2e ≤ 0, (185) = 0, ≥ 0, e ≥ 0. ∂x ∂x2 ∂t Let (186) ψ(x, t, ) = e(x, t) − (t − t0 ). If  is sufficiently small, the maximizer of ψ, (x1 , t1 ), is sufficiently close to (x0 , t0 ) such that (x1 , t1 ) ∈ ∂ˆ and e(x1 , t1 ) ≥ 0. Moreover at (x1 , t1 ) (187)

∂2ψ ≤ 0, ∂x2

∂ψ = 0, ∂x

∂ψ ≥0 ∂t

and therefore at (x1 , t1 ) (188)

∂2e ≤ 0, ∂x2

∂e = 0, ∂x

∂e > 0, e ≥ 0. ∂t

312

E. Giladi, H.B. Keller

ˆ This contradicts (183) and proves that (x1 , t1 ) ∈ ∂. Let (x0 , t0 ) denote the minimizer of e and suppose that e(x0 , t0 ) < 0. Then (x0 , t0 ) ∈ ∂ˆ and therefore (189)

∂e = 0, ∂x

∂2e ≥ 0, ∂x2

∂e ≤ 0, e < 0, ∂t

there. This contradicts (183) and proves that e ≥ 0 throughout [xl , xh ] × Ω. References 1. E. Lindel¨of: Sur l’application des m´ethodes d’approximations successives a´ l’´etude de ceratines equations diff´erentielles ordinaires. Journal de Mathematiques Pures et Appliqueees 9, 217–271 (1893) 2. E. Lindel¨of: Sur l’application des m´ethodes d’approximations successives a´ l’´etude des int´egrales r´eeles des equations diff´erentielles ordinaires. Journal de Mathematiques Pures et Appliqueees 10, 117–128 (1894) 3. E. Lelarasmee, A.E. Ruehli, A.L. Sangiovanni-Vincentelli: The waveform relaxation method for time-domain analysis of large scale integrated circuits. IEEE Trans. Computer-Aided Design 1, 131–145 (1982) 4. A.R. Newton, A.L. Sangiovanni-Vincentelli: Relaxation-based electrical simulation. SIAM J. Sci. Stat. Computing 4(3), 485–524 (1983) 5. J.K. White, A.L. Sangiovanni-Vincentelli: Relaxation Techniques for the Simulation of VLSI Circuits. Kluwer Academic Publishers, Boston, 1987 6. Ch. Lubich, A. Ostermann: Multigrid dynamic iteration for parabolic equations. BIT 27, 216–234 (1987) 7. U. Miekkala, O. Nevanlinna: Convergence of dynamic iteration methods for initial value problems. SIAM J. Sci. Stat. Computing 8, 459–482 (1987) 8. U. Miekkala, O. Nevanlinna: Sets of convergence and stability regions. BIT 27, 554–584 (1987) 9. O. Nevanlinna: Remarks on Picard-Lindelof iteration. part i. BIT 29, 328–346 (1989) 10. O. Nevanlinna: Remarks on Picard-Lindelof iteration. part ii. BIT 29, 535–562 (1989) 11. U. Miekkala: Dynamic iteration methods applied to linear DAE systems. J. Comp. Appl. Math. 25, 133–151 (1989) 12. C.W. Gear: Waveform methods for space and time parallelism. J. Comp. Appl. Math. 38, 137–147 (1991) 13. S. Vandewalle: Parallel Multigrid Waveform Relaxation for Parabolic Problems. B.G. Teubner Verlag, Stuttgart, 1993 14. S. Vandewalle, R. Piessens: On dynamic iteration methods for solving time-periodic differential equations. SIAM Numer. Anal. 30, 286–303 (1993) 15. M. Bjørhus: A note on the convergence of discretized dynamic iteration. BIT 35, 291– 296 (1995) 16. R.Jeltsch, B.Pohl: Waveform relaxation with overlapping splittings. SIAM J.Sci. Comput., 16(1):40–49, 1995 17. M.Gander, A.Stuart: Space-time continuous analysis of waveform relaxation for the heat equation. SIAM J.Sci. Comput. 19(6), 2014–2031 (1998) 18. F.Wang, J.Xu: A crosswind block iterative method for convection-dominated problems. SIAM J.Sci. Comput. 21(2), 620–645 (1999)

Space-time domain decomposition for parabolic problems

313

19. J.M.Fiard, T.A.Manteuffel, S.F.Mccormick: First-order system least squares (fosls) for convection-diffusion problems: numerical results. SIAM J.Sci. Comput. 19(6), 1958– 1979 (1998) 20. Eldar Giladi, Herbert B. Keller: Space-time domain decomposition for parabolic problems. Technical report, CRPC, Caltech, March 1997 21. M.Abramowitz, A.Stegun: Handbook of Mathematical Functions. Dover, 1970 22. N.J.DeBruijn: Asymptotic methods in analysis. Dover, 1981 23. G.F.Carrier, C.E.Pearson: Partial Differential Equations, Theory and Technique. Academic Press, 1988 24. R.V.Churchill: Four ier series and boundary value problems. McGraw-Hill, 1941