Domain Decomposition Operator Splittings for the Solution of Parabolic Equations

We study domain decomposition counterparts of the classical alternating direction implicit (ADI) and fractional step (FS) methods for solving the large linear systems arising from the implicit time stepping of parabolic equations. In the classical ADI and FS methods for parabolic equations, the elliptic operator is split along coordinate axes; they yield tridiagonal linear systems whenever a uniform grid is used and when mixed derivative terms are not present in the differential equation. Unlike coordinate-axes-based splittings, we employ domain decomposition splittings based on a partition of unity. Such splittings are applicable to problems on nonuniform meshes and even when mixed derivative terms are present in the differential equation and they require the solution of one problem on each subdomain per time step, without iteration. However, the truncation error in our proposed method deteriorates with smaller overlap amongst the subdomains unless a smaller time step is chosen. Estimates are presented for the asymptotic truncation error, along with computational results comparing the standard Crank–Nicolson method with the proposed method.


Introduction.
The numerical solution of parabolic partial differential equations by implicit time stepping procedures requires the solution of large systems of linear equations. These linear systems need only be solved approximately, provided that the inexact solutions obtained by using approximate solvers preserve the stability and local truncation error of the original scheme. Though iterative methods with preconditioners are a popular way for solving such linear systems, see [5,6,3,8,15,33,4,7], in this paper we consider only approximate solvers that do not involve iteration at each time step. Such approximate noniterative solution methods for parabolic equations include the alternating direction implicit (ADI) methods of Peaceman and Rachford [22] and Douglas and Gunn [12], the fractional step methods (FS) of Bagrinovskii and Godunov [1], Yanenko [34], and Strang [26], and also more recent one-iteration domain decomposition solvers of Kuznetsov [16,17], Meurant [21], Dryja [13], Blum, Lisky, and Rannacher [2], Dawson, Du, and Dupont [9], Laevsky [19,18] and Vabishchevich [28] and Vabishchevich and Matus [29], and Chen and Lazarov [20].
The method proposed here uses the same framework as the classical ADI and FS methods for solving parabolic equations of the form u t + Lu = f using an operator splitting L = L 1 + · · · + L q ; however, the splittings chosen here are based on domain decomposition, unlike classical splittings along coordinate directions. Furthermore, they are applicable to problems with mixed derivative terms and on nonuniform grids. The basic idea is simple. Given a smooth partition of unity {χ k } k=1,...,q subordinate to a decomposition of the domain, as described in section 4, an elliptic operator L, such as the two-dimensional Laplacian, is split as a sum of several "simpler" operators: where L i u = −∇ · χ i ∇u.
Given a splitting, an approximate solution of the parabolic equation is obtained by solving several "simpler" evolution equations of the form u t + L k u = f k ; see section 3. By comparison, classical splittings, see Richtmyer and Morton [23], of L = −Δ are of the form L 1 = − ∂ 2 ∂x 2 and L 2 = − ∂ 2 ∂y 2 . Our goal in this paper is to derive truncation error estimates of the proposed methods and to discuss their stability properties. The rest of the paper is outlined as follows. In section 2, we describe a parabolic equation and its discretization. In section 3, we describe the classical ADI and FS methods in matrix language. In section 4, we describe domain decomposition operator splittings using a partition of unity and derive asymptotic truncation error bounds and stability properties for the proposed method. A heuristic comparison is then made between the Crank-Nicolson method and the ADI method with domain decomposition splitting. Finally, in section 5, numerical results are presented.

A parabolic equation and its discretization.
We consider the following parabolic equation for u(x, y, t) on a domain Ω ⊂ R 2 : where a(x, y) is a 2 × 2 symmetric positive definite matrix function with C 2q−1 (Ω) entries, c(x, y) ≥ 0 is in C 2q (Ω), f (x, y, t) is a C 2q (Ω×[0, T ]) forcing term, and u 0 (x, y) is a C 2q (Ω) initial data. Here q is an integer to be specified later in section 4. With these smoothness assumptions, the exact solution u(x, y, t) will be in C 2q (Ω × [0, T ]). We will denote the elliptic operator by Lu = −∇ · a(x, y)∇u + c(x, y)u, and the parabolic equation by u t + Lu = f .

