Analysis of Iterative Line Spline Collocation Methods for Elliptic Partial Differential Equations.
In this paper we present the convergence analysis of iterative schemes for solving linear systems resulting from discretizing multidimensional linear second-order elliptic partial differential equations (PDEs) defined in a hyperparallelepiped $\Omega$ and subject to Dirichlet boundary conditions on some faces of $\Omega$ and Neumann on the others, using line cubic spline collocation (LCSC) methods. Specifically, we derive analytic expressions or obtain sharp bounds for the spectral radius of the corresponding Jacobi iteration matrix and from this we determine the convergence ranges and compute the optimal parameters for the extrapolated Jacobi and successive overrelaxation (SOR) methods. Experimental results are also presented.
Introduction
.
We consider the following second order linear elliptic PDE
subject to Dirichlet and/or Neumann boundary conditions
on
of\Omega (2a)
where Bu is u or @u
is a rectangular domain in R k (the space of
real variables) and ff i (! 0), fi i , fl(- 0), f and g are functions of k variables.
If Line Cubic Spline Collocation (LCSC) methods are used to solve (1a),(2a), then
the differential operator is discretized along lines in each direction independently and
then the line discretization stencils are combined into one large linear system. In
Section 2 we briefly describe such discretization procedures. We also present and
discuss the resulting coefficient linear system of collocation equations.
In [5], we were able to formulate and analyze iterative schemes for solving the
LCSC linear systems in the case of Helmholtz problems, with Dirichlet boundary
conditions and constant coefficients, that is
in\Omega (1b)
on
@\Omega (2b)
Work supported in part by NSF awards CCR 9202536 and CDA 9123502, AFOSR award
F 49620-J-0069 and ARPA-ARO award DAAH04-94-6-0010.
y Purdue University, Computer Science Department, West Lafayette, IN 47907
z Authors permanent address: University of Crete, Mathematics Department, 714 09 Heraklion,
GREECE and IACM, FORTH, 711 10 Heraklion, GREECE.
x u denotes @ 2 u=@x 2 . Unfortunately the convergence analysis presented in
[5], as it stands, can not be applied in the case of the presence of Neumann boundary
conditions on some of the faces of \Omega\Gamma In Section 3 we present a convergence analysis
of block Jacobi, Extrapolated Jacobi (EJ) and Successive Overrelaxation (SOR)
schemes for the iterative solution of the collocation equations that arise from discretizing
elliptic problems (1b) subject to Neumann (N) boundary conditions on one
or more (but not on all if
of\Omega and Dirichlet (D) ones on the others.
More specifically analytic expressions or sharp bounds for the spectral radius of the
corresponding block Jacobi iteration matrix are derived and from these we determine
the convergence ranges and compute the optimal parameters for the Extrapolated
Jacobi and SOR methods. Furthermore, based on our analysis, certain conclusions of
significant practical importance are drawn.
Finally, in Section 4 we present the results of various numerical experiments designed
for the verification of the theoretical behavior of the iterative LCSC solvers.
The experiments show very good agreement with the theoretical predictions as regards
the convergence of the iterative method used. Although the theoretical results
presented here hold for the model problem (1b), (2b), our experimental results indicate
that the behavior of the iterative LCSC solvers on the general problem (1a), (2a)
is similar.
2. The Line Cubic Spline Collocation (LCSC) Method.
In this section we briefly describe the LCSC discretization method and introduce
some notation to be used later. We start by introducing one extra point beyond
each end of the intervals [a i ,b i ]. Each of these enlarged intervals is then discretized
uniformly with step size h i by
ae
oe
The tensor product of these discretized lines,
provides an extended
uniform partition of \Omega\Gamma A collocation approximation u \Delta of u in the space S 3;\Delta of
cubic splines in k dimensions is defined by requiring that it satisfies the equation (1)
at all the mesh points of \Delta and the equation (2) on the boundary mesh points. In the
sequel the interior mesh points are denoted by (- 1
nk ), for all n 2
In the line cubic spline collocation methods we consider in this paper the collocation
approximation is made on each set L i of lines parallel to the x i axis. More
specifically, this collocation approximation u i
on each line in L i is represented
as follows
are the B-spline basis defined on \Delta i ,
that (3) is the sum of one dimensional splines in
the x i variable whose coefficients U i are functions of the other
thermore, this approximation is redundant in that there are k choices for representing
one for each coordinate direction.
2.1. The Second Order Line Cubic Spline Collocation Method.
From the assumed representation (3) of u i
\Delta and the nature of the B-spline basis
functions we conclude easily (see also [8] and [5]) that
and
at the mesh (collocation) points P i
on each line of L i . These
points are represented by the vectors
nk
indexes the L i lines.
First, we observe that the collocation equations obtained from the boundary conditions
and the differential equation at the end points of each line can be explicitly
determined as follows
Dirichlet (D) boundary conditions on P iU i
Neumann (N) boundary conditions on P iU i
Dirichlet (D) boundary conditions on P i
Neumann (N) boundary conditions on P i
where for simplicity we have assumed that all boundary conditions are homogeneous.
The rest of the unknowns are determined by solving the so called interior collocation
equations, that is, on those lines of L i not in @
The equations (11) away (2 - from the boundary are given by
\Theta U i

while for lines next to the boundary, they have similar form with appropriately modified
right sides ([5]). The full matrix form of these equations is given in Section 3.
After eliminating the predetermined (from equations (7)-(10)) boundary unknowns,
each collocation approximation u i
coefficients
' . This redundancy is handled by requiring that all these approximations agree on
the mesh points, that is
\Delta on the mesh \Delta
It has been shown ([5], [6], [8]) that the above described method leads to a second
order collocation approximation of u.
2.2. The Fourth Order Line Cubic Spline Collocation Method.
To derive the fourth order LCSC method we use high order approximations of the
derivatives D j
defined as appropriate linear combinations of the spline
interpolant and its derivatives at the mesh points [6]. Specifically, we approximate
the second derivatives in the PDE operator by the difference scheme
The collocation equations corresponding to the mesh points on a line L i away (2 -
from the boundary are written at the point P i
' as
\Theta U i

The redundancy on the coefficients is handled in the same way as in the second order
case, with the only basic difference being that the stencils now have 5, rather than 3
points along each coordinate direction.
3. Iterative Solution of the Interior Line Cubic Spline Collocation
Equations.
The LCSC equations can be written in the form
where the coefficient matrices are defined by
Y
where I denotes a unit matrix of appropriate order and
The matrix H i depends both on the discretization scheme applied and the nature
(Dirichlet (D) or Neumann (N)) of the boundary conditions on the left and right
end-points of the line direction with which the matrix H
Specifically,
ff;fi for the second order discretization and H
ff;fi for the
fourth order one. The matrices T i
are given by
where
(\Gamma2; \Gamma2) in case of (D) conditions on both ends
(\Gamma1; \Gamma2) in case of (N) conditions on the left and
(D) ones on the right end
(\Gamma2; \Gamma1) in case of (D) conditions on the left and
(N) ones on the right end
(\Gamma1; \Gamma1) in case of (N) conditions on both ends
NOTE: Since in the analysis which will follow, the existence of A \Gamma1
1 is a necessary
requirement, we assume, without loss of generality, that on at least one of the faces
of\Omega perpendicular to the x 1 -direction (D) boundary conditions have been imposed
so that the invertibility of A 1 is guaranteed even if
Equation (16), coupled with the conditions (13) leads us to
the order kK system of linear equations which, if we order according to the ordering
of unknowns U 1 (n), is given
\GammaB 3 D 3
where for we have that
Y
The matrices in (18) and (22), may be linear or quadratic functions of the matrix
\Gamma2;\Gamma2 if all the boundary conditions are Dirichlet (D). With some Neumann (N)
conditions present, this is no longer the case. Besides, in the case of the fourth order
discretization scheme, the product T i
ff;fi is not a symmetric matrix in the case where
(ff; fi) 6= (\Gamma2; \Gamma2). So, the analysis presented in Section 3 of [5] does not always apply
to the present case. We present a unified analysis for all cases based on the following
three lemmas.
Lemma 3.1. If A and B are Hermitian positive definite matrices of the same order
then (i) AB (and BA) possesses a complete set of linearly independent eigenvectors.
(ii) Further, let X 2 C n;n and oe(X) ae R, where oe(:) denotes the spectrum of a matrix.
Let also -X , -X;m and -X;M denote any, the smallest and the largest eigenvalue of
X, respectively. Then for the eigenvalues of the matrix AB there holds:
Proof. Let A 1=2 and B 1=2 be the unique square roots of A and B (see Theorem
2.7, page 22 in [10]). The matrices AB and A 1=2 B A
are similar, with (\Delta) H denoting the complex conjugate transpose of a given matrix.
Obviously, the matrix (B 1=2 A 1=2 ) H (B 1=2 A 1=2 ) is Hermitian and positive definite.
From the similarity property and the fact that A 1=2 BA 1=2 is Hermitian the validity
of the assertion in (i) follows. Now, it is easily seen that
which proves the two rightmost inequalities in (23). On the other hand, since
exists and 1
--AB;m
-B;m
\Delta-A;m
which proves that the three leftmost inequalities in (23) hold. 2
NOTE: The assertions of Lemma 3.1 hold even if one of A or B is nonnegative definite.
The only difference is that the two leftmost inequalities in (23) become equalities, that
is Indeed, the proof for the assertion in (i) is the same if
B is singular, while if A is singular one uses the similarity of the matrices AB and
. The proof for the rightmost inequalities in (23) is exactly the same as
before. For the leftmost ones we simply observe that
Lemma 3.2. Consider the real symmetric tridiagonal matrix T ff;fi of order m
given by
where (ff; Its eigenvalues - T ff;fi are given
by the expressions
(i) If (ff;
(ii) If (ff;
(iii) If (ff;
Proof. The results (25i) and (25iii) are well-known in the literature but we give
here a unified way of obtaining all three of them simultaneously. First, we note that
it can be proved that all these matrices are negative definite, with the exception of
T \Gamma1;\Gamma1 which is non-positive definite. We have ae(T ff;fi ) -k T ff;fi k1= 4 and since it
can be checked that \Gamma4 62 oe(T ff;fi ) we conclude that oe(T ff;fi) ae (\Gamma4; 0]. To determine
all the non-zero eigenvalues - of T ff;fi let denote the associated
eigenvectors. From T ff;fi z = -z, one obtains
z
The set of the above equations can be written as
z
provided one sets
The characteristic equation of (26) is
and, since we are looking for - 6= 0; \Gamma4, the two zeros of the quadratic in (28) such
that r 1 6= r 2 satisfy
Consequently the solution to (26) is given by
If we arbitrarily put z use the restriction on z 0 from (27i), the coefficients
c 1 and c 2 can be determined. Next, using the restrictions on z m+1 from (27ii), and
(29ii), one arrives at
r 2m
So, we consider the three cases of the lemma:
m+1 and
from (29ii) one obtains the m distinct expressions for -(6= 0; \Gamma4) given in
(25i).
cos
and again from (29ii) we have the m distinct
values for -(6= 0; \Gamma4) in (25ii).
Working in the same way we obtain the
expressions
\Gamma4). If we incorporate the
as well we have the expressions in (25iii). 2
Lemma 3.3. Let the n \Theta n matrices A possess complete sets
of linearly independent eigenvectors y (i;j) , with corresponding eigenvalues - (j)
k). Then the matrix A j A
1\Omega A
:\Omega A k possesses the
linearly independent eigenvectors y (j) j y (1;j1
)\Omega y (2;j2
corresponding eigenvalues -
.
Proof. First, using tensor product properties we can easily verify that Ay
In order to prove the linear independence of the y (j) 's we construct the
whose columns are the n k eigenvectors of A in the following order
(1;1)\Omega y
y
(1;1)\Omega y
(1;1)\Omega y
y
(1;2)\Omega y
y
(1;2)\Omega y
(1;2)\Omega y
(1;n)\Omega y
\Theta y (i;1) y

. From the assumed linear independence
of y (i;j) , we conclude that Y \Gamma1
k exist. This implies the linear independence of the
eigenvectors and concludes the proof of the lemma. 2
Having developed the background material required we are able to go on with the
analysis of the three methods associated with the linear system (21): block Jacobi
(J),
block Extrapolated Jacobi (EJ),
and block Successive Overrelaxation (SOR)
where
and \GammaL and \GammaM are the strictly lower and the strictly upper triangular parts of the
matrix coefficient in (23).
The block Jacobi iteration matrix J associated with the
matrix coefficient in (21), can be described as
Let H and G denote the matrices
\Theta G T
of dimensions K \Theta resectively. Consider then the matrix
apparently, is a block diagonal matrix with diagonal blocks HG and GH
of orders K and respectively. Therefore we have (see Theorem 1.12, page
in [10]) that
However, we have
and each term in the above sum can be found by using (32), (17), (22) and simple
tensor product properties to be
For the second order LCSC, we introduce the matrices S j and find their value to be
6 I
For the fourth order LCSC, we have similarly
\Gamma2;\Gamma2
\Gamma2;\Gamma2
6 I
Due to the presence of the unit matrix factors in the tensor product form (36), all
possess a linearly independent set of K common eigenvectors if and only if
each H
possess complete sets of linearly independent
eigenvectors. For the matrices in (38) this is obvious because T j
real symmetric negative (or non-positive) definite, T 1
ff;fi is real symmetric negative
definite and, in addition, oe(T j
As one can readily see, each of the matrices S j in (39) is
the product of the two real symmetric positive definite matrices 1
\Gamma2;\Gamma2 ) and
(The second matrix factor might also be nonnegative definite.)
Therefore Lemma 3.1 (or its Note) applies. For the matrix A \Gamma1
4 in (40) first observe
that the first matrix term in the brackets is similar to12 (12I
\Gamma2;\Gamma2
which, in turn, is of exactly the same form as the matrices S j in (39), except that
the second matrix factor considered previously is now always positive definite. Con-
sequently, in all possible cases of the LCSC matrices, all terms H j G j in (36) possess
a linearly independent set of K common eigenvectors. By virtue of this result and in
view of Lemma 3.3, it is implied from (36) that
where -X is used to denote any eigenvalue of the matrix X. However, from the
previous discussion it follows that - S i - 0,
-HG - 0. So, from (34) it is implied that J 2 possesses non-positive eigenvalues and
hence the block Jacobi matrix J has a purely imaginary spectrum. From the analysis
so far it becomes clear that the eigenvalues of J 2 and therefore those of J can be given
analytically in the following cases:
(1) In all the cases of the second order LCSC we are dealing with when Dirichlet
and/or Neumann boundary conditions are imposed on the faces of @
(2) In the fourth order one when only Dirichlet boundary conditions are imposed on
the faces of @
This is an immediate consequence of the fact that each matrix S j ,
(37)-(40) is a simple real rational matrix function of the matrix T j
case (1) and of the matrix T j
\Gamma2;\Gamma2 in case (2), respectively.
Having in mind the various conclusions we have arrived at in the analysis so
far, one can state the following theorem which gives analytic expressions for the
eigenvalues of the block Jacobi matrix J in (32).
Theorem 3.4. The eigenvalues - of the block Jacobi iteration matrix
defined in (32) are as follows. We have
For the second order collocation scheme the others are pure imaginary
For the fourth order scheme, where (2a) is assumed to be subjected to Dirichlet boundary
conditions only, the others are pure imaginary
\Gamma2;\Gamma2
\Gamma2;\Gamma2
\Gamma1=2
In (42) and (43) - j
are the eigenvalues of the matrix T j
obtained by the
application of Lemma 3.2. 2
NOTE: It should be pointed out that analytic expressions for the eigenvalues of
J can not be derived, in terms of the matrices T j
ff;fi involved, in the case of the fourth
order LCSC where on at least one face of
@\Omega has Neumann (N) boundary conditions
imposed. This is due to the non-commutativity of the matrices T j
\Gamma2;\Gamma2 and T j
ff;fi for
(ff; fi) 6= (\Gamma2; \Gamma2). However, one can trivially give analytic expressions in terms of the
eigenvalues of the matrices (12I +T j
also, by
virtue of Lemma 3.1, lower and upper bounds can be given for the eigenvalues of H j
in terms of the extreme eigenvalues of T j
\Gamma2;\Gamma2 , and T j
of Lemma 3.2.
To derive the spectral radii of the Jacobi matrices of Theorem 3.4 and also upper
bounds in the case of the previous Note, we introduce some notations first. Let, then,
ff;fi denote an eigenvalue of the matrices T j
as in Theorem 3.4. To
simplify the notation we omit using another index on the generic - j
ff;fi and we omit
the index j from the pair of subscripts (ff; fi). Let also
denote the smallest and the largest eigenvalues in - j
Finally, define the functions
\Gamma2;\Gamma2 )y j
in terms of the eigenvalues - j
ff;fi of the matrices T j
k. Then we have:
Theorem 3.5. (i) The spectral radii of the block Jacobi matrices corresponding
to the collocation schemes considered in Theorem 3.4 are given by the following
expressions:
Second order.
Fourth order.
(ii) Moreover the expression (48) is also a strict upper bound for the square of the
spectral radius of the Jacobi matrix in the case of the fourth order scheme corresponding
to an elliptic problem (1b) where on at least one of the faces of
@\Omega has Neumann
(N) boundary conditions imposed.
Proof: We notice that in (42) and (43) ff
and Hence the functions y
in (45) are nonnegative
and independent of each other. So are the expressions z j := z j (- j
z(-) in (46). To determine the spectral radius of the block Jacobi matrix of the
second order collocation scheme we examine the extreme values of y
k. By differentiation we have
@-
and, consequently,
For the corresponding quantity of the fourth order scheme of Theorem 3.4, working
in a similar way, we obtain
@-
which implies that
It is obvious now using (42) - (46) and (49), (51), that the results (47) and (48) follow.
This concludes the proof for part (i) of the theorem.
For part (ii) we use the analysis preceding the statement of Theorem 3.4 and refer
to the matrices in (39) and (40), especially for those indices
(ff; fi) 6= (\Gamma2; \Gamma2), and also recall the Note immediately following Theorem 3.4. It
follows from this and Lemma 3.1, its Note and Lemma 3.2, that the lower and upper
bounds for the eigenvalues of the matrices in (39) and (40) depend directly on the
extreme eigenvalues of the two
\Gamma2;\Gamma2 ) and ff j
These are readily seen to be12 (12
for the minimum and the maximum eigenvalues for the former matrix and y j
for the latter one, where c j , s j and y j are the expressions in (44) and (45).
From Lemma 3.1, its Note, and using the expressions (46), the bound for the spectral
radius of the corresponding block Jacobi matrix is readily shown to be the expression
in the right hand side of (48).2
REMARK: The analysis so far has been made on the assumption that the x 1 -
direction is somehow predetermined. However, since there are k possible choices for
the x 1 -direction in any particular case, the choice should be made in such a way as
to give the smallest possible values in (47) and (48). Consider then the quantities
min
iB @
taken over all which they are well-defined and also the quantities
min iB @
z
z
C A
considered in the same way. The indices i in (51) and (52) for which the corresponding
minima take place should be interchanged with the index 1 and therefore the x i -
direction should be taken as the x 1 -direction. If, in either case (51) or (52), more
than one index i gives the same minimum value then any such i will do. It should
be noted that after having chosen the x 1 -direction this way, the expressions in (47)
and (48) give then the smallest possible values for the spectral radius or for an upper
bound on it, as was explained, for the particular problem at hand.2
From this point on, the analysis is almost identical with that in [5] and the
interested reader is referred to it. For completeness, we simply mention some of
the theoretical results obtained in [5], which are based on the corresponding theory
developed in [10], [11] and in [1], [2], [3], [4], [9].
(i) The block Extrapolated Jacobi (EJ) method and the block SOR method corresponding
to the block Jacobi method of this paper converge for values of their
parameters varying in some open intervals whose left end is 0 and the right one is a
function or ae(J). However, the optimum SOR is always superior to the optimum EJ
and the corresponding optimal parameters for the SOR method are given by
(ii) For the efficiency of both the serial and the parallel iterative solution of the
linear system (21), a cyclic natural ordering of the unknowns U i ,
adopted according to which
The equations in (16) are reordered according to the ordering of U 1 and each block
of the auxiliary conditions is reordered according to the ordering of the unknowns
U i . It is then proved that the new coefficient matrix A is obtained by a permutation
similarity transformation of the matrix coefficient A in (21) having the same k \Theta k
block structure and, therefore, the associated block Jacobi matrix J is similar to the
previous one J . Consequently, the convergence results are identically the same so that
all the formulas in connection with eigenvalues, spectral radii, etc. of the Jacobi, the
Extrapolated Jacobi and the SOR method studied in this section remain unchanged
when these reorderings are made.
(iii) The new structure of the collocation coefficient matrix A, for the second and
fourth order scheme in 2-dimensions, is given in [5].
4. Numerical Experiments.
In this section we summarize the results of some numerical experiments that verify
the convergence properties of the iterative solution methods analyzed in Section 3. We
mention that although we present numerical data for only the O(h 2 ), 3-dimensional
case, these are very representive of problems with different dimensionality or with
discretization schemes. For experimental data on the convergence properties
of the LCSC method and the iterative solvers in the case of only Dirichlet boundary
conditions the reader is referred to [5] and [6]. The parallel implementation details
and the performance of the iterative LCSC schemes, for 2-dimensional problems, on
several multiprocessing systems can be found in [7].
We have applied the LCSC discretization techniques with uniform (NGRID by
NGRID by NGRID) meshes to approximate the known solution of the following set
of PDEs.
x
z
PDE 2: 4D 2
y
z
PDE 3: D 2
9
z
Each is defined on the unit cube and subject to one of the following types of boundary
conditions.
1: Dirichlet conditions on all faces of \Omega\Gamma
2: Neumann conditions on the faces
of\Omega and Dirichlet ones on
the rest.
3: Neumann conditions on the face
of\Omega and Dirichlet on the rest.
These PDEs and boundary conditions are combined to give 9 problems, our theory
is applicable only to 6 of these (those excluding PDE 3). The right hand side f is
selected so that the true solution is always
The linear systems from the LCSC discretization were solved by the proposed SOR
iteration method with the termination criterion being that jjU

Table
The required number of SOR iterations to solve the LCSC equations for PDE 2 and various
boundary conditions.
the interval (0; 10 \Gamma7 ). All experiments were performed, in double precision, on a
workstation.
In

Figure

1 we present the theoretically estimated (using the material developed in
Section 3 when applicable) and the experimentally determined (by systematic search),
values of the optimumSOR relaxation parameter ! opt . Specifically, the points we plot
are the experimentally observed optimal values of ! for various values of NGRID.
The lines we plot show the relation (determined using (53) and Theorem (48) ) between
opt and the discretization parameter NGRID. Since our theory can not be
directly used to determine such relation in the case of PDE 3 we plot only the experimentally
determined ! opt .
Our first observation is that there is a good agreement between our theory and the
experiments in the six cases where it applies. The theoretical values for ! opt are close
to (though always larger than) the measured ones and exhibit the same dependence
on the discretization and PDE problem parameters. Besides confirming our theory,
these experiments also show that for PDEs where our analysis is not applicable the
proposed SOR scheme still converges. Furthermore relaxing with ! opt determined by
our theory leads us to comparable rates of convergence.
Another interesting observation is that ! opt seems to go fast and asymptotically to
a number in the interval [:03; :08]. It is therefore expected that the rate of convergence
will not decrease rapidly as NGRID increases beyond 30. This is confirmed in Table
1 where we observe that increasing NGRID from 24 to 32 increases the number of
iterations only by at most 30%.

Table

1 presents the SOR iterations required to solve the discretized equations
using the optimal value for the relaxation parameter ! and a 10 \Gamma7 stopping criterion
for the problems defined by PDE 2 and various boundary conditions. We also note here
that the measured discretization errors for all the experiments confirm the expected
second order of convergence of the collocation discretization scheme. Specifically, the
measured order in all cases is in the interval [1:8; 2:1].



--R

Iterative methods with k-part splittings
The optimal solution of the extrapolation problem of a first order scheme


Iterative line cubic spline collocation methods for elliptic partial differential equations in several dimensions
Spline collocation methods for elliptic partial differential equations

Convergence of O(h 4
On complex successive overrelaxation
Matrix Iterative Analysis
Iterative Solution of Large Linear Systems
--TR
SOR iterative method;collocation methods;elliptic partial differential equations