Timestep Acceleration of Waveform Relaxation.
Dynamic iteration methods for treating certain classes of linear systems of differential equations are considered. It is shown that the discretized Picard--Lindelf (waveform relaxation) iteration can be accelerated by solving the defect equations with a larger timestep, or by using a recursive procedure based on a succession of increasing timesteps. A discussion of convergence is presented, including analysis of a discrete smoothing property maintained by symmetric multistep methods applied to linear wave equations. Numerical experiments indicate that the  method can speed convergence.
Introduction
. Much of modern chemical and physical research relies on the numerical
solution of various wave equations. Since these problems are extremely demanding of both
storage and cpu-time, new numerical methods and fast algorithms are needed to make optimal
use of advanced computers. The dynamic iteration or waveform relaxation (WR) method [9,
11] is an iterative decoupling scheme for ordinary differential equations which can facilitate
concurrent processing of large ODE systems for applications such as VLSI circuit simulation
[6, 15] and partial differential equations [1, 4].
In this article, accelerated dynamic iteration schemes are used to solve systems of linear
differential equations, with emphasis on the ordinary differential equations arising from discretization
of linear wave equations. Although our experiments use finite differences for the
spatial derivatives, other spatial discretizations could be used. For time discretization, we use
symmetric multistep methods, although other choices may also be appropriate. As is the case
for stationary iterative methods applied to spatially discretized elliptic PDEs, it is found that
finer fixed step (time) discretizations slow the convergence of the WR iteration, while large
timesteps can be used to resolve the slow modes. The idea that is explored here is to use
a coarse timestep on the defect equations to speed up convergence of the fine grid iteration.
Nevanlinna already pointed out [14, 13] that for general applications of WR it makes sense from
Department of Mathematics, University of Kansas, Lawrence, KS 66045, leimkuhl@math.ukans.edu. This
work was supported by NSF grant DMS-9303223
an efficiency standpoint to use coarser discretization in the early sweeps (when the iteration
error is large), and then incrementally to refine the time discretization near convergence. Our
point of view is rather to vary the timestep to resolve different modes present in the solution,
using two time stepsizes (or multiple stepsizes). The current article is related to recent work of
Horton and Vandewalle [4] and Horton, Vandewalle and Worley [5] which considered space-time
multigrid methods for solving parabolic equations.
The new scheme will be referred to as timestep acceleration since it relies on adjustment of
the integration timestep to accelerate the dynamic iteration. This approach shares some features
of multigrid methods. For the convenience of the reader generally familiar with multigrid
methods, we outline the algorithm in the abstract setting of solving an unspecified dynamical
system as follows:
Accelerated Waveform Relaxation. Given: fine timestep h and an approximation u 0 to
the solution with fixed stepsize h.
1. Smoothing Starting from u 0 , perform a fixed number of iterations of a smoothing waveform
relaxation iteration with timestep h.
2. Correction Compute the defect (residual) in this solution on the fine time mesh.
If the timestep sufficiently large, solve the discrete defect equation restricted
to the coarse time mesh directly (i.e. without relaxation).
Else recursively apply some number of iterations of the algorithm using stepsize H to
the defect equation restricted to the coarse time mesh.
Next, correct the solution after prolonging onto the fine time mesh.
3. Smoothing Apply a fixed number of iterations of the fine stepsize smoothing iteration.
(In x6, we present and analyze a more precisely defined version of this algorithm, TAWR.)
A major barrier to efficient solution of large scale wave equations is the need for small
timesteps. Due to the sequential character of standard ODE methods, this effectively reduces
the potential for parallel speedups. Compared to standard timestepping schemes, the method
discussed here directly addresses this problem by enabling the use of larger timesteps to recover
at least a portion of the dynamics. Another important obstacle to computation-particularly
in the case of high dimensional problems-is the necessary storage. The new method actually
exacerbates this problem since solution information at many points must be stored. However, in
waveform relaxation based on a block splitting, the storage is naturally segmented according to
the decoupling, so the scheme may be appropriate for a parallel computer based on a distributed
memory architecture.
Although standard analytical results for multigrid methods or coarse-grid acceleration are
typically developed for finite dimensional Hermitian positive definite problems, these can be
relaxed to give at least partial convergence results. In fact proving theoretical convergence
for timestep acceleration is easier than for standard multigrid due to the strong smoothing
properties of the Picard-Lindel-of operator (it is a contraction on small intervals). Analysis of
the behavior of the iteration on special linear model problems is also possible and is briefly
discussed here.
The scheme is found to work well in simple numerical experiments with linear wave equa-
tions. Although our experiments are conducted in one space dimension, nothing in principle
prevents application in higher dimensions (although many practical issues will need to be dealt
with).
2. Waveform Relaxation. Consider a second order linear system of differential equation