Discrete problem.
Let A h denote the symmetric positive definite matrix obtained by second-order finite difference discretization of the elliptic operator L on a grid on Ω with mesh size h and grid points (x i , y j ): where U ij ≈ u(x i , y j , t) denotes the discrete solution. To discretize in time, we apply an implicit scheme such as the backward Euler or Crank-Nicolson scheme, which results in a large sparse linear system to be solved.
For example, the backward Euler method yields where τ is the time step, U n = {U (x i , y j , nτ)} denotes the vector of discrete unknowns at time t n = nτ , and F n = {f (x i , y j , nτ)} is the discrete forcing term at time t n . As is well known, see [23], the above scheme is unconditionally stable in the Euclidean norm . and convergent with error: u − U be = O (τ ) + O h 2 . Thus, for transient problems τ = O(h 2 ) provides an optimal choice of time step τ for the backward Euler method.
Similarly, the Crank-Nicolson method yields which requires the solution of the linear system The Crank-Nicolson method is also unconditionally stable, see [23], with error: We note that both of the above discretizations can be written in the form for suitable choices of constants α, matrices B, and forcing terms g n . In particular, the coefficient matrices can be written in the form I + ατ A h for α = 1 or α = 1 2 , and they are symmetric positive definite. Furthermore, since it is known that the largest , it follows that the condition number of I + ατ A h is bounded by 1 + cτ h −2 , for some positive constant c independent of τ and h. Thus, for h 2 τ 1 systems (1) and (3) are ill conditioned, yet better conditioned than A h . Consequently, if iterative methods are used, preconditioners such as multigrid [4,15] or domain decomposition methods may be needed to reduce the number of iterations; see, for example, [5,6,3,8,33,7].
However, as mentioned earlier, system (4) for U n+1 needs only to be solved approximately within the local truncation error. In this paper we only seek solvers that provide approximate solutions at the cost of one iteration per time step. Such approximate solvers will modify the scheme (4) and alter its stability and truncation error. In order to distinguish this with the modified schemes, we will henceforth refer to scheme (4) as the original scheme. In sections 3 and 4, we discuss one-step approximate solvers, and estimate their truncation errors using the framework of the classical ADI and FS methods. We also describe conditions under which the stability of the original scheme is preserved.
Given the original discretization (4), we define its local truncation error T orig as follows: where u is the exact solution of the parabolic equation restricted to the grid points. Discretization (4) is said to be stable in a norm · if the following holds for some numbers c 1 , c 2 : If c 1 and c 2 are independent of h and τ , then the discretization is said to be unconditionally stable, else conditionally stable.

ADI and FS methods.
Given an evolution equation of the form u t +Lu = f , the basic idea in the ADI and FS methods is to approximately solve it using an operator splitting L = L 1 + · · · + L q , where several simpler evolution equations of the form w t + L i w = f i are solved with suitably chosen f i to provide an approximate solution w ≈ u to the original equation.
We will henceforth restrict our discussion to the matrix case, where a discretization in space of u t + Lu = f yields U t + A h U = F h , and a subsequent discretization in time yields system (4). Corresponding to an operator splitting L = L 1 + · · · + L q , where q > 1, we assume a matrix splitting of A h as a sum of symmetric positive semidefinite matrices A i : It will be assumed that the matrices A i have a "simpler" structure than A h , in the sense to be explained below.
In matrix language, given the splitting {A i }, both the ADI and FS methods approximately solve systems of the form (I + ατ A h )w = z by solving several systems of the form (I + ατ A i )w i = z i ; see sections 3.1 and 3.2. Hence, the A i must be suitably chosen so that the latter systems are easier to solve than the original system.
As an example of splittings, we briefly describe the classical splittings, see [23], of an elliptic operator along the x and y coordinate directions (when no mixed derivative terms are present). Domain decomposition splittings are described in section 4. For classical splittings, the two-dimensional Laplacian L = −Δ is split as a sum of L 1 = − ∂ 2 ∂x 2 and L 2 = − ∂ 2 ∂y 2 . In matrix terms where Δ h denotes the discrete Laplacian and D h xx and D h yy denote the discretization of ∂ 2 ∂x 2 and ∂ 2 ∂y 2 , respectively. Both A 1 = −D h xx and A 2 = −D h yy will be tridiagonal matrices after a suitable ordering of the unknowns (if the grid is uniform) and hence (I + ατ A i ) will also be tridiagonal and solvable in linear time on uniform grids! Unfortunately, such splittings are not possible if mixed derivative terms are present and tridiagonality can be lost if the grid is nonuniform. By contrast, the domain decomposition splittings described in section 4 are applicable in more general cases, though the resulting linear systems are not tridiagonal.

