Space-Time Continuous Analysis of Waveform Relaxation for the Heat Equation.
Waveform relaxation algorithms for partial differential equations (PDEs) are traditionally obtained by discretizing the PDE in space and then splitting the discrete operator using matrix splittings. For the semidiscrete heat equation one can show linear convergence on unbounded time intervals and superlinear convergence on bounded time intervals by this approach. However, the bounds depend in general on the mesh parameter and convergence rates deteriorate as one refines the mesh.Motivated by the original development of waveform relaxation in circuit simulation, where the circuits are split in the physical domain into subcircuits, we split the PDE by using overlapping domain decomposition. We prove linear convergence of the algorithm in the continuous case on an infinite time interval, at a rate depending on the size of the overlap. This result remains valid after discretization in space and the convergence rates are robust with respect to mesh refinement. The algorithm is in the class of waveform relaxation algorithms based on overlapping multisplittings. Our analysis quantifies the empirical observation by Jeltsch and Pohl [SIAM J. Sci. Comput., 16 (1995), pp. 40--49] that the convergence rate of a multisplitting algorithm depends on the overlap.Numerical results are presented which support the convergence theory.
Introduction
. The basic ideas of waveform relaxation were introduced in
the late 19th century by Picard [18] and Lindel-of [11] to study initial value problems.
There has been much recent interest in waveform relaxation as a practical parallel
method for the solution of stiff ordinary differential equations (ODEs) after the publication
of a paper by Lelarasmee and coworkers [10] in the area of circuit simulation.
There are two classical convergence results for waveform relaxation algorithms
for ODEs: (i) for linear systems of ODEs on unbounded time intervals one can show
linear convergence of the algorithm under some dissipation assumptions on the splitting
([15], [14], [4] and [9]); (ii) for nonlinear systems of ODEs (including linear ones)
on bounded time intervals one can show superlinear convergence assuming a Lipschitz
condition on the splitting function ([15], [1] and [3]). For classical relaxation
methods (Jacobi, Gauss Seidel, SOR) the above convergence results depend on the
discretization parameter if the ODE arises from a PDE which is discretized in space.
The convergence rates deteriorate as one refines the mesh.
Jeltsch and Pohl propose in [9] a multi-splitting algorithm with overlap, generalizing
the eliptic analysis of O'Leary and White in [17] to the parabolic case. They
prove results (i) and (ii) for their algorithm, but the convergence rates are mesh de-
pendent. However they show numerically that increasing the overlap accelerates the
convergence of the waveform relaxation algorithm. We quantify their numerical results
by formulating the waveform relaxation algorithm at the space-time continuous
level using overlapping domain decomposition; this approach was motivated by the
work of Bj-rhus [2]. We show linear convergence of this algorithm on unbounded
time intervals at a rate depending on the size of the overlap. This is an extension of
the first classical convergence result (i) for waveform relaxation from ODEs to PDEs.
Discretizing the algorithm, the size of the physical overlap corresponds to the overlap
of the multi-splitting algorithm analyzed by Jeltsch and Pohl. We show furthermore
that the convergence rate is robust with respect to mesh refinement, provided the
physical overlap is hold constant during the refinement process.
Giladi and Keller [8] study superlinear convergence of domain decomposition algorithms
for the convection diffusion equation on bounded time intervals, hence generalizing
the second classical waveform relaxation result (ii) from ODEs to PDEs.
It is interesting to note that, using multigrid to formulate a waveform relaxation
algorithm, Lubich and Osterman [13] prove linear convergence for the heat equation
independent of the mesh parameter.
In section 2 we consider a decomposition of the domain into two subdomains. This
section is mainly for illustrative purposes, since the analysis can be performed in great
detail. In section 3 we generalize the results to an arbitrary number of subdomains.
In section 4 we show numerical experiments which confirm the convergence results.
Although the analysis presented is restricted to the one dimensional heat equation,
the techiques applied in the proofs are more general. Future work successfully applies
these techniques to higher dimensional problems and to nonlinear parabolic equations.
2. Two Subdomains.
2.1. Continuous Case. Consider the one dimensional heat equation on the
interval [0; L],
where we assume f(x; t) to be bounded on the domain [0; L] \Theta [0; 1) and uniformly
H-older continuous on each compact subset of the domain. We assume furthermore
that the initial data u 0 (x) and the boundary data g 1 (t), g 2 (t) are piecewise continuous.
Then (2.1) has a unique bounded solution [5]. We consider in the following functions
in with the infinity norm
The maximum principle, and a corollary thereof, establishing the steady state solution
as a bound on the solution of the heat equation are instrumental in our analysis.
Theorem 2.1. (Maximum Principle) The solution u(x; t) of the heat equation
attains its maximum and minimum either on the initial line
or on the boundary at attains its maximum in the
interior, then u(x; t) must be constant.
Proof. The proof can be found in [21].
Corollary 2.2. The solution u(x; t) of the heat equation (2.1) with f(x;
and satisfies the inequality
Proof. Consider ~
solving
~
The solution ~
u of (2.3) does not depend on t and is given by the steady state solution
~
By construction we have ~
and on the boundary
in the kernel of the heat operator, we have by the maximum
principle for the heat equation ~ on the whole domain [0; L]. Hence
Likewise ~
and is in the kernel of the heat
operator. Hence
Therefore we have
Now the right hand side does not depend on t, so we can take the supremum over t,
which leads to the desired result.
To obtain a waveform relaxation algorithm, we decompose the
into two overlapping
in figure 2.1. The solution u(x; t) of (2.1) can now
x
Fig. 2.1. Decomposition into two overlapping subdomains.
be obtained from the solutions v(x; t)
which satisfy the
equations
@t
and
@w
@t
First note that
are solutions to (2.4) and (2.5).
Uniqueness follows from our analysis of a Schwarz type iteration introduced for eliptic
problems in [19] and further studied in [12] and [6]. We get
@t
and
@t
consider the
error equations
@d k+1
@t
and
@e k+1
The following Lemma establishes convergence of the Schwarz iteration on the interfaces
of the subdomains in L 1 . Using the maximum principle convergence in the
interior follows.
Lemma 2.3. On the interfaces the error of the Schwarz
iteration decays at the rate
Proof. By Corollary 2.2 we have
jjd k+2 (x; \Delta)jj 1 - x
fiL
and
Evaluating (2.11) at and (2.10) at combining the two we obtain
inequality (2.8). Inequality (2.9) is obtained similarly.
For any function g(\Delta; t) in L 1 ([a; b]; L 1 ) we introduce the norm
a-x-b
Theorem 2.4. The Schwarz iteration for the heat equation with two subdomains
converges in L 1 ([a; b]; L 1 ) at the linear rate
jjd 2k+1 (\Delta; \Delta)jj 1;1 -
jje 2k+1 (\Delta; \Delta)jj 1;1 -
Proof. Since the errors d k and e k are both in the kernel of the heat operator
they obtain, by the maximum principle, their maximum value on the boundary or on
the initial line. On the initial line and the exterior boundary both d k and e k vanish.
Hence
Using Lemma 2.3 the result follows.
2.2. Semi-Discrete Case. Consider the heat equation continuous in time, but
discretized in space using a centered second order finite difference scheme on a grid
with n grid points and
n+1 . This gives the linear system of ODEs
where the n \Theta n matrix A (n) , the vector valued function f (t) and the initial condition
are given by
. 1
We note the following property of A (n) for later use: let p := (p
A (n)
A (n)
We denote the i-th component of a vector valued function v(t) by v(i; t) and v(t) -
u(t) is understood component wise. We establish now the discrete analogs of the
Maximum Principle and Corollary 2.2:
Theorem 2.5. (Semi-Discrete Maximum Principle) Assume u(t) solves
the semi-discrete heat equation (2.14) with
are non-negative for t - 0 and u(i;
Proof. We follow Varga's proof in [20]. By Duhamel's principle the solution u(t)
is given by
Z te A (n) (t\Gammas) f(s)ds
The key is to note that the matrix e A (n) t contains only non-negative entries. To see
why write A contains only non-negative entries and I (n)
is the identity matrix of size n \Theta n. We get
e A (n)
J l
where the last expression has clearly only non-negative entries. Since the matrix
exponential in (2.18) is applied only to vectors with non-negative entries, it follows
that u(t) can not become negative.
Corollary 2.6. The solution u(t) of the semi-discrete heat equation (2.14) with
satisfies the inequality
Proof. Consider ~
solving
@~ u
@t
~
with ~
. Using the properties (2.16) and
of A (n) and the linearity of (2.20) we find that the solution ~
u of (2.20) does
not depend on t and is given by the steady state solution
~
The difference OE(j; t) := ~ satisfies the equation
and hence by the discrete maximum principle OE(j;
Thus
Likewise from /(j; t) := ~
Hence we can bound the modulus of u by
Now the right hand side does not depend on t, so we can take the supremum over t,
which leads to the desired result.
We decompose the domain into two overlapping
and\Omega 2 as in
figure 2.2. We assume for simplicity that ffL falls on the grid point a and fiL on
bc
bc
bc
x
Fig. 2.2. Decomposition in the semi-discrete case.
the grid point fiL. For notational
convenience we define
y
z
As in the continuous case, the solution u(t) of (2.14) can be obtained from the solutions
which satisfy the equations
@t
and
@w
Applying the Schwarz iteration to (2.21) and (2.22) we obtain
@t
consider the error equations
@d k+1
d
and
@e k+1
@t
e
The following Lemma establishes convergence of the Schwarz iteration on the interface
nodes of the subdomains in L 1 . Using the discrete maximum principle convergence
in the interior then follows.
Lemma 2.7. On the interface gridpoints a and b the error of the Schwarz iteration
decays at the rate
Proof. By Corollary 2.6 we have
jje
and
Evaluating (2.28) at a and (2.27) at a and combining the two we get
Now using we get the desired result. The
second inequality (2.26) is obtained similarly.
For any vector valued function h(t) in L 1
1!j!n
Theorem 2.8. The Schwarz iteration for the semi-discrete heat equation with
two subdomains converges in L 1 at the linear rate
jjd 2k+1 (\Delta; \Delta)jj 1;1 -
jje 2k+1 (\Delta; \Delta)jj 1;1 -
Proof. By Corollary 2.6 we have
Using Lemma 2.7 the result follows.
3. Arbitrary number of subdomains. We generalize the two subdomain case
described in section 2 to an arbitrary number of subdomains N . This leads to an
algorithm which can be run in parallel. Subdomains with even indices depend only
on subdomains with odd indices. Hence one can solve on all the even subdomains
in parallel in one sweep, and then on all the odd ones in the next one. Boundary
information is propagated in between sweeps.
Consider N
subdomains\Omega i of
so that all the subdomains overlap,
as in figure 3.1. We assume also that fi i - ff i+2 for so that domains
which are not adjacent do not overlap. The solution u(x; t) of (2.1) can be obtained
x
Fig. 3.1. Decomposition into N overlapping subdomains.
as in the case of two subdomains by composing the solutions v i (x;
which satisfy the equations
where we have introduced for convenience of notation the two functions v 0 and v N+1
which are constant in x and satisfy the given boundary conditions, namely v 0 (x;
(t). The system of equations (3.1), which is coupled through
the boundary, can be solved using the Schwarz iteration. We get for
@t
where again v k
and consider the error equations (compare figure 3.2)
e i+2
e i\Gamma2
e
ff i\Gamma2
ff
Fig. 3.2. Overlapping subdomains and corresponding error functions e i
@e k+1
with e k
For the following Lemma, we need some additional definitions to facilitate the
notation. We define ff and the constant functions
Lemma 3.1. The error e k+2
i of the i-th subdomain of the Schwarz iteration (3.3)
decays on the interfaces at the rate
jje k+2
jje k+2
the ratios of the overlaps are given by
Proof. By Corollary 2.2 we have
jje k+2
Since this result holds on all the
subdomains\Omega i , we can recursively apply it to the
errors on the right in (3.7), namely
Substituting these equations back into the right hand side of (3.7) and evaluating (3.7)
at leads to inequality (3.4). Evaluation at leads to inequality
(3.
This result is different from the result in the two subdomain case (Lemma 2.3),
because we cannot get the error directly as a function of the error at the same location
two steps before. The error at a given location depends on the errors at different
locations also. This leads to the two independent linear systems of inequalities,
where the inequality sign here means less than or equal for each component of the
vectors - k+2 and j k+2 . These vectors and the matrices D and E are slightly different
if the number of subdomains N is even or odd. We assume in the sequel that N is
even. The case where N is odd can be treated in a similar way. For N even we have
and the banded
and
pN
Note that the infinity norm of D and E equals one. This can be seen for example for
D by looking at the row sum of interior rows,
The boundary rows however sum up to a value less than one, namely
A similar result holds for the matrix E. Since the infinity norm of both D and E equals
one, convergence is not obvious at first glance. In the special case with two subdomains
treated in section 2 the matrices E and D degenerate to the scalar q 1 which is
strictly less than one and convergence follows. In the case of N subdomains the
information from the boundary needs to propagate inward to the interior subdomains,
before the algorithm exhibits convergence. Hence we expect that the infinity norm
of E and D raised to a certain power becomes strictly less than one. We need the
following Lemmas to prove convergence.
Lemma 3.2. Let r(A) 2 IR p denote the vector containing the row sums of the
p \Theta p square matrix A. Then
Proof. Let . Then we have
A n+1
Lemma 3.3. Let A be a real p \Theta q matrix with a ij - 0 and B be a real q \Theta r matrix
the sets I
for the product C := AB we have
I
Proof. We have, since a ik ; b kj - 0
a
Hence for fixed only if I
Lemma 3.4. D k and E k have strictly positive entries for all integer k - N \Gamma1Proof. We show the proof for the matrix D, the proof for E is similar. The row
index sets I i (D) are given by
I
even
oe
even
oe
even
oe
The column index sets are given by
even
oe
We are interested in the growth of the index sets I i (D k ) as a function of k. Once
every index set contains all the numbers 1 the matrix D k has strictly
positive entries. We show that every multiplication with D enlarges the index sets
I i (D k ) on both sides by two elements, as long as the elements 1 and are not
yet reached. The proof is done by induction: For D 2 we have using Lemma 3.3
I
even
oe
even
oe
even
oe
Now suppose that for k we obtained the sets
I
even
oe
even
oe
even
oe
Then for applying Lemma 3.3 again
I i (D k+1 )=
even
oe
even
oe
even
oe
Hence every row index set I i (D k ) grows on both sides by 2 when D k is multiplied by
D, as long as the boundary numbers 1 and are not yet reached. Now the index
set I 1 (D k has to grow most to reach the boundary number
we need for the number of iterations
1for the matrix D k to have strictly positive entries.
The infinity norm of a vector v in IR n and a matrix A in IR n\Thetan is defined by
1-j-n
1-i-n
Lemma 3.5. For all k ? N
there exists
Proof. We prove the result for D; the proof for E is similar. We have from (3.11)
and (3.12) that
By Lemma 3.4 D k has strictly positive entries for any k - N
2 . Note also that jjD k jj 1 -
1 since jjDjj 1 - 1. Now by Lemma 3.2 we have
since D k
1 and r N (D) ! 1.
Remark: It suffices for each row index set to reach one of the boundaries, either
1 or N \Gamma 1, for the infinity norm to start decaying. Hence it is enough that there are
no more index sets I i (D k ) (compare the proof of Lemma 3.4) such that 2
so that the requirement k - N \Gamma1
2 can be relaxed to k ? N \Gamma3
4 .
We now fix some k ? N \Gamma3and set
Lemma 3.6. The vectors - and j satisfy
Proof. By induction on (3.8), using that the entries of D, E, - k and j k are
non-negative, we get
Taking norms on both sides and applying Lemma 3.5 the result follows.
Theorem 3.7. The Schwarz iteration for the heat equation with N subdomains
converges in the infinity norm in time and space. We have
jje 2km+1
jje 2km+1
Proof. We use again the maximum principle. Since the error e k
i is in the kernel of
the heat operator, by the maximum principle e k
attains its maximum on the initial
line or on the boundary. On the initial line e k
jje 2km+1
jje 2km+1
Using Lemma 3.6 the result follows.
Note that the bound for the rate of convergence in Theorem 3.7 is not explicit.
This is unavoidable for the level of generality employed. But, if we assume for simplicity
that the overlaps are all of the same size then we can get more explicit rates
of convergence. We set r
The matrices D and E then simplify to
~
pr
pr
pr
pr
pr
pr
pr
In this case we can bound the spectral norm of ~
D and ~
E by an explicit expression
less than one. We use common notation for the spectral norm, namely
Lemma 3.8. The spectral norms of ~
D and ~
are bounded by
Proof. We prove the bound for ~
D. The bound for ~
can be obtained similarly.
We can estimate the spectral norm of ~
D by letting ~
and F has only O(N) nonzero entries and these are equal to 1. In fact jjF jj
Using that the eigenvalues of J are given by
the spectral norm of ~
D can be estimated by
since
Lemma 3.9. Assume that all the N subdomains overlap at the same ratio r 2
(0; 0:5]. Then the vectors - and j satisfy
Proof. The proof follows as in Lemma 3.6.
Note that r = 0:5, which minimizes the upper bound in Lemma 3.9, corresponds
to the maximum possible overlap in this setting, namely fi Figure 3.2.
Theorem 3.10. The Schwarz iteration for the heat equation with N subdomains
that overlap at the same ratio r 2 (0; 0:5] converges in the infinity norm in time and
space. Specifically we have
jje 2k
jje 2k
Proof. From the proof of Theorem 3.7 we have
jje 2k+1
jje 2k+1
Since the infinity norm is bounded by the spectral norm we get
jje 2k+1
jje 2k+1
Using Lemma 3.9 the result follows.
The results derived above for the continuous heat equation remain valid as in the
two subdomain case, when the heat equation is discretized. Details of this analysis
can be found in [7].
4. Numerical Experiments. We perform numerical experiments to measure
the actual convergence rate of the algorithm for the example problem
To solve the semi-discrete heat equation, we use the Backward Euler method in time.
The first experiment is done splitting the
3] into the two
3] for three pairs of values
0:52)g. Figure 4.1 shows the convergence of
the algorithm on the grid point b for 0:01. The solid line is
the predicted convergence rate according to Theorem 2.8 and the dashed line is the
measured one. The measured error displayed is the difference between the numerical
solution on the whole domain and the solution obtained from the domain decomposition
algorithm. As initial guess for the iteration we used the initial condition constant
in time. We also checked the robustness of the method by refining the time step and
obtained similar results.
We solved the same problem (4.1) using eight subdomains which overlap by 35%.

Figure

4.2 shows the decay of the infinity norm of - k . The dashed line shows the
measured decay rate and the solid line the predicted one. Note that in the initial
phase of the iteration the error stagnates, since information has to be propagated
across domains.

Acknowledgments

We thank Gene Golub for showing us how to prove Lemma
3.8 and Olavi Nevanlinna, Morten Bj-rhus and Sigitas Keras for many interesting
discussions.
Iteration k
at
grid
point
Fig. 4.1. Theoretical and measured decay rate of the error for two subdomains and three
different sizes of the overlap
Iteration k
Infinity
norm
of
xi^k
Fig. 4.2. Theoretical and measured decay rate of the error in the case of eight subdomains



--R

The Use of Runge-Kutta Formulae in

A Note on the Convergence of Discretized Dynamic Iteration
Parallel and Sequential Methods for Ordinary Differential Equations
The One-Dimensional Heat Equation
Domain Decomposition Algorithms
Analysis of Parallel Algorithms for Time Dependent Partial Differential Equa- tions

Waveform relaxation with overlapping splittings
The
Sur l'application des m'ethodes d'approximations successives 'a l"etude des int'egrales r'eeles des 'equations diff'erentielles ordinaires.
On the Schwarz Alternating Method

Convergence of Dynamic Iteration Methods for Initial Value Problems
Remarks on Picard-Lindel-of Iterations
Domain Decomposition and Iterations in Parabolic Problems

Sur l'application des m'ethodes d'approximations successives 'a l"etude de certaines 'equations diff'erentielles ordinaires

Matrix Iterative Analysis
Partial Differential Equations of Applied Mathematics
--TR

--CTR
D. Firsov , S. H. Lui, Domain decomposition methods in image denoising using Gaussian curvature, Journal of Computational and Applied Mathematics, v.193 n.2, p.460-473, 1 September 2006
S.-H. Lui, On monotone iteration and Schwarz methods for nonlinear parabolic PDEs, Journal of Computational and Applied Mathematics, v.161 n.2, p.449-468, 15 December
overlapping Schwarz;multisplitting;domain decomposition;waveform relaxation