(1)
where the eigenvalues of the matrix A are assumed to lie in the left half plane. A special
case that we will frequently refer to is the 1-D wave equation U discretized with finite
differences on the unit square with periodic boundary conditions:
. 0
A is called the discrete Laplacian. Here is a vector of approximations
at nodes x i ,
Another potential application is to the Schr-odinger equation. Discretizing, for example,
with finite differences leads to
d
dt
(2)
where A is the discrete Laplacian. If v(x; t) is the potential energy function of the corresponding
classical system, we have V In simplified settings
v(x; t) is time-independent, hence so is V .
The waveform relaxation method for (1) is based on a splitting results
in ODE IVPs
For example, we might choose A+ to be the diagonal of A (Jacobi splitting), a block-diagonal
part of A (block-Jacobi splitting), the lower triangular part of A (Gauss-Seidel splitting), etc.
Much work involving the discrete Laplacian in elliptic PDEs is based on Gauss-Seidel splitting
in red-black ordering. Another useful splitting is the damped Jacobi splitting where
with D the diagonal of A. The extreme case is called the Picard splitting.
When referring to (1) and (3) we will generally limit discussion to the case where A and
are symmetric negative semidefinite matrices.
The WR iteration proceeds as follows: starting from a given initial waveform u
(which may be constant), we solve (3) with as a forced linear system for u 1 over some
time interval, say [0; T ]. (This interval is referred to as the window). The function u 1 then
yields a forcing for the next iteration or sweep, and the process repeats. In practice, the systems
are solved numerically over the entire interval, and the storage of the resulting discrete
approximation is an important drawback of the method which may place severe limitations
on the size of the time window. On the other hand, we gain in two ways: first, the systems
we solve at each iteration can be decoupled into problems of reduced dimension, and second,
the decoupled problems can often be solved on separate processors of a parallel computer. An
alternative approach would be based on solving the linear equations that result at each step
of a standard discretization using a parallel algorithm, however, depending on the computer
architecture employed, the flexibility in the choice of window size may reduce the overall communication
cost, e.g. by eliminating some of the time spent in initializing the transfer of data
between processors.
Preliminary convergence results for WR appear in the paper by Lelarasmee et al [9].
Miekkala and Nevanlinna [11, 12] and Nevanlinna [13] have developed an extensive theory for
studying waveform relaxation for linear systems. Lubich and Osterman proposed to combine
the WR method with spatial multigrid schemes [10]. Recent work by Horton and Vandewalle
[4] and by Horton, Vandewalle and Worley [5] has shown that a careful implementation of (spa-
tial) multigrid-WR methods for parabolic PDEs can provide excellent parallel speedups. The
use of waveform relaxation for solving hyperbolic partial differential equations and relations to
domain decomposition were explored by Bj-rhus [1].
3. Mathematical Background. In this section we state some elementary results concerning
the iteration (3). The reader is directed to the papers of Nevanlinna and Miekkala for
basic theory.
The waveform relaxation method for (1) can be viewed as an iteration u
As shown in [11], implying superlinear convergence. On the other hand,
for stiff dissipative linear systems, it makes sense to allow T ! 1 in which case meaningful
spectral information is obtained [11]. Since the solution to the equations (1) and (2) does not
in general lie in L 2 ([0; 1)), this approach must be modified. A reasonable practical approach
is that taken in [13] where an exponential weighting function e \Gammafft is inserted into the usual L 2
norm. For ff ? 0, the space L 2
ff is normed by
kuk ff :=
-Z 1je \Gammafft u(t)j 2 dt
and ae ff refers to spectral radius in that space.
If we take the Laplace transform of (3), we obtain:
u(0).
The following results are proved by Nevanlinna and Miekkala [11]:
ae ff
Rez-ff
which follows from the Paley-Wiener theorem (the second expression follows from a maximum
principle after a suitable remapping of the domain) and
which follows from Parseval's identity.
We now provide some simple estimates for the response of the iteration operator in weighted
2-norm.
First, consider the behavior of the solution operator of (1) in the
weighted space, where d=dt. Examining the spectral radius of the normal matrix L(z)
along the line Rez = ff, we find the eigenvalues are:
are the eigenvalues of A. Hence
By maximizing these functions over y, we can compute the moduli of the eigenvalues of solution
operator in the weighted space.
Theorem 3.1. Define
Then
ae ff (L
In particular eigenvalues near zero have the strongest influence. When A is the discrete
Laplacian, or any symmetric negative semidefinite matrix with an eigenvalue at zero, we have
ae
We can use this to estimate the norm of the iteration matrix, since
where the -
are determined from Theorem 3.1 with - i the eigenvalues of A+ rather than A.
Asymptotic estimates for the relation between spectral radius of S(z) and
ff are given in [8].
Let us consider the wave equations with periodic boundary conditions on the square as a
model problem. We will use damped Jacobi iteration with parameter ! 2 (0; 1]. In this case,
A, are all diagonalized by the discrete Fourier transform, so we arrive readily at
the eigenvalues
e
and
In this case the spectral radius can be readily computed. We have
\Deltax
pp
pp
and
ae ff
We are interested in moderate weights ff which we define to mean ff !
ae(A). (Intuitively, this
corresponds to looking in the time domain on intervals greater than the smallest period of the
motion.)
In the standard theory, one uses the value of ! in the damped Jacobi splitting to enhance a
smoothing property: a damping in the iteration of the modes corresponding to larger eigenvalues.
However, the important consideration for timestep acceleration is not the way in which the
smoother acts on fast "spatial" modes but rather the response of the smoother to high frequency
forcings. In fact, the real smoothing property we are interested in has to do with the shape of
the graph of ae(S(ff + iy)) as a function of y. For example, when a damped Jacobi splitting is
applied to solve the semidiscrete wave equation, we find that the spectral radius of S achieves
its maximum on Rez = ff at the point (if ff !
(or at
). The maximum is typically achieved well away from
For the Picard splitting, it is easy to see rather that the maximum occurs at In
this case, we say that the iteration has a smoothing property with respect to high frequency
forcings. It is not necessary to use a slowly converging splitting such as the Picard splitting
to obtain a good smoothing property. A typical feature of a good splitting for this purpose
is that A+ would have an eigenvalue at or near the origin. Thus the smoothing property of a
block-Jacobi splitting of the discrete Laplacian would improve with the block size.
To illustrate this smoothing concept, consider the time-dependent Schr-odinger equation
(2). The iteration becomes
d
dt
We could again use a Jacobi or damped Jacobi splitting, but in practice, a more useful choice
might be
or, more simply
After time discretization, these choices will lead to equations at each timestep which can be
efficiently solved by, for example, using a parallel implementation of the fast Fourier transform
Still another possibility is to work directly in the Fourier coefficients. Let QAQ
where
and
so that the equations become
d
dt
\Psi:
We can then apply a Jacobi splitting to this problem. One finds that the diagonal of QV (t)Q H
is dI where
The computation can be implemented efficiently using the FFT.
The symbol of the WR iteration operator R for the Schr-odinger equation with V
constant is
and, trivially,
ff
Using the Laplacian+potential splitting, or one of its cousins can be expected to yield a
good smoothing property with respect to high frequency forcings.
4. Discretization. In this section, we focus on (1) and apply a discrete transform as in
[12] to analyze the symmetric multistep methods commonly used for integrating oscillatory
problems.
Multistep methods construct an approximating sequence fu n g to fu(t n )g at successive time
points use fu k
n g to refer to the numerical solution generated at the kth sweep of
waveform relaxation. Symmetric multistep methods for -
take the form:
with sequences, for which ff
These methods are used for integration of second order oscillatory
problems. An important feature of this class of methods is their time-reversibility. Note that
the multistep methods require k starting values.
In discretizing dissipative problems it is sensible to replace the space L 2 by l 2
h with norm
. For our investigations, we use the weighted space with norm
je \Gammanhff u
which can be viewed as a discretization of the L 2
ff norm.
Following the usual practice we define operators a and b on sequences by
We also use the symbols a and b to refer to the corresponding characteristic polynomials:
To preserve the intuitive correspondence of results from the continuous-time to discrete
worlds, define a discrete transform which takes fu n g to - u(z) by -
tially a discretization of the Laplace transform, and equivalent to the discrete Laplace or i
transform).
Applying (5) to the linear problem (3) and computing the discrete transform, we find
where
and OE includes the effects due to the k starting values. We are going to assume that these starting
values are exact (for the unsplit discrete problem) so they do not effect the convergence of the
iteration.
For example, if
where
In the general case, discrete versions of the Paley-Wiener theorem and Parseval's identity
give (after modifying results in [12] to take into account the exponential weight):
ae h;ff (S h
Rez-ff
ae(S h (z))
and
In order for the discretized operator to be bounded, we evidently need to require
for any - 2
We now consider an example. Ignoring rounding error, the popular leapfrog method for
second order systems is equivalent to St-ormer's rule (also known as the Verlet method), a
symmetric 2-step method with ff 1. Applying this
scheme to the WR iteration for the linear problem and taking the discrete transform gives
The function
is an an O(h 2 z 4 ) approximation to z 2 . The poles of the transformed discrete iteration operator
z =h
cosh
with - an eigenvalue of A+ . Explicit multistep scheme are always conditionally stable meaning
that the stability of the schemes will depend on the stepsize being restricted roughly in inverse
relation to the square root of the spectral radius of A+ . For the St-ormer method, the stability
condition is that - ! 0 and \Gammah 2 - != 2, which is also the condition that the poles of S h remain
on the imaginary axis. The function Im cosh(-) is monotone in the real variable - on [\Gamma1; 1],
hence the ordering of the poles is preserved along the imaginary axis.
Another popular second order method is the (implicit) trapezoidal rule which has transform
4.1. Decay of the discrete symbol. Theory due to Miekkala and Nevanlinna [12] compares
the convergence of the discrete iteration in l 2
h to that of the continuous time iteration
for dissipative problems and for methods that are not weakly stable. We need to modify this
mechanism to cover convergence for stable methods for second order differential equations in
the weighted spaces rather than L 2 and l 2
h . In what follows, it is assumed that the k starting
values are held fixed as we iterate. These could also be obtained by some convergent process,
but this does not seem a meaningful generalization.
In the case of the discretized iteration we need to examine the images of segments ff
. The situation for representative and in Figure 1 we see
the images of this line for the trapezoidal and St-ormer discretizations, for various values of the
stepsize.
Putting real ff into each of the functions P t:r:
h and z 2 one can show that for sufficiently
small hff
which means, somewhat surprisingly, that in the neighborhood of y = 0, the St-ormer discretization
actually leads to a slightly more stable overall iteration than that generated by the
trapezoidal rule.
h=.2
h=.3
h=.4
h=.1 h=.2
Fig. 1. Approximation of z 2 by St-ormer (dashed lines) and trapezoidal rule (dotted lines)
The situation for large h is more dramatic. For nonstiff problems with eigenvalues - very
near the origin in the complex plane, large steps should be possible and one might suppose
that the St-ormer and trapezoidal rule discretizations would behave similarly with respect to
WR convergence. In fact, this is not the case and it turns out that the St-ormer method yields
a much more stable WR iteration than the trapezoidal rule over comparable time intervals.

Figure

2 shows the image of 1
h and z 2 .

Figure
also indicates that results such as Proposition 9 of [10] and Theorem 3.4 of [12] which bound the
spectral radius of the discrete iteration in terms of that of the continuous iteration for A-stable
multistep methods will not typically hold in our setting.
The problem with generalizing the results of [12] is that they were based on the strengthened
stability assumption that the stability region includes a disk on the negative real axis touching
the origin, whereas many of the symmetric methods we consider (e.g. St-ormer's rule) do not
satisfy this condition.
We will use the exponential weight to correct for the weak stability of the method. For
define the fl-stability
region\Omega fl of the method as the set of all - 2 -
C such that the roots
lie in the disk jij - e fl and are simple on the boundary. The iteration
operator S h is bounded in l 2
h;ff if
2.5 -2 -1.5
-113s.r. h=2
t.r. h=2
Fig. 2. Large time-step comparison of images of 1
Now observe that
We can directly relate the spectral radius of the discrete iteration to that of the continuous
time iteration. In fact,
ae ff;h (S h
Rez-ff
ae
Rez-ff
ae
s
Let the notation bdyW be used to indicate the boundary of the set W . Since f
r
bdy\Omega hff , we have, analogous to a result in [12].
Theorem 4.1. Suppose
int\Omega ffh , then
ae ff;h (S h
int\Omega ffh
and
@\Omega ffh g:
Theorem 4.2. If the dynamic iteration converges in L 2
ff and the symmetric multistep
method is irreducible and convergent then the discretized iteration converges in l 2
ff;h for sufficiently
small h and
ae ff;h (S h
We will outline a proof of this result since the reasoning is somewhat different than that
used in [12].
be the k zeros of a, with being the principle root counted with
multiplicity two. For simplicity, assume that these zeros all lie on the unit circle S 1 and that
they are ordered counterclockwise about the unit circle, thus fi
(It would not be difficult to treat the case where some zeros lie inside the unit circle.) From
consistency, we must have that fi double root, while all of the other roots are simple.
We can view ja(e hff w)j 2 as a function of w on S 1 . For at the
zeros of a; for h sufficiently small and ff ? 0, it has local minima located near the points
. We can expand a in Taylor's series about the fi i to obtain
Only multiple root of a, hence a 0 1. This means that, for
e hff w in the vicinity of fi i , must have a(e hff must
also hold at the local minimum.
Using this, we can prove a small lemma which shows that the spectral radius is determined
for small h by the approximation property of the principle root.
Lemma 4.3. For h sufficiently small,
ae(S
ae(S
and a similar result holds for kS h k.
Proof: By symmetry, ' sufficiently small,
the global minimum of ja(e hff e ihy )j 2 on I h must occur at one of its local minima over that
interval or at the endpoints. Since b(e hz ) can be uniformly bounded in any bounded region, it
is straightforward to see that the quantity - h (y) defined by
s
satisfies
min
and hence that
ae(S
Due to symmetry, the behavior of ae is the same on the intervals ['
other words, we need only look in the latter subinterval for the maximizing value. 2.
Given -
consistency implies that
0-y-y
Choose a large enough rectangular neighborhood N of the origin so that e.g. ae(S(z)) ! ae ff (S)=2
for z outside N . Now j- h (' 2 =h)j - h \Gamma1=2 , thus for h sufficiently small, the curve \Gamma h := f- h (y) :
leaves N . After leaving N , it cannot reenter N (or j- h (y) 2 j would have another
local minimum). Within N , the curve \Gamma h will approximate the line segment ff + iy to O(h).
Thus for h sufficiently small, the maximum value of ae(S h (z)) will occur when - h (y) lies within
N , and since this point lies within O(h) of ff iy we can see that asymptotically, the spectral
radius of S h in the weighted space can differ by only O(h) from that of S. This concludes the
proof of the theorem.
5. Aliasing effects. Consider the transformed discrete iterator S h
on the vertical line ff R. The degree to which an eigenvalue \Gamma! 2 of A+ has an impact
on the solution at frequency - depends inversely on the separation between P h (ff + i-) and
\Gamma! 2 . Those frequencies - for which P h (ff lies far from the spectrum of A+ will be only
weakly propagated by the iteration.
For any multistep method, the function P h (ff + i-) is actually periodic in - with period
2-=h. This aliasing effect means that high frequencies can be excited with large stepsizes.
Frequencies give the same response. Actually, the situation is even somewhat
worse due to the symmetry about the real axis: the response to \Gamma- will be the same as
the response to - 0 . Of course, if there are no frequencies present in the forcing function above
say -=h then these anomalous excitations do no harm.
We will illustrate with the wave equation example. We first look at the response of the
discrete solution operator for the (unsplit) spatially discretized wave equation along
the line 1 iy. The curves shown in Figure 3 show the spectral radii (hence also the norm) of
versus y for :25. The maximum value is achieved near expected
from Theorems 3.1 and 4.2. An increase in the stepsize h provides accuracy for small y while
introducing some extraneous excitation at 2-k=h, k 2 Z.
solid h=1
-. h=.5
. h=.25
Fig. 3. Spectral radius of L(1+iy), wave equation, N=32.
We next examine the spectral radius of the transformed discrete iteration operator S(z)
for the the Jacobi splitting of the spatially discretized wave equation along the line
iy. The curves shown in Figure 4 show the spectral radii versus y for h in the progression
As h decreases, the spectral radius has increasing maxima achieved at
increasing values of y.
solid h=1
-. h=.5
. h=.25
Fig. 4. Spectral radius of S(1+iy), wave equation, N=32.
Using a large stepsize to resolve the small y response will not apparently improve the
convergence of the small step Jacobi iteration, since the large stepsize solution operator does
not even act on the high frequencies (where the spectral radius is large). The exception to this
will be in the case that h is so small that the "coarse" grid is not coarse at all (in which case,
little is gained through iteration). Moreover, unless ae(S h (z)) is small outside of an interval
about the origin of length roughly 2-=H , the artifacts introduced at the high frequencies on
the coarse grid will not be damped out.
To see a substantial improvement, our iteration operator should be designed to achieve its
maximum at or near As mentioned previously, for the wave equation, a natural (if slow)
choice is standard Picard iteration.
If we turn to the Schr-odinger equation and consider for example the splitting (4) for V
constant. In case V is not too large, we would expect here that the maximum of ae(R(z)) is
achieved near has an eigenvalue near the origin) and that substantial improvement
may be possible by exploiting a coarse time step solution.
6. Timestep Acceleration of WR. The examples of the previous section suggest that
an approach in which different time meshes are used at each sweep could be successful. The
goal, as for standard multigrid is to iterate on successively coarser grids, thus resolving those
components of the residual that are most difficult to obtain on the fine grid. We envision
ultimately combining spatial multigrid with this timestep acceleration scheme. For the formulation
and analysis of standard multigrid methods in the context of elliptic PDEs, the reader
is referred to [2].
We will now define the steps of the algorithm described in the introduction. Let b and a
represent the operators which define the discretization. We use h to represent the fine timestep,
and H to represent the coarse timestep. Normally, in solving elliptic PDEs, we use In
our case, this choice may or may not be appropriate; for the purposes of discussing an algorithm,
we assume that the stepsize changes by a common factor at each iteration, but this is perhaps
not essential in practice. Let I H
h and I h
prolongation operators,
respectively, which act between the fine and coarse time-meshes, thus I H
H;ff and
I h
h;ff .
Note that whether we wish to solve problems with or without forcing, description for a
forced problem permits an easy recursive definition.
Algorithm TAWR(h). Given: fine timestep h, a sequence ff n g 2 l 2
h;ff , an approximation
u 0 to the solution with fixed stepsize h, and a splitting , the following algorithm
solves
subject to k prescribed starting values ' i ,
1. Small-timestep Pre-Smoothing. Starting from u 0 , perform - sweeps of WR iteration
with timestep
where u l+1
2. Large-timestep Correction. Compute the defect fd n g satisfying d
and
from the formula
If the timestep sufficiently large, solve
directly (i.e. without relaxation) using zeros for starting values.
Else apply - iterations of TAWR(H) to (9), using zeros for starting values.
Next, correct:
3. Small-timestep Post-Smoothing Apply - iterations of the fine mesh smoothing operation
(8).
Notes:
ffl For this is the V-cycle, for - 2, it is called the W-cycle.
ffl Different numbers of smoothing steps could be used in the pre- and post-smoothings.
ffl To solve the problem using timestep acceleration, we first compute ff n g :=
)g.
7. Convergence Analysis. In this section we present an elementary general convergence
result regarding two-grid acceleration. This result could be easily extended to the full
timestep acceleration iteration. The iteration operator in the two stepsize case can be written
as S -
h where S h represents the smoothing sweep and C h;H represents the coarse-grid
correction. In general it is enough to understand
h . The operator C h;H can be
where we have denoted the prolongation and restriction operators by p and r, respectively. On
the other hand,
It is enough to show that
is O(h \Gamma2 ) in l 2
ff;h , while the norm of
is O(h 2 OE(-)), where OE tends to zero as - ! 0. In fact we anticipate that the situation is often
rather better than this result would indicate, but this approach allows us to state a quite general
convergence result.
Based on the relation R \Gamma1
h , the fact that b is a bounded operator, and the
theorems of the last section, we have:
Lemma 7.1. Suppose consistent, stable linear multistep is used and the restriction and prolongation
operators are bounded operators. Then kM sufficiently
small.
The proof follows since (i) kL \Gamma1
O(h), (ii) the same thing holds for h
replaced by H and qh, and (iii) the restriction and prolongation operators are bounded.
2.
For the smoothing, we have
We therefore have
This converges to zero provided
ae ff;h (S h
Thus we can state:
Theorem 7.2. If the undiscretized smoothing iteration is convergent (ae ff (S) ! 1), a consis-
tent, stable linear multistep is used, the restriction and prolongation are bounded operators between
l 2
h;ff and l 2
H;ff and enough smoothing iterations are performed, then the timestep-accelerated
waveform relaxation algorithm converges.
Because of the strong contractivity of the Picard operator on small time intervals, it would
be straightforward to extend this result to the full multiple mesh recursive acceleration scheme.
On the other hand, besides proving asymptotic convergence, this simplified approach provides
no practical estimates of convergence.
7.1. Treatment of Model Problems. A key observation is that two modes are coupled
via restriction. It is possible to write a formula for the "symbol" of the iteration operator as a
matrix operating on the pair of coupled modes e hnz and e hn(z+i-) . As an example, taking
the operators
\Theta 111
(full weighting restriction) and linear interpolation), and assuming any
symmetric multistep method (a; b), then we find the action of M on the pair of modes is given
by
where
and
(\Omega is the Kronecker product).
Now for Jacobi or Picard iteration on the wave equation, for example, the matrix -
M is
easily reduced to a diagonal matrix of 2 \Theta 2 blocks, so the asymptotic convergence behavior
can be determined relatively easily. For red-black Gauss-Seidel iteration on the square, we get
a further pairing of the spatial modes, so -
M actually is reduced to 4 \Theta 4 blocks. By studying
the spectra of these blocks, various ODE discretizations could be compared for their effect on
asymptotic rate of convergence, as could other choices of restriction/prolongation.
Note that besides the restriction and prolongation having an adjoint relationship
if we choose the second order, two-step discretization
then we find that also rR h This appears to be the only consistent, stable two-step
scheme for which this property holds with the given r and p. These resemble the conditions
for "variational form", however the operator R and its discretization are not self-adjoint in our
setting, so we do not have the space decomposition
l 2
and the standard theoretical results cannot be directly applied.
8. Numerical Experiments. We performed experiments using the two-grid iteration
on linear wave equations. We found that the performance improvements were sensitive to
many factors, including timestep, time window length, and splitting. Unfortunately, we cannot
expect to have complete flexibility in the choice of the time interval or "window" as this may
be determined from a storage or communication limitation. Similarly, the timestep is typically
chosen for accuracy reasons.
Consider the standard 1D wave equation (1), N=16, using Jacobi iteration for the smooth-
ing. We used the discretization (10) together with full weighting restriction and piecewise linear
interpolation. We did not anticipate very good behavior since the smoothing property is relatively
weak for this splitting (fast modes are not very strongly damped). Indeed, this is what
we observed. For most values of the stepsize, the 2-grid acceleration improved convergence,
but not by very much. (In some cases performance was even slightly degraded.) In each of
the Figures the 2-norm of the error is graphed as a function of the sweep number s and the
timestep n. In Figure 5 the error in Jacobi WR is indicated for stepsize
shows the mild improvement in the error when a coarse grid correction is applied at each Jacobi
WR sweep.
We next examined a modified wave equation of the form
where A is the discrete Laplacian and - is a scalar parameter. We used "Laplacian splitting"
into A and -I . It is easy to see that this splitting possesses a strong "smoothing property". We
first chose means that we had a substantial perturbation of
timestep
sweep
error
Fig. 5. Errors in Jacobi WR, h=.025, 40 steps without correction.103050515255001500timestep
sweep
error
Fig. 6. Errors with coarse grid correction, h=.025, 40 steps, showing poor acceleration. The benefit of coarse
grid correction is diminished by the poor smoothing property of the Jacobi smoother.
the discrete Laplacian. Initial data excited the first two eigenfunctions of the Laplacian (slow
modes), although this choice was not critical to the results we obtained. Twenty timesteps of
size were used. In this case, the coarse grid corrections offer substantial improvement, as
shown in Figure 7. The left figure shows the log 10 error versus timestep and sweep number; on
the right we have shown the log 10 ratio of the errors with and without the two-grid acceleration.
The improvement evidenced here is as much as a factor of 10 5 in 20 sweeps, or a little under
a factor of 2 per sweep on average.
The improvement is evident until the error reaches the level of roundoff. At the weaker
perturbation the effect somewhat diminished (Figure 8). If we instead increased the
strength of the applied field 100), the coarse grid correction continued to offer substantial
acceleration; this is indicated in Figure 9. At larger or smaller timesteps, the improvement
slightly diminished. A linear acceleration effect was observed on longer time intervals (Figure
timestep iteration sweep
log
-5
timestep iteration sweep
log
error
ratio
Fig. 7. (a) log error and (b) log error ratio, steps . The splitting provides a strong
smoothing property, and a substantial improvement is possible with the two-grid iteration.515020
timestep iteration sweep
log
error
ratio515020
timestep iteration sweep
log
error
Fig. 8. (a) log error and (b) log error ratio,
Although these experiments suggest that time-mesh coarsening accelerations hold promise
for improving the parallel waveform relaxation algorithm, they certainly do not settle all the
issues. In particular, we do not have an easy and robust mechanism for determining what
timestep iteration sweep
log
-5
-1timestep iteration sweep
log
error
ratio
Fig. 9. (a) log error and (b) log error ratio,
-5
timestep iteration sweep
log
error
ratio1030020
-55timestep iteration sweep
log
error
Fig. 10. (a) log error and (b) log error ratio, steps . On longer time intervals, the
convergence acceleration factor (ratio of errors with and without acceleration) becomes approximately linear in
the sweep.
splittings will benefit from acceleration, or for determining various parameters such as number
of smoothing sweeps, optimal coarsening, etc. We also have not yet experimented with the use
of more than two levels of time-mesh acceleration.

Acknowledgements

: The author is indebted to Pawel Szeptycki for several helpful discussions
during the early stages of this work. Stefan Vandewalle read a preliminary version of
the manuscript and contributed several very useful comments. The computers of the Kansas
Institute for Theoretical and Computational Science were used for the numerical experiments.



--R


New York
Solving Ordinary Differential Equations
A spacetime multigrid method for parabolic PDEs
An algorithm with polylog parallel complexity for solving parabolic partial differential equations
Zukowski D.

Estimating waveform relaxation convergence
The waveform relaxation method for time-domain analysis of large scale integrated circuits
Multigrid dynamic iteration for parabolic equations
Convergence of dynamic iteration methods for initial value problems
Sets of convergence and stability regions
Remarks on Picard Lindel-of iteration
Power bounded prolongations and Picard-Lindel-of iteration
Partitioning algorithms and parallel implementation of waveform relaxation algorithms for circuit simulation
--TR

--CTR
D. Guibert , D. Tromeur-Dervout, Parallel adaptive time domain decomposition for stiff systems of ODEs/DAEs, Computers and Structures, v.85 n.9, p.553-562, May, 2007
multigrid methods;waveform relaxation;wave equation