ADI method.
This method was originally proposed by Peaceman and Rachford [22] and Douglas [11,10], and was further developed by several authors. We follow here the general formulation described by Douglas and Gunn [12], where additional references may be found. The reader is also referred to [30,31,32,14,24,27]. In matrix language, given a splitting A h = A 1 + · · · + A q and an approximate solution W n ≈ U n to the discrete solution at time t n , an approximate solution W n+1 ≈ U n+1 of the linear system 3. Define the ADI approximate solution at time t n+1 = (n + 1)τ to be We note that in order to compute W n+1 in the ADI scheme, we need to solve a total of q linear systems, each of the form I + ατ A i of size n × n (the size of A h ). The following stability results for the ADI method can be found in [12]. For q ≥ 3 in the noncommuting case, examples are known of positive semidefinite splittings for which the ADI method can loose unconditional stability; see [12]. Nevertheless, the above stability results are a bit pessimistic and instability is only rarely encountered in practice.
Next, we describe the truncation error terms for the ADI method, derived in [12]. The magnitude of these truncation terms will be estimated in section 4.3 for domain decomposition splittings. THEOREM 3.2. Given the discrete problem with local truncation error T orig (see (5)), the ADI solution W n+1 solves the following perturbed problem: which introduces the following additional terms in the local truncation error: where u is the exact solution of the parabolic equation restricted to the grid points.
Proof. The proof is given in [12]; however, for completeness we include it here. From step 2 of the ADI algorithm, we obtain that W n+1 By recursively applying this, we deduce that Now, from step 1 of the ADI algorithm, we deduce that Multiplying (9) by (I + ατ A 1 ) and substituting (10) for the resulting right-hand side, we obtain Noting that Explicit estimates of these truncation error terms for the case of domain decomposition splittings are given in section 4.3.

