Integration preconditioners for di erential operators in spectral -methods E. A. Coutsias
T. Hagstromy
J. S. Hesthavenz
D. Torresx
Abstract
1 Introduction
We present simple, banded preconditioners that transform linear ordinary di erential operators with polynomial coecients into banded form. These are applicable to a wide class of Galerkin approximation problems, including expansions in terms of all the classical orthogonal polynomials. The preconditioners are in fact the n-th order integration operators for the polynomial families employed in the Galerkin approximation, with n the order of the di erential operator. The resulting matrix problems are algorithmically simpler, as well as better conditioned than the original forms. The good conditioning allows the extension of our ideas even to problems with arbitrary, nonsingular coecients as well as to certain quasilinear problems by the use of iterative methods. We also present extensions to partial di erential operators with polynomial coecients by considering preconditioners in the form of tensor products of appropriate combinations of integration operators. The origin of the tridiagonal integration operators for arbitrary classical orthogonal polynomial families is shown to lie with the Gauss contiguity relations for Hypergeometric functions. Key words: spectral methods, orthogonal polynomials, boundary value problems. AMS subject classi cations: 65Q05, 65L60, 65P20, 76-08, 33A45, 33A50, 33A65.
The classical orthogonal polynomial bases most commonly used in Numerical Analysis originate as eigenfunctions of singular Sturm-Liouville problems. The derivatives of such polynomials form an orthogonal basis as well; in fact they are also classical orthogonal polynomials. As a result of the Gauss contiguity relations for Hypergeometric Functions [1], elements of the original basis have a simple expression in terms of elements of the derivative basis, involving at most three terms of contiguous degrees. Consequently, although the matrices representing di erentiation to various orders in terms of an orthogonal polynomial basis are almost full upper triangular matrices (only exception are the Hermite polynomials), those representing integration are banded, of bandwidth 2n + 1, with n the order of integration. This, together with the fact that multiplication by a monomial is also a tridiagonal matrix leads to simple, banded preconditioners which simultaneously band classes of matrices of the form PDm , with P the operator of multiplication by a polynomial p(x) and Dm the operator of m-fold di erentiation. This allows the formulation of ef cient algorithms for the solution of di erential equations with polynomial coecients in any of the classical polynomial bases. In Sec. 2 we establish the above facts and discuss how they can be used in the context of the Lanczos -method to construct a spectrally accurate solution of a di erential equation with general polynomial (and rational) coecients in O(M) operations with M the truncation order. The appeal of the method is due not only to the eciency of solution, but also to the excellent conditioning of the resulting matrix problems [6]. This latter property, established under certain general assumptions in Sec. 3 permits the extension to problems with arbitrary nonsingular coecients as well as certain nonlinear problems by the use of iterative methods, and we present examples in Sec. 4. In Sec. 5 we show how the ideas can be extended to higher dimensional problems. The pre-conditioning scheme presented here has the advantage that it can be extended easily to treat problems in multidimensions. Additional dimensions are included in the banded matrix formulation through tensor products, which amounts to re-
Dept. of Mathematics and Statistics, University of New Mexico, Albuquerque, NM 87131. E-mail:
[email protected] y Dept. of Mathematics and Statistics, University of New Mexico, Albuquerque, NM 87131. E-mail:
[email protected] z Association EURATOM - Ris National Laboratory, Optics and Fluid Dynamics Department, P.O. Box 49, DK-4000, Roskilde, Denmark. E-mail:
[email protected] x Dept. of Mathematics and Statistics, University of New Mexico, Albuquerque, NM 87131. E-mail:
[email protected] This space left blank for copyright notice.
1
ICOSAHOM 95
2 placing entries in banded matrices by blocks which are themselves banded. A natural way of deriving such blockbanded forms for the higher dimensional case is by interpreting integration and di erentiation as change of basis transformations among related polynomial families. Effectively, then, we consider expansions in each variable in terms of a basis of derivative polynomials of order equal to the maximum order derivatives in that variable present in the given equation. The numerical solution of several simple test cases is presented. As an alternative to nitedi erence based time stepping schemes for time dependent problems, we apply our method to problems in one space dimension plus time. Functions are expanded in a double spectral expansion, in both the space and time variables. The problems are solved by inversion of the block-banded matrix resulting from the application of integration preconditioners of appropriate order in each variable. The operation count for a problem discretized to a N M resolution is min(O(M 3 N); O(N 3M)) for a single solution but it is lower if the same matrix is inverted several times, say as a result for solving with di erent forcing functions etc, in which case the cost is min(O(M 2 N); O(N 2M)). These numbers re ect the alternatives available when deriving the block-banded forms for the di erential operator. Analogous estimates also hold for higher dimensional problems. Finally, in Sec. 6 we present some concluding remarks and further connections with previous works. In a previous paper [6] we examined the use of the spectral integration operators as postconditioners for linear differential operators with polynomialcoecients in arbitrary bases of classical orthogonal polynomials. In that form, the method was a generalization of the method of treating a di erential equation by expanding the highest derivative in terms of Chebyshev polynomials originally introduced by Clenshaw [4]. The integration preconditioner for the Chebyshev polynomials is also analyzed by Greengard [11], while the recurrence relation for the derivatives of the Jacobi polynomials has also appeared in a recent review by Fornberg [7].
2 The method
will assume that the functions under consideration possess sucient di erentiability properties over (a; b) and can be expressed as a series involving the Qk . See [3] for a discussion of the convergence properties and the introduction of relevant function spaces.
2.1 Integration operators and derivative bases R
We write hQl ; Qm iw = ab Qm Ql w(x)dx = hm lm , with hmn the norming constants of the Qm [1]. We set Qn(x) = P Knm xm . It is well known that all orthogonal poly0 nomial families share a three-term recurrence of the form [1]: 1 X (2) Qk+l ak+l;k = xQk ; k = 0; 1; l=
1
This follows since, if l < k 1, deg(xQl ) = l + 1 < k and, by orthogonality, hxQk ; Ql iw = hQk ; xQl iw = 0. By matching powers we easily see that [14]: an;n 1 = KnK 1;n 1 ; an;n+1 = hhn+1 an+1;n; n;n n K Kn+1;n n;n 1 an;n = K (3) K n;n
n+1;n+1
the second relationship following since hn+1 an+1;n = hxQn; Qn+1iw = hxQn+1; Qniw = hnan;n+1: In many important cases, including the classical orthogonal polynomials (i.e. Jacobi, Chebyshev, Legendre, Gegenbauer, Hermite and Laguerre polynomials) there also holds a relation of the form (4)
X 1
l=
1
Q0k+l bk+l;k = Qk ; k = 0; 1; :
Here, as well as in (2), we introduced Q 1 = 0. In fact:
Lemma 2.1 Suppose that the eigenfunctions of the SL problem (1) are polynomials. Then their derivatives, fQ0k g1 1 also constitute an orthogonal family with respect to weight p(x) where it is assumed that p(x) > 0 for x 2 (a; b) and p(a) = p(b) = 0. Moreover, the expression of Qk (k > 1) in terms of the Q0k involves at most only fQ0k gkk+11, i.e. it has the form (4).
Throughout, we assume that we are working with a family of polynomials fQk g1 0 which are orthogonal and complete over the interval (a; b) (here a and/or b can be in nite) with respect to the weight w(x). In the cases of interest, these Proof: Integrating by parts we see Zb Zb are the eigenfunctions of a Sturm-Liouville (SL) problem, b 0 0 0 Qk Ql p(x)dx = Qk Ql pja Qk (Q0l p)0 dx (1) (p(x)Q0k )0 + k w(x)Qk = 0; a a Zb so that the Q0k form an orthogonal family as well, with = k Qk Ql wdx = k hk kl weight p(x) which satis es p(x) ! 0 as x ! a; b. We a
Integration Preconditioners
3
since p(x) vanishes at the end-points, and the orthogonal- and ! ity of the Q0k follows from that of the Qk . We introduce nX1 0 0 1 K 1 1 ; 0 hk 1 = k hk for the norming constants of the Qk , indexing al;l K : bn;n = n + 1 an;n n(n + 1) according to the degree of the corresponding polynomial, 1;1 l =1 k +1 namely deg(Q0k ) = k 1. Now, P since the fQ0l g1 are independent, we have that Qk = k1 +1 Q0l bl;k and in or- Proof: Since the Q0k are orthogonal, they must satisfy a der to establish (4) we must show that hQk ; Q0l ip = 0 for relation of form (2): l = 1; : : :; k 2. Clearly 1 X Z (7) Q0k+l+1 a0k+l;k = xQ0k+1 ; k = 0; 1; b b 0 0 0 0 0 0 hQk ; Ql ip = Qk Ql pxja (xQk Ql p + Qk x(Ql p) ) dx l= 1 a Pn K 0 xm = Pn(m + Since Q0n+1 (x) = = hQ0k ; xQ0lip + l hQk ; xQliw = 0 m=0 nm 0 0 = 1)Kn+1;m+1 xm , i.e. Kn;m (m + 1)Kn+1;m+1 , it follows where the rst integral vanishes due to the orthogonality from (3) that of the Q0k since deg(xQ0l ) = l k 2 < k 1 = deg(Q0k ), + 1a while the second integral vanishes since deg(Qk ) = k > a0n+1;n = nn + 2 n+2;n+1 ; k 1 l + 1 = deg(xQl ). Note: for the classical orthogonal polynomials, the derivan n+1hn+1 a tives of arbitrary order are also classical orthogonal polya0n 1;n = n + 1 nhn n+1;n; nomials. This follows easily from general properties of the n + 1 Kn+2;n+1 : Hypergeometric functions F(a; b; c; z) [1]. For instance, for 0 = n Kn+1;n a n;n the Jacobi polynomials we have n + 1 Kn+1;n+1 n + 2 Kn+2;n+2 Di erentiating (2) and using (7) there results: Pn( ; )(1 2z) = ( +n!1)n F( n; n + + + 1; + 1; z) : 1 X Q0k+l ak+l;k a0k+l 1;k 1 ; Qk = Since l= 1 d F(a; b; c; z) = ab F(a + 1; b + 1; c + 1; z); (5) dz and the claim follows after some algebra. c Note: The computation of the bi;j from the above expresit follows that sions is not very practical, and was only given to establish the connection to the ai;j . A more direct calculation in d P ( ; )(x) = n + + + 1 P ( +1; +1)(x) : (6) dx n n 1 the case of the Jacobi polynomials follows from their con2 nection to the hypergeometric functions [14]. Indeed, difSimilarly, for the Laguerre and Hermite polynomials we ferentiating (2) for the Jacobis, multiplying (6) by x and have expressing the latter in terms of the Pk( +1; +1) using (2) d d ( ) we nd ( +1) dx Ln (x) = Ln 1 (x) ; dx Hn = 2nHn 1: n + + + 1 a( +1; +1) ; ) ( ; ) b(n ; +1;n = an+1;n In all cases, di erentiation can be seen as a change of basis, n + + + 2 n;n 1 always within the set of classical orthogonal polynomials. ; ) = a( ; ) a( +1; +1) ; We let A, B be the coecient matrices in (2,4) respectively. b(n;n n;n n 1;n 1 The relation between the two sets of coecients is found + + 1 ( +1; +1) from the following: b(n ; 1;n) = a(n ; 1;n) n + n + + an 2;n 1 : Lemma 2.2 The coecients bm;n in (4) are found from In Appendix A we show that similar relations for the those in (2) as derivatives resulting in a tridiagonal integration operator (as well as a tridiagonal monomial multiplication) hold for bn+1;n = n +1 1 an+1;n ; all hypergeometric and con uent hypergeometric functions as a result of the Gauss contiguity relations, and this can (n 1) be used to give an alternative direct derivation of the ren an 1;n; bn 1;n = 1 n currence coecients ai;j and bi;j . n
1
ICOSAHOM 95
4 The nonzero elements of the matrices A, B for the classical orthogonal polynomials are given in Table 1, together with other relevant quantities, using the standard notation [1]. The relations for the Gegenbauer polynomials Cn( ) can be constructed from those of the general Jacobis since + 1) (2 + n + 1) ( ; ) Cn( ) = ( ( + n + 1) (2 + 1) Pn where = = 1=2.
2.2 Reduction to banded form via integration preconditioning Our main result can be stated as follows:
n rows which after left multiplication by B[nn] are null, are replaced by row vectors associated with the -constraints. These alter the matrix but do not a ect the order of complexity of the solution algorithm, apart from an increase in the bandwidth which remains R + n. We establish the above result in a series of lemmas. We de ne the function evaluation functional at the point x, qx, as the row vector qx = (Q0(x); Q1(x); ) : Also we introduce qx(k) = Q(0k)(x); Q(1k)(x); ; the operator of evaluating the k-th derivative. Clearly, the
Theorem 2.1 Consider the orthogonal polynomial family fQnk g1 form an independent set of polynomials which, fQng1 and assume that the Qk satisfy (4) for some matrix for thek classical orthogonal polynomials, can be shown ( )
0
B = (bij ). Then, the matrix representing the di erential operator
(8)
L=
n X
k=0
pk (x)Dk
with deg(pk ) = k , the degree of the polynomial coecient of Dk , becomes banded upon left multiplication by B[nn] with bandwidth R = maxk (2k +2(n k) +1) where 0 k n.
(In the sequel we use L both for the operator and for its matrix representation. Also, arrays are indexed from 0 to N, the maximum order of truncation. An array G whose rst k rows have been replaced by zeroes is denoted by G[k] . Similarly f^[i] will denote the vector f^ with its rst i components set to zero.) Example: The following ordinary di erential equation arises when one solves the 2-D Helmholtz equation on an annulus: 2 k 1 2 u(x) D + x + a D (x + a)2 I u(x) = f The inner radius of the annulus is a 1 > 0 ,k is an integer (representing the Fourier mode), and the range of x is 1 x 1. Multiplying through by the factor (x+a)2 we arrive at an equation with polynomial coecients which can be transformed to nine-diagonal form via left multiplication 2 by the matrix B[2] [5]. This example is further discussed in the next section, where conditioning is considered. Note: In conjunction with the Lanczos -method [10], or alternatively by making use of proper subspaces where di erentiation is invertible [6], the above idea can be incorporated into the design of algorithms for the ecient and accurate solution of di erential equations with polynomial coecients. When the -method is used, the rst
to be orthogonal with a weight related simply to w(x). If we write f^(0) f^ = (f^0 ; f^1; )T for the vector of the expansion coecients of a function f(x) in the given basis and f^(k) = (f^0(k) ; f^1(k); )T for the vector of the expansion coecients of the function f (k) (x) (denoting the k'th derivative of the funcion f(x)) we have that f (k) (x) = qx(k)f^(0) = qx f^(k) . Then the relations (2), (4) can be written as (9) xqx = qxA and qx = qx(1)B = qx(1) B[1] where A, B are the recursion coecient matrices for (2), (4) respectively. The matrix B has the form 2 b0;0 b0;1 0 0 : : : 3 66 b1;0 b1;1 b1;2 0 : : : 77 7 6 (10) B = 66 0 b2;1 b2;2 b2;3 .0 77 : 64 0 0 b3;2 b3;3 . . 75 0 0 0 ... ... The very rst row of B will always be set to zero. Thus the elements b0;0 and b0;1 are irrelevant and will never be used. The basic recursion (4) remains unchanged regardless of how the rst row of B is de ned since Q(1) 0 (x) = 0 for polynomial families. We have Lemma 2.3 The operator Pn of multiplication by a polynomial pn(x) is expressed in the basis Qn by a banded matrix, of bandwidth 2n + 1. Proof: Consider an expansion f(x) = qxf^(0) . Multiplying by x, one obtains xf(x) = xqxf^(0) where xqx denotes
Integration Preconditioners xqx = (xQ0(x); xQ1(x); ). Invoking (9) one can substitute for xqx to obtain xf(x) = qxAf^(0) . Thus A is the matrix which transforms the vector of expansion coecients for f(x) into the vector of expansion coecients for xf(x). The matrix A is tridiagonal since the recurrence relation (2) is a three-term recurrence relation. Extending this argument, it is evident that An is the matrix that transforms the vector of expansion coecients for f(x) into the vector of expansion coecients for xnf(x). Since A was tridiagonal An is obviously banded with bandwidth 2n + 1. Thus multiplication by a polynomial pn (x) becomes the operation of multiplication by pn(A). Note: If the family has a simple convolution, as is the case for the Chebyshev polynomials for which 2Tm Tn = Tm+n + Tjm nj , then it is convenient to expand the polynomial pn(x) in terms of the basis thus simplifying the constructionPofk=the matrix Pn. Speci cally, one rst expands pn(x) = k=0n ck Tk . For purposes of implementation, let us consider a truncated Chebyshev expansion of order N, that is f(x) = TxN f^(0) where TxN is the row vector TxN = (T0 (x); T1(x); ; TN (x)). The product is pn(x)f(x) = P k=n c T T N f^(0) where T T N represents the vector k x k k x k=0 Tk TxN = (Tk (x)T0 (x); Tk (x)T1 (x); ; Tk (x)TN (x)). The vector Tk TxN can be expressed in the form (Tk TxN ) = TxN ATk (accurate up to the N+1'th coecient) where ATk is the N + 1 by N + 1 matrix 20 0 1 0 03 .. 7 66 ... 0 1 0 . . . . . . . 66 0 1 0 . . . . . . . . . ... 777 66 . . . . . . 0 77 1 2 0 77 6 ATk 2 66 . . . . . 1 77 66 0. . 1 . . . 7 66 .. . . . . . . 0 77 .. 5 ... ... ... 4 ... . 0 0 1 0 0 where ak0 = 1, a0k = 1=2 etc, which follows from 2Tm Tn = Tm+n + Tjm nj . Note AT0 reduces to the idena bandwidth of 2k + 1. tity matrix. As expected, P =nAc TTk has N T f Thus pn (x)f(x) = kk=0 k k x ^(0) can be written as P P =n c A =n pn(x)f(x) = TxN kk=0 ck ATk f^(0) where kk=0 k Tk is sum of matrices ATk each with bandwidth less or equal to 2n + 1. This equation illustrates that multiplication by pn(x) translates into multiplication by a banded matrix for Chebyshev polynomials. We also have the related obvious consequence of Leibnitz's rule: Lemma 2.4 The commutator of the operator Dk with the operator P of multiplication by a polynomial p(x) is given
5 by
k P; Dk PDk Dk P = X ( 1)m mk Dk m P m : (
)
m=1
The properties of the B[kk] are established in lemmas (2.52.6): Lemma 2.5 Let f 0 (x) = Pk f^k(1) Qk (x); then B[1] f^(1) = (0) and the rst element of f^(0) , f^0(0) , is undetermined. f^[1] Proof: We have that f(x) = qx f^(0) and f 0 (x) = qx(1) f^(0) = qxf^(1) . Also, by assumption, qx(1)B = qx. Combining we nd qx(1) f^(0) B f^(1) = 0: Since the Q0k are independent and orthogonal, the relation claimed follows. Clearly, the rst element of f^ remains undetermined, since Q00 0. Corollary 2.1 A similar relationship holds for the n-th derivative coecients of f(x). Recall that qx(1)B = qx from (9). Similarly, one has qx(i+1) B = qx(i) . One can generate, using repeated applications of qx(i+1) B = qx(i) and recursive substitutions,
qx(n)B p = qx(n p) where p n. Setting n = k and p = k, in (11) one obtains qxk B k = qx. Arguing as in Lemma 2.5 it is seen that f (k) (x) = qx(k)f^(0) = qxf^(k) = qx(k)B k f (k) using qxk B k = qx for the last equality. Using the second (11)
and last expression in the above equation, one obtains
k ^(k) qx(k) B k f^(k) f^(0) = 0 ! f^[(0) k] B[k] f = 0 and, again, this relation gives the coecients of f(x) in terms of those of its k-th derivative, with the components f^0 through f^k 1 remaining unde ned. Recall that the rst k components of qx(k) are zero, since the k'th derivative annihilates the polynomials of order rst through k 1. We then have Lemma 2.6 If D is the matrix of di erentiation in the family, then B[1] D = I[1] , the identity with its rst row set
to zero. Proof: Indeed, by de nition, Df^(0) = f^(1) , so that
qxf^(1) = qx Df^(0) from which qx(1) B[1] (f^(1) Df^(0) ) = 0 using (9). Thus qx(1) (B[1] f^(1) B[1] Df^(0) ) = 0 or (0) B[1] Df^(0) ) = 0 by Lemma 2.5. One can now qx(1) (f^[1] write (I[1] B[1] D)f^(0) = 0 thus concluding the proof. It then follows that:
ICOSAHOM 95
6
Corollary 2.2 If Dn is the matrix representing n-fold difl n ferentiation in the family Qn , then B[l++n] Dn = B[ll+n] . Proof: By the de nition of the di erentiation matrix D.
(12) qx(k)f^(0) = qx Dk f^(0) Setting p = n in (11), one has qx(n)B n = qx. Substituting this expression for qx into (12), one obtains (13) qx(k)f^(0) = qx(n)B n Dk f^(0) : Setting p = n k 0 in (11), one has qx(n)B n k = qx(k). Substituting this expression for qx(k) in (13), one obtains (14) qx(n) B n k f^(0) = qx(n) B n Dk f^(0) : Now the rst n components of qx(n) are zero. However, the remaining components of qx(n) are composed of orthogonal polynomials and are therefore independent. Thus, one has B n k f^(0) = B n Dk f^(0) except for the rst n components. Since f (0) is an arbitrary vector, one has B n k = B n Dk except for the rst n rows. This equation can be stated concisely in the form B[nn] k = B[nn] Dk or letting n = n + l and k = n, in the form B[ll++nn] Dn = B[ll+n] . We now give the proof of our main assertion: Proof of theorem: Following lemma (2.4) we rewrite the di erential operator L as k n X X ( 1)m mk Dk m Pk(m) L = k=0 m=0 = =
n n X rX r=0
n X
D
k=r
(
1)k r
k
(k r ) r Pk
Dr Sr
r=0 with deg(Sr ) = r = maxrkn(k k + r). by cor.(2.2) B[nn] Dr = B[nn] r , n X B[nn] L = B[nn] r Sr r=0
and the bandwidth is obviously 2 0max ((n r) + r ) + 1 = R : rn
Then, since
Replacing the rst n-rows (containing zeroes) by appropriate constraint coecients, originating from the boundary, initial or other conditions imposed on the solution of the ODE Lu = f, will transform the matrix into a banded form with n additional (generally nonzero) rows. This matrix still factors similarly to a banded matrix, with bandwidth R + n.
3 Conditioning and Convergence It is well known that spectral di erentiation operators suffer the dual defects of full upper triangular matrix representations with very poor conditioning. Above we have shown how integration preconditioners band spectral differentiation matrices. In this section we show that they produce well-conditioned linear systems, and exploit that fact to give good error estimates for general ordinary di erential equations. In the following section we will show that the favorable conditioning properties lead to the rapid convergence of iterative methods for spectral approximations to general variable coecient and nonlinear equations. In particular we study the preconditioned, discrete system: (15) TN v^ = b; (16)
1 0 nX n j @I n + B n Sj A v^ = Bnn FN ; 1
[ ]
j =0
[ ]
[ ]
under the simplifying assumption that the lead coecient is 1. Here TN v^ = T v and the matrices Sj are Galerkin approximations to multiplication by the polynomials, sj (x). Throughout we assume that the polynomial family is one of the symmetric Jacobi (Gegenbauer) families scaled so that: supk kQk k! = < 1: (17) Q inf k kQk k! We denote by AN the coecient matrix of this system and assume that FN is computed by interpolation of f at the Gauss or Gauss-Lobatto points associated with the truncated orthogonal system. In [6] a post-conditioning scheme based on the integration operators is analyzed. The main di erence here is the inclusion of the T -conditions in the matrix. As the T conditions typically involve point evaluations of functions and their derivatives, we cannot expect TN to be bounded in l2 . Therefore, we introduce the space hr of in nite vectors satisfying: (18)
kZ k2hr =
1 X l=0
jZl j2(l + 1)2r < 1:
For nite truncations, the norm is simply de ned by the truncated sum. It is associated with the inner product (Y; Z)hr = Y T Dr2 Z, where Dr is the diagonal matrix whose jth diagonal entry is j r . Given any matrix C we have: (19) kC khr = kDr CDr 1 k2: Note that for integer r 0, and Z the expansion coecients of a function, z, we have, for positive constants G0,
Integration Preconditioners
7
G1, G 1: G0hz; (L + 1)r z i kZ k2hr G1hz; (L + 1)r z i (20) G 1kz k2Hr;! ; where L = DpD. (Here, h; i is the unweighted L2 inner product.) We now state our assumption on the T conditions: Assumption 3.1 There exists r0 0 and a constant, CT , such that, for all N , r r0 , and N + 1-vectors V , jTN V j CT kV khr : Moreover, for some integer r r0 , r n 1. We illustrate Assumption 3.1 with the standard examples of Chebyshev approximations and Dirichlet or Neumann boundary conditions. In the former case, X TN V = ( 1)l Vl ; (21) l
where l = l or l = 0, depending on the boundary. Choosing r0 > 1=2 we have, for all r r0: (22)
jTN V j
1 X l=0
(l + 1)
r
!=
1 2
2 0
kV khr :
jTN V j
1 X
(l + 1)4 2r0
Proof: We have:
1 X
kZ ZN k2hr = (25)
(k + 1)2r jZk j2
k=N +1
(N + 2)2(r
s)
1 X
(k + 1)2sjZk j2 ;
k=N +1
which implies the stated estimate. We note that by (20) we can replace the right-hand side by a multiple of the Hs;! norm of z. This inequality is, then, stronger than can be obtained for the Hr;! norm of z = zN . (See [3, Ch. 9].) However, we must generally use a larger r than Sobolev's inequality would require to bound T . The example of Neumann conditions illustrates this.
Theorem 3.1 a. For any r > 0 the operators B ln : hr ! hr , l = kB[ln] U khr gl;r kU khr l :
Choosing r0 > 5=2 we have, for all r r0: (24)
kZ ZN khr (N + 2)r s kZ khs :
[ ]
l
1 2
in nite vector obtained by setting all but the rst N + 1 components of Z to 0, we have:
1; : : :; n, are compact and, for some constants, gl;r ,
For Neumann conditions, X TN V = ( 1)l l2Vl : (23)
!=
Lemma 3.2 For any r 0, s > r, Z 2 hs and ZN the
kV khr :
b. For any r > 0, kSj khr , j = 0; : : :; n 1 are uniformly bounded in N .
Proof: The matrices B[ln] and Sj are banded. Bounds on
their l2 norms are given in [6] and, by Lemma 3.1, extend to their hr norms. As the compactness is proved by apWe now state a number of results concerning the indi- proximating B[ln] by its nite truncations, it can also be vidual operators in (16). In many cases the proofs can be found in [6, Sec. 4] or constructed in obvious analogy with extended. We also have an estimate for the entries of B[ln] : proofs given there, We will avoid repeating the details of j(B[ln] )km j k l : the arguments in [6]. The primary additional facts we will (26) use are: Therefore, Lemma 3.1 For any r > 0 and integer k there exists a 1 constant K(k; r) independent of N such that, for any N X kB[ln] U k2hr = (k + 1)2r j(B[ln] U)k j2 N matrix, C , with bandwidth k, l=0
kC khr K(k; r)kC k2: Proof: We need only consider the matrix Dr CDr 1 . Its r r
nonzero elements are multiplied by (i + 1) (j + 1) with i j k. Clearly this factor is uniformly bounded above by a function of k and r, completing the proof.
k=n
(27)
gl;r
1 X
(k + 1)2(r l) jUk j2;
k=0
from which the desired estimate follows. We are now in a position to prove:
ICOSAHOM 95
8
Theorem 3.2 Suppose that w = 0 is the only solution of the homogeneous system, Lw = 0, T w = 0 and that T satis es Assumption 3.1. Let r be an integer satisfying r r0 and r n 1. Further suppose that, for some p 1, f 2 C p ((a; b)). Then: a. There exist constants, C0 and C1, and an integer, N0 such that for any N > N0 and vector, y, with kykr = 1:
C0 kAN ykr C1 : b. The matrix AN I approaches a compact operator on hr . c. There exists a constant, C , and an integer N0, such that, for all N > N0 : h(u vN ); (L + 1)r (u vN )i CN 2 kf k2!;p ; = p (1=2) max(0; r n): Proof: The upper bound on AN and the compactness of AN I follow directly from Theorem 3.1. The lower bound follows from the analysis in [6, Sec. 4], which we outline here for completeness. First, de ne the compact operator, K~ : Phr ! hr , in thePfollowing way. Given Z 2 hr let z = k Zk Qk , sj z = k Zk(j )Qk . Then ~ i= (28) (KZ)
(
(TPz)i Zi ; i = 0; : : :; n 1 ( jn=01 B[nn] j Z(j ) )i ; i n:
~ = 0, it can be easily shown If, for some Z 6= 0, (I + K)Z that a nontrivial solution of the homogeneous problem exists, violating the hypotheses of the theorem. By the Riesz~ 1 khr is bounded. We next show Schauder theory, k(I + K) ~ In particular, let that AN I approximates K. ((A I)Z ) ; i = 0; : : :; N; N i (29) (K N Z)i = 0; N i > N: Here, ZN is the N + 1-vector containing the rst N + 1 components of Z. Let > 0 be given. Given any vector, Z, kZ khr = 1, and positive integer M, let ZM now denote the in nite vector obtained by setting all but the rst M + 1 components of Z to zero. Now, for M = M() suciently large, we have, for all N, ~ (30) kK(Z ZM )khr < 2 ; kK N (Z ZM )khr < 2 : Moreover, if N > M + n + q, where q is the maximum ~ M = K N ZM . Therefore, degree of the polynomials, sj , KZ for N > M() + n + q, (31) kK~ K N khr < :
Choosing suciently small, and, hence, N0 suciently large, we conclude using the Banach lemma that (AN ) 1 is uniformly bounded for N > N0 . Standard ode theory implies the existence of a solution, u 2 C p+n ((a; b)). Set e = u vN , EN = UN VN , F = FN FN and uN the polynomial whose expansion coecients are given by UN . We then have: (32) AN EN = RN ; TN UN Tu RN = B n (F) + WN ; (33) [n];N (34) WN =
nX1 j =0
(B[nn];Nj Sj;N UN (B[nn] j Sj U)N ):
From our bounds on AN and AN1 and Assumption 3.1 we conclude, (35) kEN khr C1 kRN khr ;
0
(36) kRN khr C kB[nn];N (F)khr + kU UN khr ; where in U UN , UN denotes the in nite vector obtained by extending the nite vector by 0. By Lemma 3.2 and (20) we have: (37) kU UN khr G2N r p n kf kHp;! : Using (20) and the properties of the integration operators we obtain: kB[nn];N (F)khr G4 kfN fkHl;! ; (38) (39) l = max(0; r n) 1: By the results of Bernardi and Maday on interpolation error, [2], we have: (40) kfN fkHl;! N (l=2) p kf kHp;! : Therefore, since (41) kU VN khr kEN khr + kU UN khr ; combining these inequalities and applying (20) yields the error estimate. This completes the proof. We note that these results fall short of those proved in [6] for the post-conditioning scheme, both in terms of the restrictive assumptions on the coecients and in the convergence rates. We hope to improve these in future work. Numerical experiments show that the results on conditioning hold for a number of important operators with variable lead coecients. In the following tables we display the condition numbers in the norms h0, h1 and h2, of truncated approximations to the Dirichlet problem for the cylindrical Poisson, cylindrical Helmholtz and helical Poisson operators:
Integration Preconditioners Truncation 8 16 32 64 128 256
Poisson 12.6 20.6 31.9 51.2 90.0 156.
9 Helmholtz 1264 3788 9550 19580 39320 78700
Helical 229 300 399 537 737 1020
Table 1: h0 Condition Numbers: Preconditioned Truncation 8 16 32 64 128 256
Poisson 11.4 11.7 11.9 12.0 12.1 12.1
Helmholtz 1321 2777 3844 4084 4133 4154
Helical 1110 1108 1108 1108 1108 1108
Table 2: h1 Condition Numbers: Preconditioned Poisson Operator: @ 2 + r @ k2 ; r2 @r 2 @r Helmholtz Operator:
@ 2 + r @ k 2 ; r2 r2 @r 2 @r
Truncation 8 16 32 64 128 256
Poisson 3:774 103 5:533 104 8:427 105 1:314 107 2:075 108 3:298 109
Helmholtz 9:584 59:98 922:4 14; 400 2:273 105 3:613 106
Helical 5:886 104 8:948 105 1:374 107 2:148 108 3:393 109 5:393 1010
Table 4: h0 Condition Numbers: Unpreconditioned Truncation 8 16 32 64 128 256
Poisson 4:293 103 4:978 104 6:480 105 9:225 106 1:387 108 2:150 109
Helmholtz 80:81 66:20 867:4 12; 400 1:865 105 2:891 106
Helical 9:098 104 1:102 106 1:450 107 2:070 108 3:114 109 4:827 1010
Table 5: h1 Condition Numbers: Unpreconditioned Clearly, the condition numbers grow with N in the h0 norm but remain bounded in h1 and h2. This is consistent with the analysis above. Finally, we display the growth of the condition numbers in these norms for the unpreconditioned systems, demonstrating dramatic e ects of the preconditioning.
4 General Variable Coecients
A well-known diculty with spectral methods is that multiplication by arbitrary functions is represented by full ma2 @ Therefore, direct solution of the linear systems fol2 3 2 2 2 4 4 2 4 2 @ ( r +r ) @r2 +( r +r) @r k (1+2 r + r ): trices. lowing from the Galerkin approximation to general di erential equations may be expensive. Similar considerations apply to the solution of nonlinear equations. On the other Here, k = 3, = :001, = 1:5, and 1 r 3. Helical Operator:
Truncation 8 16 32 64 128 256
Poisson 47.9 47.9 47.9 47.9 47.9 47.9
Helmholtz 3,568 7,678 10,070 10,570 10,610 10,620
Helical 10,870 10,700 10,700 10,700 10,700 10,700
Table 3: h2 Condition Numbers: Preconditioned
Truncation 8 16 32 64 128 256
Poisson 6:409 103 6:469 104 7:802 105 1:068 107 1:575 108 2:417 109
Helmholtz 2; 269 906:4 3; 887 53; 750 7:931 105 1:217 107
Helical 2:335 105 2:362 106 2:885 107 3:961 108 5:844 109 8:968 1010
Table 6: h2 Condition Numbers: Unpreconditioned
ICOSAHOM 95
10 hand, multiplication in point space may be accomplished in operations, where N is the number of points where O(N) the product is required. This fact is exploited by pseudospectral methods. In this section we show how to combine the integration preconditioners with pseudospectral approximations to multiplication by smooth functions to iteratively compute approximate solutions to variable coecient and nonlinear equations. For families with a fast interpolation algorithm, such as the Chebyshev family, the complexity of the algorithm will be O(N lnN), where N is spectral truncation order. We consider:
(46)
1 0 nX @I n + Bnn j Cj;N A v^ = Bnn F^N : 1
[ ]
j =0
[ ]
[ ]
Let AN denote the coecient matrix of the system above. Note that its rst n rows contain approximations to the constraints and its nal N + 1 n rows contain the approximation to the di erential equation. Using Lemma 4.1, we can prove the following result on conditioning and convergence. As the proof is essentially identical to the proof of Theorem 3.2, we omit it here. Theorem 4.1 Suppose that w = 0 is the only solution of 1 0 the homogeneous system, Lw = 0, T w = 0 and that T nX1 satis es Assumption 3.1. Let r be an integer satisfying (42) Lu @Dn + Dj cj (x)A u = f; x 2 (a; b); r r and r n 1 and suppose that Assumption 4.1 0 j =0 holds for this choice of r. Further suppose that, for some p 1, f 2 C p ((a; b)). Then: subject to the constraints, a. There exist constants, C0 and C1, and an integer, N0 (43) T u = d: such that for any N > N0 and vector, y, with kykr = 1 1 : Here, the functions cj are assumed to be smooth (C for C0 kAN ykr C1 : convenience). As before, we approximate u by a nite expansion, b. The matrix AN I approaches a compact operator on N X hr . (44) u v = v^i Qi (x): i=0 c. There exists a constant, C , and an integer N0, such Multiplication by cj (x) is approximated, in spectral space, that, for all N > N0 : by the following recipe: rst, evaluate the expansion at Second, mulsome interpolation points, xk , k = 0; : : :; N. h(u vN ); (L + 1)r (u vN )iC N 2 kf k2!;p ; tiply at the interpolation points by cj (xk ). Finally, use the = p (1=2) max(0; r n): new data at the interpolation points to construct expansion coecients. We denote by Cj;N the matrix representing this process. Note that Cj;N is usually never formed. 4.1 Solution by Iteration For our examples, its action is computed by fast trans operations. We may choose N > N to Although Theorem 4.1 establishes the good conditioning forms in O(N ln N) avoid aliasing errors, but always N N for some xed of the discretization matrices and the rapid convergence
as N ! 1. Interpolation points will be at Gauss or of the approximations for smooth3f, the matrix AN is full Gauss-Lobatto points associated with the family. We as- so its factorization requires O(N ) operations. However, sume the following result on the uniform boundedness of if a fast transform is available, multiplication by AN can be carried out in O(N lnN) operations. In this section we the matrices, Cj;N : exploit this feature along with the conclusions of Theorem Assumption 4.1 There exists a constant, G, and an in- 4.1 to develop an ecient iterative solution algorithm. teger r > r0, r n 1, such that: Here we consider Broyden's method, due to its ease of implementation for both linear and nonlinear problems and sup kCj;N kr G: to the availability of convergence results which are directly j;N applicable to our problem [15]. For completeness we list We expect that this assumption can be proven under ap- the algorithm as we use it, which involves only storage of propriate assumptions on the functions cj and the choice and computation with a small number of vectors of dimenof nodes. sion N + 1 [16]: Our nal speci cation of the discrete system is: Broyden's Method for the Linear System, AN v^ = Y^N : 1. Initialize: v^1 = s1 = r0 = Y^N , 1 = hs1 ; s1i. (45) TN v^ = d;
Integration Preconditioners
p
2. Until hrk ; rk i < do: rk = Y^N AN v^k ; zk(1)+1 = rk ; For j = 2; : : :; k do: ) zk(j+1 = (I + j 11sj sj 1)zk(j+11);
k+1 = hsk ; zk(k+1) i: Choose k+1 2 (0; 2) such that k k+1 k+1 6= 0;
k z (k) ; sk+1 = k k+1 k+1 k+1 1
k+1 = k+1 hsk+1 ; sk+1 i; v^k+1 = v^k + sk+1 : Here, h; i denotes some inner product and is the outer product of vectors de ned by the inner product. We choose k+1 = 1 unless k k+1 is small. Note that if p iterations are needed the total work is O(pN ln N + p2 N). Hwang and Kelley [15] have shown that if an operator, A, is such that A I is compact, then Broyden's method as described above produces a q-superlinearly convergent sequence of iterates. By part b of Theorem 4.1, this applies (uniformly in N) to our operators AN if we use the hr inner product with r r0 . Therefore, for any > 0, the number of iterations required to produce a residual with hr norm smaller than is bounded independent of N. Hence, the system can be solved in O(N ln N) operations. To illustrate this result we use Chebyshev expansions to solve: (47) D2 u + sin x u = f(x); x 2 ( 1; 1): The function f and the Dirichlet boundary conditions are chosen so that, u = p1 e (x x0 )2 = ; (48) is an exact solution. No dealiasing was used. The results tabulated are for = 5 10 4 and x0 = 1=2 and, for N 128, = 10 14. (The solution for N = 64 was so large that an absolute residual of 10 14 was unattainable. In that case only we use = 5 10 14.) Clearly, the number of iterations is independent of N while the error rapidly decreases. We note that the results presented are for the l2 inner product. The convergence does not follow directly from the theory discussed above, because the T -conditions are unbounded in this norm. However, due to their low rank, this did not harm the convergence. Tests with the h1 inner product show similar behavior. We expect that the properties of AN will lead to rapid convergence of other iterative schemes. For GMRES, this follows from [17].
11 N No. of Its. 64 33 128 23 256 23 512 23
Max. Error 1:1 103 1:2 2:5 10 5 1:0 10 11
Table 7: Linear Test Problem
4.2 Nonlinear Problems
The method may also be generalized to solve semilinear equations. In particular we consider: (49) Dn u + F(u; Du; : : :; Dn 1u; x) = 0; x 2 (a; b); with nonlinear constraints: (50) T (u) = 0: The discrete equations are formulated in spectral space by approximating F via the same recipe as above. That is, ^ v ) is computed by rst evaluating v and its derivatives in F(^ point space, then evaluating F at these points, and nally ^ Similarly, TN interpolating the point values to obtain F. is computed by evaluating T (v). The discrete system after preconditioning is given by: (51) TN (^v ) = 0; ^ v ) = 0: (52) I[n] v^ + B[nn] F(^ Broyden's method may be applied to the nonlinear discrete problem by simply de ning the residuals, rk , using equations (51-52). Generally, a good initial approximation, v^0, is needed. Linearizing about a smooth solution, the discrete system has the same properties as for the linear variable coecient equations discussed above. Therefore, the results of Hwang and Kelley [15] imply local, qsuperlinear convergence of the iterates with the number of iterations required to attain a given tolerance bounded independent of N. That is, the nonlinear discrete problem can be solved in O(N lnN) operations for a suciently good initial approximation and assuming that I is a suf ciently good approximation to the Jacobian. (Of course, di erent initial approximations to the Jacobian, which are low rank perturbations to I, could also be used.) To illustrate these results we solve the well-known reaction-di usion equation: (53) D2 u + eu = 0; x 2 ( 1; 1); u(1) = 0: For < c two solutions exist. Here, :87 < c < :88. In our example we chose = :87 and an initial guess of
ICOSAHOM 95
12
v^ = 0. The exact solution which we are approximating is consider expansions of the form given by: X f(x) = f Q ; 0 0 1 1 r !2 AA; u = ln @A @1 tanh A (54) where the expansion basis is formed as a product of (pos2x sibly di erent) m orthogonal polynomial bases and i = (i1 ; : : :; m) is a multi-index, (55) A = 2:801710482773216533343 : ::: Q (x) Qi1 ;1 (x1) Qim ;m (xm ) : We solved the problem for N = 16; 32; 64; 128;256; 512 P with = 10 14. Again, no dealiasing was employed. In We now let L = k Lk be a linear di erential operator in all cases the iterates converged in 67 to 70 iterations, con- the xi ; i = 1; : : :; m with polynomial coecients in the rming the N-independence of the iterative scheme. As independent variables. As before, rational function coefthe solution is smooth, the error was already 1:8 10 9 cients can also be allowed, provided no singularities are for N = 16 and on the order of 10 14 for ner discretiza- present in the domain and we reduce to polynomials by tions. multiplying by the least common multiple of the denominators. The Lk have the form i
i
i
i
5 Higher dimensional problems
The complexity of the spectral di erentiation operator has made the direct use of spectral methods in more than one dimension impractical. Standard treatments of commonly occuring problems, such as the Poisson [12] and Helmholtz [13] problems have been approached through diagonalization methods. These techniques perform well but have the disadvantage that an expensive matrix multiplication must be performed to transform from eigenvectors of the operators back to physical variables which is necessary, e.g. for the solution of nonlinear problems. Also, the treatment of time dependent problems is typically pursued through a nite-di erence discretization in time, which inrtoduces stability problems and limits the time-accuracy of the method. As an exeption to the latter we must mention the work of Tal-Ezer et al. [18] who employ a Chebyshev discretization of the time-evolution operator. As discussed in the monograph by Canuto et al. [3], the extension to multidimensions is open for several interesting problems in more than two space dimensions. The simplicity and generality of the integration preconditioner method can be exploited to produce block-banded forms and improve the conditioning of problems in higher dimensions treated by spectral -methods. The straightforward extension is based on the use of preconditioners constructed by tensor products of integration operators in each variable, and it allows for the use of di erent basis functions in each variable. We consider a problem in a rectangle in Rm . For simplicity we will only consider boundary conditions of Dirichlet type. We let x = (x1 ; : : :; xm) be a coordinate system such that the sides of the domain under consideration are parallel to coordinate planes. We will
Lk =
m O i=1
Lki;
with the Lki a linear di erential operator with polynomial coecients in the variable xi , i.e. an operator of the form assumed in Theorem (2.1). Let ni = max order(Lki ) : k Then the extension of theorem (2.1) to multidimensions can be stated as follows:
Theorem 5.1 The Galerkin representation of the di erential operator L in the basis Q is transformed into blockbanded form via left multiplication by the operator B=
m O i=1
B[nnii ];i ;
where 1B[1];i is the operator of integration for the family fQkigk=0 . The resullting operator has the form
BL =
m XO k i=1
B[nnii ];i Lki :
The proof follows from repeated application of theorem (2.1). We illustrate the use of theorem (5.1) by some simple examples. We limit the discussion to two space dimensions or one space dimension plus time, as we are basically interested in demonstrating the tensor product technique. Questions of conditioning and ecient implementation by the use of sparse matrix solvers will be pursued elsewhere.
Integration Preconditioners
13
In the following examples we use exclusively Chebyshev polynomial expansions, again for simplicity of exposition. Other bases could have been employed in principle, and the only added complication would have been the loss of the fast cosine transform. Thus, the preconditioners employed will be tensor products based on powers of the Chebyshev integration operator 20 0 0 0 ::: 0 3 1 66 1 0 2 0 ::: 0 7 66 0 41 0 14 0 : : : 0 777 6 .. .. 77 ... ... ... . 7 : B[1] = 666 . . . 1 .. 777 1 66 .. 0 2i 2i 64 .. . . . . . . ... 75 . 1 0 0 ::: 2M Here B[1] is a M M matrix.
Truncation 88 16 16 32 32
Abs error 1:2 10 2 4:6 10 6 1:0 10 14
Table 8: Wave equation - Exact solution: e(2(x t))2 with one redundant condition at the point ( 1; 1). Table 5 lists the absolute error for the uni-direction wave equation for various truncations. In this, as in the two subsequent examples, a homogeneous solution is chosen, and the boundary conditions are constructed by evaluating that function at the appropriate boundaries.
Example 2: the Laplace equation in a rectangle
Consider now the problem @ 2 u + @ 2 u = 0 ; (x; y) 2 [ 1; 1]2 ; Example 1: the uni-directional wave equation (58) @x2 @y2 We consider the problem @u + @u = 0 ; (x; t) 2 [ 1; 1]2 ; where u = u(x; y), f = f(x; y) and the boundary condi(56) tions are given as @t @x where u = u(x; t) and the boundary conditions are given (59) u(x; 1) = h (x) ; u(1; y) = v (y) : as For Laplace`s equation L1 = L1;1L1;2 and L2 = L2;1 L2;2 (57) u(x; 1) = h(x) ; u( 1; t) = v(t) : L1;1 = I,N L1;2 = Dx2 , L2;1 = Dy2 , and L2;2 = I and Here L1 = L1;1L1;2 and L2 = L2;1L2;2 assume the forms become is B[2] B[2] . L1;1 = I, L1;2 = Dt , NL2;1 = Dx , and L2;2 = I. The integrator For Laplace's equation BL is the sum of the following integrator for L is B[1] B[1] . two matrices The matrix BL for the wave equation is
0 3 I[1] 0 77 [1] [1] 2 I [1] [1] B[1] 4 0 0 777 4 .. 7 ... ... ... . 7: .. 77 I[1] I[1] B . 77 [1] 2i 2i .. 75 ... ... . I [1] 0 ::: B[1] 2M Note that the entries are (N +1)(N +1) blocks. The same will apply to the matrix operators given in both subsequent examples. The tau conditions are
2 0 0 66 I B 66 0 I 66 .. 66 . 66 ... 66 . 4 ..
0
M X
i=0 N X j =0
0 0
::: ::: :::
( 1)i uij = vj ; j = 0; : : :; N;
( 1)j uij = hi ; i = 0; : : :; M;
2 66 66 66 64
0 0
0 0
0 0 0 I[2] 46 0
(2 )(2 +2)
0
0
0
I[2] 4
0 0 I[2] 2(13) 0
I[2] i i
20 66 00 66 0 +6 66 0 66 0 40
0 0 0 I[2] 2(24) 0 . ..
0 0 I[2] 46 0
I[2] i i
2( )( +2)
0
0 0 0 I[2] 68 0 .. .
0 0 0 0 0 0 0 0 0 0 0 0 2 0 B[2] 0 0 0 0 2 0 0 B[2] 0 0 0 2 0 0 0 0 0 B[2] 2 0 0 0 0 B[2] 0 2 0 0 0 0 0 B[2]
0 0
0
0
0
0
0
The tau conditions can now be inserted: nX =N n=0
umn ( 1)n = h^ m ;
0 0
i
I[2]
i
(2 +2)(2 +4)
0
0 0 0 .. .
3 77 77 77 : 77 75
3 77 77 77 75
ICOSAHOM 95
14 Truncation 88 16 16 32 32
Abs error (k = 2) 3:2 10 5 5:8 10 14 4:3 10 14
Abs error (k = 8) 84:6 8:5 10 2 1:8 10 11
Truncation 88 16 16 32 32
Table 9: Laplace's equation - Exact solution: eky sin(kx) nX =N
(c = 5; = :5) :12 2:7 10 3 2:1 10 7
Table 10: Advection-di usion with x20 = 0:8 ; t0 = 1:05 Exact solution: exp( (x x40 (tc(tt0)t0)) )=(t t0) 21
umn = h^ +m ;
20 66 00 66 0 60 6 66 0 66 0 40
Again, the number of tau conditions exceeds the number of zero rows in BL. However, the tau conditions are not all independent, and four of them need to be discarded, corresponding to redundant speci cations at the four corners of the domain. This leaves 2(M + 1) + 2(N + 1) 4 conditions which matches the number of zero rows in BL. Table 5 lists the computed errors using this matrix formulation of Poisson's equation.
Example 3: the advection-di usion equation
3 7 7 7 7 7 7 0 7+
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 B[1] 0 0 0 0 0 0 0 B[1] 0 0 0 0 0 0 0 B[1] 0 0 0 0 0 0 0 B[1] 0 0 0 0 0 0 0 B[1] 0 0 0 0 0 0 0 B[1] 0 0 0 0 0 0 0 0
n=0 mX =M umn ( 1)m = v^n ; m=0 mX =M umn = v^n+ : m=0
2 66 66 66 64
0 0
0 0
0 0 0 I[1] 46 0
(2 )(2 +2)
0
0
0
I[1] 4
0 0 I[1] 2(13) 0
I[1] i i
0 0 0 I[1] 2(24) 0 . ..
0 0 I[1] 46 0
I[1] i i
2( )( +2)
0
0 68 0 .. .
77 75
0 0 ...
0 0 0 0
I[1] i
I[1]
i
(2 +2)(2 +4)
0
3 77 77 77 : 75
The tau conditions that must be imposed are:
Finally we consider the problem
@u + c @u = @ 2 u ; (x; t) 2 [ 1; 1]2 ; (60) @t @x @x2 where u = u(x; y), f = f(x; y) and the boundary conditions are given as u(x; 1) = h(x) ; u(1; t) = v (t) :
(61)
(c = 2; = 1:3) 5:6 10 3 3:6 10 5 1:5 10 9
nX =N
n=0 mX =M
umn ( 1)n = h^ m ; umn ( 1)m = v^n ;
m=0 mX =M
umn = v^n+ ;
For this equation L1 = L1;1 L1;2 and L2 = L2;1L2;2 bem=0 come L1;1 = I, L1;2 =N Dt , L2;1 = cDx Dx2 , and L2;2 = I and, again, there are two redundant tau conditions at the and integrator is B[1] B[2] . points x = 1 ; t = 1. The composite matrix BL is the sum of the following Table 5 gives the absolute errors computed for the three matrices 2 3 advection-di usion equation. 0 0 0 0 0 0
66 66 66 66 66 64 0
0
0 0
cB[1] 4
0 0 0
0 0
cB[1] 4
0 0 0
0 0 0
cB[1]
0 0
cB[1]
0
cB[1]
0 0
0 0
0 0
0 0
6
0 8
cB[1] 6
0
10
0 0 0 0
cB[1] 8
0 . .. 0
0 0 0 0 0
cB[1] 10
0 ...
0 0 . .. 0
77 77 77 77 77 75
6 Conclusions The methods discussed in this article are quite useful in deriving ecient, spectrally accurate algorithms for the treatment of initial-boundary value problems in simple geometries with more than one nonperiodic directions. Separation of variables e.g. for the Laplace operator, leads
Integration Preconditioners to equations, which can be easily transformed to form (8). As one may require the repeated solution of such equations high accuracy and eciency are clearly essential. The bad conditioning of spectral di erentiation operators is avoided by the integration preconditioning method, and this permits the treatment of problems at very high order of truncation that may otherwise be impractical. More complex geometries may be accessible as well: if a rational map to a rectangle is available, then the essential features of the method are preserved. Even if that is not feasible, the good conditioning of the resulting problems allows ecient iterative treatments to be applicable. We must remark, however, that the preconditioners discussed in this note, although quite general, might prove inappropriate for certain problems with singular behavior. The speci c structure of a given di erential operator might lead to simpler preconditioners and to more natural reduced forms for the system. An example is offered by the Laplace operator in disk geometry; indeed, in solving 4u = f in 0 r 2 ; 0 2, using a Fourier/Chebyshev expansion in the azimuthal and radial directions respectively (with 1 x = r 1 1), hwe are led 2to the2i equation^ for2 the n-th Fourier mode ((x + 1)D) n u^ = (x + 1) fn. The method discussed above would lead to a pentadiagonal, ill-conditioned operator. However, closer examination of the matrix elements reveals that under left-multiplication by a certain tridiagonal preconditioner [19] (see also [3]) we get a tridiagonal system which can be solved quite naturally by using techniques developed for the study of 3-term recurrence relations [9], and diculties relating to the coordinate singularity at x = 1 are easily bypassed. We note that Tuckerman [19] gives a theorem on the transformation of matrices into banded form through left multiplication by preconditioners whose form depends on certain properties of the matrix elements. As is also mentioned in [3], preconditioners that lead to banded form have not been readily available, and have had to be searched for in a case-by-case basis. The main appeal of the method presented here is its generality, achieved through the construction of the preconditioner from the basic recursions of a family, and its identi cation with integration operators. Indeed, the preconditioner depends only on the basis used and the order of the di erential operator L, not on its special explicit form, which can be quite complex. Also, if the coecients (or the solution) exhibit rapid variation over small neighborhoods, a rational coordinate mapping can be introduced to handle the situation with no substantial increase in algorithmic complexity while avoiding the need for considering very high-order truncations. In [6]
15 we employed a variant of the present method, using integration postconditioning, to eciently resolve shock-layer behavior through a low-order rational map. Naturally, as the Poisson equation in the disk suggests, problems with an underlying singularity may necessitate exploiting further properties of a given problem and special, tailor-made methods may need to be invented in place of the generalpurpose technique presented here.
Acknowledgment
Partially supported by DOE Grant DE-FG03-92ER25128; part of this work was performed at Ris National Laboratory, DK{4000 Roskilde, Denmark. The work of the second author was partially supported by NSF Grant DMS9304406. The authors would like to acknowledge many helpful conversations with J.P. Lynov.
ICOSAHOM 95
16
References
[15] D.M. Hwang and C.T. Kelley, Convergence of Broyden's method in Banach spaces, SIAM J. Optim., 2 (1992), pp. 505-532. [16] C.T. Kelley, Iterative Methods for Linear and Nonlinear Equations, SIAM, Philadelphia, 1995, to appear. [17] C.T. Kelley and Z.Q. Xue, Collective compactness and mesh independence of GMRES, (1995), in preparation. [18] H. Tal-Ezer, J.M. Carcione and D. Koslo , An accu-
[1] M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions, Dover, New York (1965). [2] C. Bernardi and Y. Maday, Polynomial interpolation results in Sobolev space, J. Comput. and Appl. Math., 43 (1992), pp. 53-82. [3] C. Canuto, M.Y. Hussaini, A. Quarteroni and T.A. Zang, Spectral Methods in Fluid Dynamics, SpringerVerlag, New York (1988). rate and ecient scheme for wave propagation in lin[4] C.W. Clenshaw, The numerical solution of linear difear viscoelastic media, Geophys. 55 (1990), pp. 1366ferential equations in Chebyshev series, Proc. Camb. 1379. Phil. Soc., 53 (1957), pp. 134-149. [19] L.S. Tuckerman, Transformations of matrices into [5] E.A. Coutsias, F.R. Hansen, T. Huld, G. Knorr and banded form, J. Comp. Phys. 84 (1989), pp. 360-376. J.P. Lynov. Spectral Methods in Numerical Plasma Simulation, Phys. Scripta 40 (1989), pp. 270-279. A Hypergeometric Recurrence [6] E.A. Coutsias, T. Hagstrom and D. Torres, An alRelations gorithm for the ecient spectral solution of linear ordinary di erential equations with rational function coecients, NASA Technical Memorandum 106762,
[7] [8] [9] [10] [11] [12] [13]
[14]
Hypergeometric functions satisfy the di erential equation 2 ICOMP-94-25, (1994). (62) z(1 z) ddzF2 + [c (a + b + 1)z] dF abF = 0 dz B. Fornberg, A review of pseudospectral methods for solving partial di erential equations, Acta Numerica where a, b, and c are parameters and F represents F(a; b; c; z). We show that the hypergeometric family (1994), pp. 203-267. S F(a + k; b k; c; z) satis es the following recurrence k D.G. Fox and I.B. Parker (1968), Chebyshev Polyno- relationships: mials in Numerical Analysis, Oxford University Press, (63) F(a; b; c; z) = F 0(a + 1; b 1; c; z) London. + F 0 (a; b; c; z) + F 0 (a 1; b + 1; c; z); W. Gautschi, Computational aspects of 3-term recurrence relations, SIAM Rev. 9 (1967), pp. 24-82. [2(1 z) 1]F(a; b; c; z) = ~F(a + 1; b 1; c; z) D. Gottlieb and S. Orszag, Numerical Analysis of (64) ~ b; c; z) + ~F(a 1; b + 1; c; z): + F(a; Spectral Methods, SIAM, Philadelphia (1977). ~ Here F 0 = dF dz and , , , ~ , , and ~, are all coecients L. Greengard, Spectral integration and two-point that depend on a,b, and c. boundary value problems, SIAM J. Numer. Anal., 28, The properties of the hypergeometric function that we (1991), 1071{1080. use are D.B. Haidvogel and T. Zang, The accurate solution of ab d Poisson's equation by expansion in Chebyshev polyno- (65) dz F(a; b; c; z) = c F(a + 1; b + 1; c + 1; z) mials, J. Comp. Phys. 30 (1979), pp. 167-180. (b a)(1 z)F(a; b; c; z) (c a)F(a 1; b; c; z) P. Haldenwang, G. Labrosse, S. Abboudi and M. Deville, Chebyshev 3-D spectral and 2-D pseudospectral (66) +(c b)F(a; b 1; c; z) = 0 ; solvers for the Helmholtz equation, J. Comp. Phys. 55 and the Gauss contiguity relations (1984), pp. 115-128. D. Hochstadt, Special Functions of Mathematical (67) (c a 1)F(a; b; c; z) + aF (a + 1; b; c; z) (c 1)F(a; b; c 1; z) = 0 ; Physics, Dover, NY (1975).
Integration Preconditioners (68)
(b a)F(a; b; c; z) + aF(a + 1; b; c; z) bF(a; b + 1; c; z) = 0 :
17 (70) ( 1)( ; ; z) + ( + 1; ; z) ( 1)( ; 1; z) = 0 :
The following equations are required in the derivation of Using Equation (70) evaluated at ( + 1; + 1), ( + 1; ), equation (63): Equation (68) evaluated the points (a; b; c), and ( ; ), and Equation (69), one can derive (a+1; b; c), and (a; b+1; c); Equation (65) evaluated at the points (a+1; b 1; c 1), (a; b; c 1), and (a 1; b+1; c 1); (71) ( ; ; z) = 0 ( ; ; z) + 0( 1; ; z): 1 and Equation (67). After involved algebra, one obtains the rst recurrence relation (63) with A similar recurrence relation can be derived for con uent hypergeometric functions of the second kind. These functions satify the following two equations = (b 1)(ba(ca)(bb) a 1) ; d ( ; ; z) = ( + 1; + 1; z); (72) (a + b + 1) 2c dz = (b a 1)(b a + 1) ; (73) ( ; ; z) ( + 1; ; z) ( ; 1; z) = 0: b(c a) Using (73) evaluated at points ( +1; +1) and ( +1; )
= (a 1)(b a)(b a + 1) : and Equation (72), one can generate The derivation of the second recurrence relation (64) requires the following equations: Equation (68) evaluated (74) ( ; ; z) = 1 1 0 ( 1; ; z) + 0 ( ; ; z): at the points (a 1; b; c) and (a; b 1; c) and Equation (66). After some involved algebra, one can generate the second In addition to the recurrence equations developed above, recurrence relation (64) with both types of con uent hypergeometrics satisfy a recurrence of the form 2(c b)a ~ = (b a)(b a 1) ; zf( ; ; z) = c1 f( 1; ; z) + c2 f( ; ; z) + c3f( + 1; ; z): ~ = (a + b 1)(a + b + 1 2c) ; (b a + 1)(b a 1) The rst con uent hypergeometric functions satisfy a)b (75) ( )( 1; ; z) + (2 + z)( ; ; z)
~ = (b 2(c a)(b a + 1) : ( + 1; ; z) = 0:
B Con uent Hypergeometric Recurrence Relations Con uent hypergeometric functions satisfy the di erential equation 2 z ddzu2 + ( z) du dz u = 0
This equation can be rearranged into the form (76)
z( ; ; z) = ( )( 1; ; z) + ( 2 )( ; ; z) + ( + 1; ; z):
The con uent hypergeometric functions of the second kind satisfy
( 1; ; z) (2 + z) ( ; ; z) We show that the con uent hypergeometric functions (77) + ( + 1) ( + 1; ; z) = 0; satisfy recurrence relations analogous to the recurrence relations for hypergeometric functions. There are two types of con uent hypergeometric functions. Each one is treated which can be rearranged into the form separately. z ( ; ; z) = ( 1; ; z) We start with the rst con uent hypergeometric func- (78) (2 ) ( ; ; z) + ( + 1) ( + 1; ; z): tion. It satis es the following equations: d ( ; ; z) = ( + 1; + 1; z); (69) dz
ICOSAHOM 95
18
Family Chebyshev Tk Q0 1 Q1
x
ak 1;k ak;k ak+1;k bk 1;k
x k
1 2 0 1 2 1
hk
2(k 1) 0 1 2(k + 1) (1 x2 ) 1=2 (1 x2 )1=2 ( 1; 1) =2(; k = 0)
k
k2
bk;k bk+1;k w(x) p(x) (a; b)
Legendre Pk 1 2k + 1 0 k+1 2k + 1 1 2k + 1 0 1 2k + 1 1 (1 x2 ) ( 1; 1) 2 2k + 1 k(k + 1)
Gegenbauer Ck( ) 1 2x k + 2 1 2(k + ) 0 k+1 2(k + ) 1 2(k + ) 0 1 2(k + ) (1 x2 ) 12 (1 x2 ) + 12 ( 1; 1) 21 2 (k + 2 ) k!(k + )[ ( )]2 k(k + 2 )
Jacobi Pk( ; ) 1 1 (( ) + ( + + 2)x) 2 2(k + )(k + ) (2k + + +2 1)(22k + + ) ( ) (2k + + + 2)(2k + + ) 2(k + 1)(k + + + 1) (2k + + + 2)(2k + + + 1)
ak 1;k k+ + 2ak;k + ak+1;k k+1 (1 x) (1 + x) (1 x2 )w(x) ( 1; 1) 2 + +1 (k + + 1) (k + + 1) (2k + + + 1)k! (k + + + 1) k(k + + + 1)
Hypergeometric F (a + k; b k; c) 2(c a k)(b k) (b a 2k)(b a 2k + 1) (a + b 1)(a + b + 1 2c) (b a 2k + 1)(b a 2k 1) 2(c b + k)(a + k) (b a 2k)(b a 2k 1) ak 1;k k+a 1 2ak;k a+b 1 ak+1;k k+1 b (1 + x)c 1 (1 x)a+b c (1 + x)c (1 x)a+b c+1 ( 1; 1)
(a + k)(b k)
Table 11: Recursions for polynomial families of Hypergeometric type (a0;k = 0; for the Tk , a1;0 = b1;0 = 1)
Family Hermite Hk Q0 1 Q1 2x ak 1;k ak;k ak+1;k bk 1;k bk;k bk+1;k w(x) p(x) (a; b) hk k
k
0 1 2 0 0 1 2(k + 1)
e x2 e x2 ( 1; 1) p 2k k!
2k
Laguerre L(k ) 1 1+ x (k + ) 2k + + 1 (k + 1) 0 1 1
Con uent ( rst) (a; b) a b+k b 2a 2k a+k 1 0
Con uent (second) (a; b) 1 b 2a 2k (a + k)(a b + k + 1) 1 1 a k 1 0
x e x x +1 e x (0; 1) ( + k + 1) k! k
xb 1 e x xb e x (0; 1)
xb 1 e x xbe x (0; 1)
(a + k )
(a + k )
b a k a+k 1
Table 12: Recursions for Con uent Hypergeometric functions