FS methods.
These methods were originally developed by Bagrinovskii and Godunov [1] and Yanenko [34] and further developed by Strang [26] and other authors. Given a matrix splitting A h = A 1 + · · · + A q , the FS methods approximate (I + ατ A h ) by a product of several matrices of the form (I + ατ A i ). For example, the first-order FS method uses the approximation which can be formally verified by multiplying all the terms, bearing in mind that the matrices A i may not commute. Given (12), an approximate solution W n+1 ≈ U n+1 of can be obtained by solving We summarize below the first-order fractional step algorithm to approximately solve 1. Let z 1 = g n − BW n and for i = 1, . . . , q solve (14) is a second-order approximation (locally) in τ of (13), it is only firstorder accurate globally, due to accumulation of errors; see [23]. The following result concerns the stability of the first-order FS method. THEOREM 3.3. If A i are symmetric and positive semidefinite, and if B ≤ 1 in the Euclidean norm . , then the first-order FS method is stable with W n+1 ≤ W n + c g n for some c > 0 and independent of τ and h.
Proof. Since W n+1 = (I +ατ A q ) −1 · · · (I +ατ A 1 ) −1 (g n − BW n ) and since B ≤ 1 in the spectral norm, we only need verify that in the spectral norm. But this follows easily: since each of the terms (I + ατ A i ) are symmetric positive definite with eigenvalues greater than one, the spectral norms of their inverses are each bounded by one.
The following result concerns the local truncation error of the first-order FS method.
THEOREM 3.4. The first-order FS approximate solution W n+1 of The local truncation error T fs of the first-order FS method has the following terms in addition to the terms in the local truncation error T orig when exact solvers are used: Here u n+1 is the exact solution of the parabolic equation at time (n + 1)τ restricted to the grid points.
Proof. Since the product the proof immediately follows by replacing (I + ατ A h ) by (I + ατ A 1 ) · · · (I + ατ A q ) in the original discretization. For a given splitting {A i }, the above local truncation error terms are O(τ 2 ) (though the constant in the leading-order term τ 2 can vary for different choices of splittings). However, the global error will be O(τ ), see [23], due to accumulation of errors, and consequently, the above first-order FS method is suitable only for globally first-order schemes such as backward Euler but not for globally second-order discretizations such as Crank-Nicolson.
For globally second-order accurate schemes such as Crank-Nicolson, second-order fractional step approximations can be obtained by using Strang splitting [26]. We outline the main idea for a splitting involving two matrices: The aim is to approximate e −ατ A h by a locally third-order accurate approximation as follows: Then, as in the Crank-Nicolson method, a locally third-order approximation in τ of the terms e −ατ Ai can be used to replace each of the exponentials Applying this to each of the terms in (17) yields the implementation of Strang splitting. We state without proof the following. THEOREM 3.5. Strang splitting provides unconditionally stable approximations when A i are symmetric positive semidefinite: A i ≥ 0. For any fixed choice of positive semidefinite matrices A 1 and A 2 , the truncation error for Strang splitting is thirdorder accurate locally in τ , and globally second-order accurate in τ .
Since Strang splitting provides a globally second-order accurate approximation in time, it can be applied to inexactly solve the Crank-Nicolson scheme. However, it is more expensive and complicated to generalize when there are more than two matrices in the splitting, and it requires more linear systems to be solved than the first-order FS method. In the applications we consider, the ADI methods which were described earlier are preferable and generate smaller truncation errors, though they are not guaranteed to be unconditionally stable without additional assumptions.
As noted in the introduction, given a smooth partition of unity {χ k (x, y)} k=1,...,q subordinate to a collection of open subdomains {Ω * k } k=1,...,q which cover Ω, we split an elliptic operator L as a sum of several "simpler" operators as follows: In what follows, we will denote the coefficients of L k by a k (x, y) = χ k (x, y)a(x, y) and c k (x, y) = χ k (x, y)c(x, y). For numerical purposes, we could replace a smooth partition of unity by one which is piecewise smooth. In such a case, it may be necessary to define the operators L k weakly using a variational approach. The matrices A k are obtained by discretization of L k .

Partition of unity. Given an overlapping collection of open subregions
..,q whose union contains Ω, a partition of unity subordinate to this covering is a family of q nonnegative, C ∞ (Ω) smooth functions {χ k (x, y)} k=1,...,q whose sum equals one and such that χ k (x, y) vanishes outside Ω * k : For numerical purposes, we will consider χ k (x, y) which are not necessarily C ∞ (Ω) but are continuous and piecewise smooth. We will refer to such a partition of unity as piecewise smooth. For more information on partitions of unity, see Sternberg [25].
For example, the subdomains Ω k can be chosen to be the coarse elements in a coarse grid triangulation of Ω with mesh size H; see Figure 1. 2. Next, enlarge each subdomain Ω k to Ω * k to contain all points in Ω within a distance β of Ω k , where β > 0; see Figure 1. Then {Ω * k } will form an overlapping covering of Ω: Ω = ∪ q k=1 Ω * k . We note here that the subdomains {Ω * k } can often be grouped into a small number of colors so that any two subdomains of the same color are disjoint, see Figure 1, for a coloring of 16 overlapping subdomains into four colors (q = 4, only color 1 is indicated).
The above constructed {χ k (x, y)} will be continuous and piecewise smooth (with discontinuous derivatives across ∪ k ∂Ω * k ). To obtain C ∞ (Ω) smooth partitions of unity, we refer the reader to [25].
Before we proceed further, we define the Hölder norms for integer k ≥ 0 as follows: where we have used the multi-index notation For smooth partitions of unity with overlap β, the following result holds.  (19) in the Hölder C m (Ω) norm, where c m > 0 is a constant independent of k and β but dependent on the geometry of the region.
3. By repeated applications of the chain rule, it follows that Using the definition of the Hölder norm it follows that where c m is defined by c m ≡ χ k C m (Ω) . We omit further details.
For numerical implementation, smooth partitions of unity can be replaced by piecewise smooth ones based on the distance functions w k (x, y). See section 5 for an example of an alternative partition of unity for rectangular subdomains.

Domain decomposed matrix splittings.
As mentioned before, given a partition of unity {χ k } k=1,...,q as in section 4.1, a domain decomposition splitting of the operator Lu = −∇ · a(x, y)∇u − c(x, y)u is obtained as L k u = −∇ · (a k (x, y)∇u) + c k (x, y)u(x, y). By construction, q k=1 L k = L, and each L k are self-adjoint. A k are then defined as suitable discretizations of the differential operators L k , and so A h = q k=1 A k . We note that if c(x) ≥ 0, then each L k and A k will be positive semidefinite. The A k constructed above can then be used in the ADI and FS methods described in section 3, where the solution of (I +ατ A k )w = z corresponds to solving a problem on subregion Ω * k . Additionally, by coloring, problems on disjoint subregions of the same color can be solved in parallel.

Truncation errors of the ADI and FS methods for domain decomposition splittings. As described in sections 3.1 and 3.2, when the ADI and FS methods are used to approximately solve a discretization
the resulting approximate solutions W n+1 adi and W n+1 fs solve the perturbed problems (7) and (15), respectively. These perturbed problems correspond to new discretizations of the original parabolic equation u t + Lu = f with additional truncation terms of the form τ m A σ1 · · · A σm (W n+1 − W n ) and τ m A σ1 · · · A σm W n+1 for the ADI and FS methods, respectively. Our goal in this section is to estimate the magnitude of these additional terms for the case of domain decomposition splittings.
As is well known, see for instance Richtmyer and Morton [23], once the local truncation errors T adi and T fs are estimated, the global errors E adi (t) = u(t) − W adi (t) and E fs (t) = u(t) − W fs (t) will follow by Lax's convergence theorem, if the schemes are stable. We note, however, that the global error will be larger than the local truncation errors by a factor of τ −1 due to the accumulation of truncation errors.
Since the local truncation errors T orig for the case of backward Euler or Crank-Nicolson with exact solvers are known, our estimates will focus on the new truncation terms A σ1 · · · A σm W and A σ1 · · · A σm (W n+1 − W n ). For simplicity, we will assume that c(x, y) = 0 and that a(x, y) is a smooth scalar function. The analysis for the cases when c(x, y) ≥ 0 and for smooth matrix functions a(x, y) is similar. Our analysis will require the following smoothness assumptions for the functions a(x, y), W (x, y, t), and {χ k (x, y)}; a(x, y) should be in C 2q−1 (Ω), W (x, y, .) should be in C 2q 0 (Ω), and {χ k (x, y)} should be in C 2q−1 (Ω). With some abuse of notation, W will also refer to W = {W (x i , y j , .)}, the vector of nodal values of W (x, y, .) on the grid points of Ω.
Our estimates will be based on two simple observations, the details of which will be provided subsequently.
1. The finite difference approximations A σ1 · · · A σm W (which occurs in T fs ) and A σ1 · · · A σm (W n+1 −W n ) (which occurs in T adi ) can be shown to be "similar" to their continuous counterparts L σ1 · · · L σm W and L σ1 · · · L σm (W n+1 −W n ), respectively. This result is obtained by repeated applications of the integral form of the mean value theorem to the difference quotients. 2. Using step 1, the expressions A σ1 · · · A σm W and A σ1 · · · A σm (W n+1 − W n ) can be estimated in terms of the maximum values of the derivatives of W and the derivatives of the coefficients {a k (x, y)}. The derivatives of {a k (x, y)} can be estimated in terms of the derivatives of a(x, y) and of the partition of unity functions {χ k (x, y)}. Using the bounds (19) for the partition of unity, we can then obtain that where C(a) is a positive constant that depends on the coefficients a(x, y) and its derivatives. Similar estimates will be valid for A σ1 · · · A σm W n+1 − W n but with an additional factor of τ due to the difference in time W n+1 − W n . Before we proceed, we introduce some notation. In addition to multi-index derivatives, ∂ α = ∂ |α| ∂x α 1 ∂y α 2 , we will use multi-index integrals Iterated integrals are defined as follows: We also define I x h, 1 2 as

Similar integrals are defined in y.
We now outline the details of step 1. By repeatedly applying the integral form of the mean value theorem (or the fundamental theorem of calculus), one can express all the difference quotients in A σ1 · · · A σm W in terms of integrals of derivatives of W and integrals of the derivatives of the coefficients a σ k (x, y). We have the following result. LEMMA 4.2. Let W ∈ C 2m 0 (Ω) and let the expansion of L σ1 · · · L σm W be of the form where each C α1,...,αm+1 = 0 or 1. Then, for W (x, y) the product A σ1 · · · A σm W satisfies where each of the I * h terms (with at most 2m such terms) results from a suitable choice of single or multiple integrals I h .
Proof. The proof will follow by induction on m. For m = 1, the standard five- This follows trivially by repeatedly applying the fundamental theorem of calculus to the difference quotients in the finite difference discretization. We can restate equation (21) using the multi-index derivative and integral notation as But this is just an expression for L k W with multi-index integrals applied to its terms. Consequently, our proof for the case m = 1 is complete. For m ≥ 2, we use the induction hypothesis and assume that the result holds true for m − 1. We will then show that the result holds for m. To do this, define V = A σ2 · · · A σm W . Then, by the induction hypothesis Now, A σ1 V can be calculated by using (22) which involves first-and second-order partial derivatives of V (x, y) locally.
The resulting expression can be simplified by the following observation. The partial derivative operators ∂ α "commute" with all the integral operators I x h , I y h , etc. The commutativity of ∂ x and I y h follows from the differentiation rule for integrals depending on a parameter. The commutativity of ∂ x and I x h can be seen by the following example: When higher-order partial derivatives and iterated integrals are involved, they can be handled by repeated applications of the above. Consequently, all the partial derivatives can be moved inside the integrals. When all the partial derivatives are moved inside the integrals, we obtain an expression similar to L σ1 · · · L σm W but containing integrals of the terms. This gives the desired result for m.
This completes the sketch of the proof of step 1. We obtain step 2 immediately from Lemma 4.2 by applying the mean value property of integrals.
Applying this to all the integrals I h in Lemma 4.2, we obtain that This gives us the desired result. Using Leibnitz's rule for partial differentiation of a k (x, y) = χ k (x, y)a k (x, y) and applying bounds (19) for the partition of unity, we deduce that Here C(a) is a constant that depends only on the coefficients a(x, y) and some of its derivatives.
By substituting these results in Lemma 4.3, we obtain the following.

Similarly, the discretization
This completes the details of step 2. We are finally able to estimate the local truncation errors and the global errors of the ADI method with domain decomposition splitting.
where T orig is the local truncation error of the original scheme (4) when exact solvers are used. Whenever the ADI scheme is stable, the error E adi (t) = u(t) − W adi (t) satisfies Proof. By (8), the local truncation error satisfies Considering the Euclidean norm · of the above expression, and estimating the additional terms using Corollary 4.4, we obtain where n is the number of grid points in Ω. Multiplying the entire expression by h and using that h · is equivalent to · L 2 (Ω) , we can replace h T adi and h T orig by T adi L 2 (Ω) and T orig L 2 (Ω) , respectively. The desired estimate for the local truncation error then follows by noting that √ nh is proportional to the square root of the area of Ω, which can be absorbed in C(a).
The error follows from the local truncation error by an application of the Lax convergence theorem.
Remark 1. Each of the terms (∂ α1 a σ1 ) · · · (∂ αm a σm ) (∂ αm+1 W ) will be zero in most of Ω, more specifically, outside ∩ m k=1 Ω * σ k , and our proof has not taken this into account. This is due to the compact support of the coefficients a σ k (x, y) in Ω * σ k . However, we expect that an estimate taking this into account will be much more complicated. Consequently, our theoretical estimate for the truncation errors are pessimistic when compared with actual numerical results from section 5.
Remark 2. For a fixed choice of subdomains {Ω * k } and partition of unity {χ k }, the overlap β is fixed. Consequently, the above truncation error bounds for the ADI method with domain decomposition splitting is asymptotically second order in time, i.e., as τ → 0 Remark 3. As the number of subdomains increase, and their diameters H and overlap β decrease (i.e., β, H → 0), the accuracy of the proposed ADI method deteriorates. However, since all the terms involving β are multiplied by powers of τ , we can choose a smaller time step τ adi for the ADI method so that the resulting global error is O(h 2 ), the same as for the Crank-Nicolson method with exact solvers. A simple calculation yields that substituting However, since τ adi < h, the ADI method will require more time steps than the Crank-Nicolson method. For a heuristic comparison of the complexity of the ADI method with the standard Crank-Nicolson method, see section 4.4. The truncation errors of the first-order FS method can be analyzed similarly. They are larger than the corresponding errors for the ADI method by a factor of τ −1 .
THEOREM 4.6. The local truncation error T fs of the FS method with domain decomposition splitting satisfies where u is the exact solution of the partial differential equation. The global error satisfies Proof. The proof follows by an application of Lemma 4.3 to the truncation error of the first-order FS method.

A comparison of the Crank-Nicolson method with the proposed ADI method.
In this section, we heuristically compare the work complexity (total floating point operations) of the Crank-Nicolson method using exact solvers with that of the ADI method using domain decomposition splitting. Our comparison will be based on calculating the cost of each method for computing the solution up to time t f and to the same specified accuracy. In order to obtain an explicit comparison, we make the following assumptions: 1. An algebraic direct solver is used to solve the linear systems occurring in both the Crank-Nicolson method and in the ADI method with domain decomposition splitting. The cost of the direct solver to solve a linear system in n unknowns is assumed to be C α n α , where 1 ≤ α ≤ 3. 2. The number of subdomains in the ADI method is chosen to be N d , for some positive integer N , where d is the space dimension (d = 2 or d = 3). The diameter H and overlap β of each subdomain Ω * k is chosen to be proportional to 3. The number of unknowns in each subdomain Ω * k is assumed to be n/N d M , where M > 1 is some magnification factor that enlarges the number of unknowns in each subdomain. 4. The time step τ cn of the Crank-Nicolson method is chosen to be τ cn = h. In order to ensure that both the Crank-Nicolson method and the ADI method have similar errors, see Remark 3 from section 4.3, the time step τ adi of the ADI method is chosen to be For the above choices, both methods have global errors satisfying Finally, since we chose β ∼ N −1 in assumption 2, the ADI time step becomes For both methods, the cost of computing the solution up to time t f satisfies Cost = (Number of time steps) × (Cost per time step) .
For the Crank-Nicolson method we obtain where n is the number of space unknowns. For the ADI method with domain decomposition splitting, we obtain where M > 1 is the constant magnification factor and q is the number of colors in the domain decomposition splitting; see section 4.1. Asymptotically, when 1 N d < n, we obtain that or equivalently if α > 1 + 2q−1 qd . Thus, in two dimensions (d = 2), the ADI method will cost asymptotically less than Crank-Nicolson if the direct algebraic solver costs C α n α with α > 2 − 1 2q . In three dimensions, the ADI method will be preferable if α > 5 3 − 1 3q . Remark 4. Numerical tests in section 5 indicate that a more reasonable time step for the ADI method is τ adi ∼ hβ, which leads to similar global errors for the ADI and Crank-Nicolson methods. For such choices of time steps, the ADI method will be competitive with the Crank-Nicolson method if the complexity of the direct solvers is C α n α where where Ω ⊂ R d .

Numerical results.
We report here on the results of numerical tests conducted using the generalized ADI scheme and the first-order FS method described in section 3, using the domain decomposition splittings of section 4. Our goal in these The standard five-point finite difference scheme was used to discretize in space on a uniform grid of mesh size h. The Crank-Nicolson method was used to discretize in time with time steps τ cn as indicated in the tables. In the tables, we tabulate the scaled Euclidean norm (scaled according to area) of the global error at time t f = 1 for the Crank-Nicolson solution using exact solvers denoted by CN, and Crank-Nicolson with inexact ADI-based solvers denoted by ADI, and Crank-Nicolson with inexact first-order fractional step solvers denoted by FS.
Partitions of unity: The subdomains Ω * k chosen in our experiments were rectangular. For the case of overlapping subregions, we constructed piecewise smooth partitions of unity as follows. For a model rectangle Ω * k = [0, a] × [0, b] we define ω k (x, y) = sin(πx/a) sin(πy/b) inside Ω * k and zero outside it. Then, {χ k (x, y)} is derived from {ω k (x, y)} as described in section 4.1. Table 1 contains the results for the limiting case of zero overlap, i.e., β = 0. This case does not strictly fit in the theoretical framework developed in the paper. We consider two nonoverlapping subdomains Ω * 1 = Ω 1 = [0, 1 2 ] × [0, 1] and Ω * 2 = Ω 2 = [ 1 2 , 1] × [0, 1], and we define χ i (x, y) = χ Ωi (x, y), for i = 1, 2, i.e., the characteristic or indicator functions of the subdomains. In this case, the operators L k can only be defined weakly, in a variational sense. However, the finite difference approximations were well defined without additional modification since the grid we chose was aligned with the interface. Though the overlap β = 0, the numerical results indicate that the domain-decomposition-based ADI solution has, approximately, only thrice as large errors compared to the solution of Crank-Nicolson with exact solvers. The FS solution, however, had about 100 times larger error than the CN solution. Table 2 contains the results for two overlapping subdomains: Ω * 1 = [0, 3 4 ] × [0, 1] and Ω * 2 = [ 1 4 , 1] × [0, 1]. In this case the overlap β = 1/4, and we note that the domain-decomposition-based ADI solution has error very close to the CN solution (using exact solvers). The FS solution had approximately 50 times larger error than this CN solution. Table 3 contains the results for the case of 16 overlapping subdomains, as shown in Figure 1. A 4 × 4 rectangular partition was first formed and each subrectangle was enlarged to include overlap of β = 1/12 (1/3 of the subdomain width). The subdomains were grouped into q = 4 colors, and one of the four colored subregions is shown in Figure 1. In solving (I + ατ A k )w = z, parallelization is possible with different solvers on disjoint subregions of the kth color. The results indicate that the domain-decomposition-based ADI solution has approximately seven times larger error than the CN solution with exact solvers. This was expected due to the larger β −1 term in the truncation error formula of (23). The FS solution had approximately 26 times larger error than this CN solution. Table 4 tabulates the global error for a fixed grid of size h −1 = 64 on which 16, 64, and 256 subdomains were chosen. According to the asymptotic error bounds of equation (23), the global error should remain approximately O(h 2 ) if τ adi = hβ 2q−1 q . We chose τ adi = hβ, which is more optimistic than the theoretical bounds (23). Thus, as the number of subdomains increased, the cost of calculating the solution also increased proportionally to the inverse of the subdomain width. The results indicate that the error in the ADI solution remained approximately constant and close to the error for the CN solution with exact solvers.

Discussion.
The numerical results indicate that in the case of two subdomains with reasonable overlap, the global error of the ADI method with domain decomposition splitting is similar to errors for the Crank-Nicolson method with exact solvers. However, for the case of many subdomains, the terms β −1 in the error (23) causes a deterioration in the global error. If a smaller time step τ adi = hβ is used in the ADI method, then the resulting global error of the ADI solution remains close to the error for CN with exact solvers and τ cn = h.
Our heuristic study of the complexity of these methods in section 4.4 leads us to note that if algebraic direct solvers of complexity C α n α are used to solve all linear systems, then the ADI method with domain decomposition splitting can be competitive with the Crank-Nicolson method provided α > 1 + 2q−1 qd , where d is the space dimension.
Finally, we note that though the ADI method can become unstable, when q ≥ 3 colors are used, such instability was not observed for the range of mesh, subdomain, and time step sizes used in our numerical experiments.