- Email: [email protected]
A convergence analysis for the superconsistent Chebyshev method
Applied Numerical Mathematics 58 (2008) 88–100 www.elsevier.com/locate/apnum
A convergence analysis for the superconsistent Chebyshev method Lorella Fatone a,∗ , Daniele Funaro a , Gang Joon Yoon b a Dipartimento di Matematica Pura e Applicata, Università di Modena e Reggio Emilia, Via Campi 213/b, 41100 Modena, Italy b School of Mathematics, Korea Institute for Advanced Study, 207-43 Cheongnyangni 2-dong,
Dongdaemun-gu, Seoul 130-722, Republic of Korea Available online 12 December 2006
Abstract The superconsistent collocation method is based on collocation nodes which are different from those used to represent the solution. The two grids are chosen in such a way that the continuous and the discrete operators coincide on a space as larger as possible (superconsistency). There are many documented situations in which this technique provides excellent numerical results. Unfortunately very little theory has been developed. Here, a theoretical convergence analysis for the superconsistent discretization of the second derivative operator, when the representation grid is the set of Chebyshev Gauss–Lobatto nodes is carried out. To this end, a suitable quadrature formula is introduced and studied. © 2006 IMACS. Published by Elsevier B.V. All rights reserved. PACS: 02.60.Lj; 02.70.Jn Keywords: Convergence; Chebyshev polynomials; Consistency; Collocation grids
1. Introduction An incredible variety of numerical schemes, for the solution of a wide range of problems, can be generated by using two grids at the same time. The first grid is used to describe the approximated solution. As a matter of fact, the computed values are directly related to it. Through a sort of duality process, the second grid is used instead to enforce the equations to be satisfied. There, the residual is required to vanish. This approach allows for the creation of new schemes or the improvement of existing ones. Actually, we may interpret most of the well-known schemes as methods in which the two grids coincide. In order to appreciate the performances of this technique, the two grids must be somehow related. The first one is in general chosen for convenience (because we want the solution to be expressed in that way). The second one should be taken with the aim of obtaining the best from the approximation scheme. A very effective way of choosing the second grid is to build it by asking the discrete operator to be as much “accurate” as possible. For example, when the first grid is fixed, the superconsistent method indicates how to obtain the second grid, by requiring the discretized operator to coincide with the exact operator in a finite dimensional space as larger as possible.
* Corresponding author.
E-mail addresses: [email protected] (L. Fatone), [email protected] (D. Funaro), [email protected] (G.J. Yoon). 0168-9274/$30.00 © 2006 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2006.11.001
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
89
The applications of this idea are unlimited. For the Legendre spectral collocation method in two dimensions we address the reader to [3], where the superiority of the superconsistent method is shown, mainly experimentally, especially for transport-dominated advection–diffusion equations. A survey of possible other applications is provided in [5]. Superconsistent collocation techniques are also considered for integral type equations in [6]. Usually, the extra costs of computing the residual grid and implementing the corresponding discrete operator, are well compensated by the advantages (in term of accuracy and computer time) that one can achieve. So far, except for the results in [7] (dealing with finite-differences for pure hyperbolic problems), very little theory has been developed, despite the large amount of successful experimental material available. This is also true for the simplest examples. This paper, which is quite technical, is aimed to finally provide some theoretical proof of convergence. Some of these results were anticipated in [8]. In that paper Legendre and Chebyshev collocation methods on special (superconsistent) grids were considered for a simple one-dimensional second-order equation. The main interest was to study the properties of suitable preconditioners for the spectral discretization matrices. The technique was successively applied to one-dimensional transport-dominated equations, showing also in that situation the reliability of the method and the effectiveness of the new preconditioners. In order to study the symmetric Chebyshev case, a suitable quadrature formula was introduced. It was noted that such a formula is related to an interesting set of orthogonal polynomials of Sobolev type and conjectures were made about its properties. As a matter of fact, the knowledge of the behavior of nodes and weights of the formula is of crucial relevance for the development of a convergence theory for the approximation of the differential equations examined. The purpose of this paper is to prove these conjectures and add further insight. This will be used to carry out a rigorous convergence analysis for the Chebyshev superconsistent collocation method. Unfortunately, the theory here presented is able to cover only a very simple case, which is marginal with respect to the potentiality of the superconsistent method. In this circumstance, the improvements obtainable by adopting the two-grids approach are irrelevant if compared with the standard Chebyshev collocation method (which can be seen as another two-grids method, having coincident grids). However, together with the results developed for the quadrature formula (that certainly have a general interest for other types of applications), there is now the hope that this preliminary analysis, carried out on a very basic problem, could open the path to more serious generalizations. d2 In order to start, let us briefly recall some general assumptions. Let L be the second derivative operator: L = − dx 2. Then consider the following Dirichlet boundary-value problem: Lu = −u = f
in ]−1, 1[,
with u(−1) = u(1) = 0,
(1)
where f is a given function. The standard Chebyshev collocation method to approximate u is constructed by taking the points of the grid x to be the zeros of Tn . Then, one looks for un ∈ P0n such that: −un (xi ) = f (xi )
for i = 1, . . . , n − 1,
(2)
is the space of polynomials of degree less than or equal to n, vanishing at x = ±1. where In the superconsistent method (see [4]), one introduces another collocation grid z and looks for a new solution vn ∈ P0n in such a way that: P0n
−vn (zi ) = f (zi )
for i = 1, . . . , n − 1.
(3)
Each node zi , i = 1, . . . , n − 1, of the new grid z is required to be solution to the equation: (zi ) = 0, (Lχn+1 )(zi ) = −χn+1
i = 1, . . . , n − 1,
(4)
where
χn+1 = 1 − x 2 Tn ∈ P0n+1 .
(5)
Actually, it is easy to realize that Eq. (4) has exactly n − 1 distinct zeros. This choice is made for the following reasons. If one expands vn in the usual Lagrange basis related to the grid x ∪ {±1}, a (n − 1) × (n − 1) discretization matrix An−1 corresponding to (3) is obtained. Thus, we observe that both the exact operator L and the discrete operator, when applied to χn+1 , are vanishing at the nodes z. The first property is true due to (4), the second one is a consequence of the fact that the matrix An−1 only acts on the values attained at the nodes x , that in the case of the function χn+1 are zero. This means that the kernel of the difference between
90
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
the exact and discrete operators, not only contains the space P0n , but also the enlarged space P0n+1 . This justifies the name superconsistent. Note that we only enlarged the approximation space of one dimension, but this seems to be sufficient for achieving good results in applications. As documented by experiments in [8], in terms of accuracy, the superconsistent method turns out to be superior to the standard collocation method, so that the error |u − vn | is generally less than the error |u − un |. Actually, the Lebesgue constant associated to the distribution of nodes z ∪ {±1} is less than the one corresponding to the nodes x ∪ {±1}. The convergence, for n → +∞, of the discrete solution un , obtained by collocation at the Chebyshev nodes, to the exact solution u of problem (1), is proven in [1], where estimates for the spectral decay of the error are given. Such a proof, based on variational techniques, makes use of the exactness of Gaussian formulas up to degree 2n − 1. In order to analyze the convergence also for the discrete solution vn , it is necessary to introduce and study other suitable high-order integration formulas based on the grid z. This is what will be done in the coming sections. 2. Study of a quadrature formula of Sobolev type are even (odd) functions. Let us start by observing that, if n is odd (even), then χn+1 and χn+1 It should be noted that the grid z is different from the grid x . As a matter of fact, in [8] the following estimates regarding the location of the nodes are given: ⎧ ⎨ xj < zj < ξj +1 , j = 1, . . . , [n/2], j = [n/2] + 1, . . . , n − 1, ξj < zj < xj , (6) ⎩ z = 0, n even and j = n/2, j −1)π , j = 1, . . . , n. However, it is where [∗] is the integer part of ∗ and ξ are the zeros of Tn , i.e.: ξj = − cos (2j 2n interesting to point out that there may be cases in which the grid x and its corresponding superconsistent grid coincide. In fact, in [8], the following theorem was proven. 2
d Theorem 1. Concerning the differential operator L = − dx 2 , with homogeneous Dirichlet boundary conditions, one has that z = x if and only if the points xi , i = 1, . . . , n − 1, are the zeros of the derivative Pn of the Legendre polynomial Pn .
With the help of some trigonometry, another characterization of the grid z can be given. Assuming that, for i = 1, . . . , n − 1, each point of z is written as zi = cos θi , for a certain angle θi , then the following lemma holds. Lemma 2. For each node zi = cos θi , i = 1, . . . , n − 1, the following equation holds: tan nθi =
n cos θi sin θi n2 sin2 θi + 1
,
i = 1, . . . , n − 1.
(7)
Proof. The Chebyshev polynomial Tn (x) = cos nθ , with x = cos θ , satisfies the Sturm–Liouville problem: 1 − x 2 Tn − xTn + n2 Tn = 0. The differentiation of Eq. (8) gives the two relations: 1 − x 2 Tn = − n2 Tn + xTn , 1 − x 2 Tn = − n2 + 1 Tn + xTn . Therefore, from (4) we have: −χn+1 (zi ) = 0 ⇐⇒ n2 + 1 Tn (zi ) + zi Tn (zi ) = 0,
(8)
(9) (10) i = 1, . . . , n − 1.
(11)
Since: n sin nθ , x = cos θ, sin θ n sin nθ cos θ n2 cos nθ Tn (x) = − , 3 sin θ sin2 θ Tn (x) =
(12) x = cos θ,
(13)
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
each θi , i = 1, . . . , n − 1, is solution of: n n2 + 1 sin nθi sin2 θi + n sin nθi cos2 θi − n2 cos nθi cos θi sin θi = 0,
91
(14)
that is:
sin nθi n2 sin2 θi + 1 = n cos nθi cos θi sin θi ,
(15)
2
from which we arrive to (7).
Due to the above considerations, one can see that the difference between two consecutive nodes behaves like n−1 in the central part of the interval [−1, 1], and like n−2 near the endpoints ±1. Then, let us recall some other lemmas, previously obtained in [8], characterizing the superconsistent Chebyshev nodes. Lemma 3. Let χn+1 = (1 − x 2 )Tn . Then, the following orthogonality relation holds: 1
χn+1 qn ω dx = 0,
−1
1 with ω(x) = √ , 1 − x2
(16)
for all polynomials qn of degree n with qn (±1) = 0 and qn (±1) = 0. This lemma provides a link between the Chebyshev superconsistent grid z ∪ {±1} and the so-called Sobolevtype orthogonal polynomials (see, for example, [9]). In fact, given a function f , the set of nodes z is related to the integration formula: 1 f ω dx ≈
n
f (zi )ωi + f (−1)ω˜ 0 + f (1)ω˜ n ,
(17)
i=0
−1
where z0 = −1, zn = 1. The weights ωi , i = 0, . . . , n, and ω˜ 0 , ω˜ n , are defined as follows: 1 ωi =
lˆi ω dx,
i = 0, . . . , n,
(18)
−1
1 ω˜ 0 =
l˜0 ω dx,
1 ω˜ n =
−1
l˜n ω dx,
(19)
−1
where the Lagrange polynomial basis used in (18), (19) is defined as: lˆi ∈ Pn+2 with lˆi (zj ) = δij l˜0 ∈ Pn+2 with l˜0 (zi ) = 0, l˜n ∈ Pn+2 with l˜n (zi ) = 0,
and lˆi (±1) = 0, i = 0, . . . , n, l˜0 (−1) = 1, i = 0, . . . , n,
l˜n (−1) = 0,
(20) l˜0 (1) = 0, l˜n (1) = 1.
(21) (22)
The n + 3 polynomials satisfying conditions (20)–(22) are unique. It is possible to rewrite them with the help of the polynomial of degree n + 1: Ω(x) = (x − z0 ) · · · (x − zn ) = cn x 2 − 1 χn+1 (x),
with cn = −
1 , 2n−1 n2 (n + 1)
(23)
and defining: i (x) =
Ω(x) , − zi )
Ω (zi )(x
i = 0, . . . , n.
(24)
92
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
Thus, it is not hard to show that: ˆl0 (x) = (x − 1) 1 0 (−1)(x + 1) − 1 (x + 3) 0 (x), 2 4 (1 − x 2 ) i (x), i = 1, . . . , n − 1, lˆi (x) = (1 − zi2 ) ˆln (x) = −(x + 1) 1 n (1)(x − 1) + 1 (x − 3) n (x), 2 4 1 l˜0 (x) = 1 − x 2 0 (x), 2 1 l˜n (x) = − 1 − x 2 n (x). 2 Furthermore, since the zeros z0 , z1 , . . . , zn are symmetric with respect to x = 0, that is: zi = −zn−i ,
i = 0, . . . , n,
we come to: lˆi (−x) = lˆn−i (x),
i = 0, . . . , n,
(25)
(26) and l˜0 (−x) = −l˜n (x).
(27)
Consequently, as expected, one has ωi = ωn−i , i = 0, . . . , n, and ω˜ 0 = −ω˜ n . In [8] the following result has been proven. Lemma 4. The quadrature formula (17) is exact if f is a polynomial of degree less than or equal to 2n − 1. Some other properties of Chebyshev polynomials will be useful in the following. In particular for Tn , it is not difficult to get the relation (see [10] or [11]): n (28) 1 − x 2 Tn = {Tn−1 − Tn+1 }. 2 Differentiating the above equation yields the relation: n n + 2xTn − Tn+1 . (29) 1 − x 2 Tn = Tn−1 2 2 On the other hand, by recalling (12), we can write Tn (x) = nUn−1 (x), with x = cos θ , where Un−1 is the Chebyshev polynomial of the second kind, defined by: 1 sin nθ , x = cos θ. (30) Un−1 (x) = Tn (x) = n sin θ Moreover, from (8) one has: n Tn = (xUn−1 − nTn ). (31) 1 − x2 By taking the second derivative in (28), we arrive to: n = {Tn−1 − Tn+1 } χn+1 2
n (n − 1)xUn−2 − (n − 1)2 Tn−1 − (n + 1)xUn + (n + 1)2 Tn+1 . (32) = 2(1 − x 2 ) In addition, the following relations will be also useful (see again [10]): Tn = Un − xUn−1 ,
(33)
U2n+1 = 2Un Tn+1 ,
(34)
xU2n−1 = 1 + 2T2 + · · · + 2T2n−2 + T2n ,
(35)
U2n = 1 + 2T2 + · · · + 2T2n .
(36)
Now the following theorem can be proven.
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
93
Theorem 5. The weights ωi , i = 0, . . . , n in (18) and ω˜ 0 in (19) are positive (as a consequence, ω˜ n is negative). , the corresponding weight ω is positive. Let Q be Proof. First of all, let us show that for any zi = 0, zero of χn+1 i i the polynomial of degree n − 3 defined by
Qi (x) =
(x) χn+1
x 2 − zi2
(37)
.
Then, let fi be the polynomial of degree 2n − 2 defined by: 2 fi (x) = 1 − x 2 Q2i (x).
(38)
It is easy to realize that fi is a non-negative even function such that fi (±1) = 0, fi (zi ) = fi (zn−i ) > 0 and fi (zk ) = 0 for k = i, n − i. The quadrature formula ensures: 1 fi ω dx = ωi fi (zi ) + ωn−i fi (zn−i ) = 2ωi fi (zi ) > 0,
(39)
−1
showing that the corresponding weight ωi must be positive. The next step is to show the positiveness of ω˜ 0 . By definition (see (19)), one has: 1 ω˜ 0 =
1
l˜0 ω dx =
(−1) 4χn+1
−1
1
1 − x 2 (1 − x)ω dx. χn+1
(40)
−1
From (10) we can write: 1 ω˜ 0 = 4χn+1 (−1)
1
2 n + 1 Tn + xTn 1 − x 2 (x − 1)ω dx.
(41)
−1
Since the derivatives of Chebyshev polynomials form a set of orthogonal functions with respect to (1 − x 2 )ω, the following orthogonality relation holds: 1
1 − x 2 ω dx = 0, Tn rn−1
∀rn−1 ∈ Pn−1 ,
(42)
−1
which leads to: ω˜ 0 =
1 (−1) 4χn+1
1
xTn 1 − x 2 (x − 1)ω dx.
(43)
−1
The Sturm–Liouville equation (8) and the orthogonality relations of Chebyshev polynomials (including (42)) bring for n > 2 to: 1 ω˜ 0 = (−1) 4χn+1 1 = 4χn+1 (−1) 1 = 4χn+1 (−1)
1
x xTn − n2 Tn (x − 1)ω dx
−1
1
x 2 Tn (x − 1)ω dx
−1
1 −1
Tn (1 − x)
1 1 − x 2 dx
+ −1
Tn (x
− 1)ω dx
94
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
1 = 4χn+1 (−1)
1
Tn (x − 1)ω dx.
(44)
−1
In order to compute ω˜ 0 explicitly, we first observe that (see [10]): Tn (1) = n2 ,
Tn (1) =
n4 − n2 , 3
and Tn (1) =
n2 (n2 − 1)(n2 − 4) . 15
(45)
Then, using (10), we get:
χn+1 (−1) = (−1)n χn+1 (1) 2(−1)n 2 2 = (−1)n n2 + 1 Tn (1) + Tn (1) = n 2n + 1 . 3 On the other hand, (34)–(36) ensure that: 1
Tn (x − 1)ω dx = (−1)n nπ.
(46)
(47)
−1
Consequently, one recovers the following value of ω˜ 0 : ω˜ 0 =
3π = −ω˜ n . 8(2n3 + n)
(48)
)2 the It is now the time to check that the remaining weights ω0 = ωn are positive. To this end, we apply to f = (χn+1 quadrature formula (17) obtaining:
1
f ω dx = 2ω0 f (1) + ω˜ 0 f (−1) − ω˜ 0 f (1).
(49)
−1
Due to the fact that f (1) > 0 and f (−1) = −f (1) < 0, considering that the integral above is positive and that ω˜ 0 > 0, we deduce the positivity of ω0 . At the same time we entail that ωn is also positive. is allowed to have zm = 0 as a root. Therefore, There is one more case left, corresponding to n = 2m, where χn+1 we have to deal with the weight ωm . When m = 1, it is easy to see that: 1 ω1 =
2 1 − x 2 ω(x) dx > 0.
(50)
−1
Now let us consider the case when m 2. Let us introduce a polynomial f of degree n + 3 = 2m + 3 as: 2 f (x) = 1 − x 2 Q(x),
(51)
where now Q is the polynomial given by: Q(x) =
(x) χ2m+1
. x By applying to f , defined in (51), the quadrature formula (17), we obtain the equation: 1 −1
2 χ2m+1 ω dx = Q(0)ωm . 1 − x2 x
(52)
(53)
Then, considering (23), it is easy to see that: Q(0) = lim Q(x) = x→0
m−1 (−1)m−1 2 zk , cn i=1
(54)
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
where the coefficient cn is provided in (23). Therefore one gets the sign of Q(0): sgn Q(0) = (−1)m .
95
(55)
On the other hand, Eq. (32) and the orthogonality of Ui with respect to (1 − x 2 )ω give: 1 −1
2 χ2m+1 ω dx = m 1 − x2 x
1
−1
1 (2m + 1)2 T2m+1 − (2m − 1)2 T2m−1 1 − x 2 2 dx x 1
= m 4m2 + 1
−1
1 T2m+1 − T2m−1 1 − x 2 2 dx, x
(56)
where we used the relation: T2m+1 + T2m−1 = 2xT2m
(57)
and that the weighted integral of this last quantity is zero. Going ahead, the relations (33), (34) and the orthogonality condition (42) bring to: 1 −1
1 T2m+1 − T2m−1 1 − x 2 2 dx = x
1
−1
1 U2m+1 − U2m−1 1 − x 2 2 dx x
1
1 Um Tm+1 − Um−1 Tm 1 − x 2 2 dx. x
=2 −1
(58)
When m = 2i is even, we have: 1 −1
1 Um Tm+1 − Um−1 Tm 1 − x 2 2 dx = x
1
−1
1 =
1 T2i+1 U2i 1 − x 2 2 dx − x
1
−1
1 2 1 − x 2 2 dx + U2i
−1
U2i−1 1 − x2 T2i ω dx x
1 T2i2 ω dx > 0.
(59)
−1
When m = 2i + 1 is odd, the same arguments bring to: 1 −1
1 Um Tm+1 − Um−1 Tm 1 − x 2 2 dx = x
1
−1
U2i+1 T2i+2 ω dx − 1 − x2 x
1
−1
U2i
1 T2i+1 1 − x 2 2 dx < 0. x
Hence we have shown that 1 χ 2 2 2m+1 ω dx = (−1)m , 1−x sgn x
(60)
(61)
−1
which, from (53), implies ωm > 0. This completes the proof.
2
In addition, the following results can be proven. Lemma 6. For each n 2, one has the two inequalities: n 1 1
2 2 C 2 2 φ ω dx φ (zi )ωi + φ (−1)ω˜ 0 + φ (1)ω˜ n 2 φ 2 ω dx, n −1
i=0
−1
for all φ polynomials of degree n, where the positive constant C does not depend on n.
(62)
96
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
Proof. First, let us show the second inequality. Since φ can be written as a multiple of Tn plus lower degree terms, due to the exactness of the quadrature formula up to degree 2n − 1, it is enough to prove the inequality for φ = Tn . Therefore, it is necessary to show that: n
Tn2 (zi )ωi
+ 2Tn (−1)Tn (−1)ω˜ 0
+ 2Tn (1)Tn (1)ω˜ n
i=0
1 B
Tn2 ω dx,
(63)
−1
where B does not depend on n. Since Tn (x)Tn (x) is always an odd function and ω˜ n = −ω˜ 0 , the part containing the boundary terms on the left-hand side of (63) is less than zero: 2Tn (−1)Tn (−1)ω˜ 0 + 2Tn (1)Tn (1)ω˜ n = 4Tn (−1)Tn (−1)ω˜ 0 < 0.
(64)
Thus, we just have to deal with the summation. Let us note that, by taking f = 1 in (17) , we have: 1 ω dx = π =
n
(65)
ωi .
i=0
−1
Moreover, since |Tn (x)| 1, ∀x ∈ [−1, 1] and ωi > 0, i = 0, . . . , n, the sum in (63) can be easily bounded. As a matter of fact: n
Tn2 (zi )ωi
i=0
n
i=0
π ωi = π = 2 · = 2 2
1 Tn2 ω dx.
(66)
−1
Therefore (63) is true by choosing B = 2. In order to show the first inequality in (62), due to the exactness of the quadrature formula (17) up to degree 2n − 1, )2 (1 − x 2 ) ∈ P . Therefore, it is necessary to have: we now choose φ 2 = (χn+1 2n C n
1
2 (χn+1 )2 1 − x 2 ω dx 4 χn+1 (1) ω˜ 0 .
(67)
−1
Such an inequality is correct since Eqs. (46)–(48) ensure that the right-hand side of (67) behaves like n5 . On the other hand, the left-hand side of (67) behaves like n6 . In fact Eq. (10) entails: 1
(χn+1 )2
1 − x 2 ω dx =
−1
1
2 2 n + 1 (Tn )2 + (xTn )2 + 2x n2 + 1 Tn Tn 1 − x 2 ω dx,
(68)
−1
and by taking only the first integral in (68) one has: 1
2 2 2 n + 1 (Tn )2 1 − x 2 ω dx = n2 n2 + 1
−1
that grows as n6 .
π
2 π , sin2 (nθ ) dθ = n2 n2 + 1 2
(69)
0
2
We continue by collecting other results about our quadrature formula. In particular the asymptotic behavior of the weights ω0 (= ωn ) and ω˜ 0 (= −ω˜ n ) is investigated. Let us start by defining a function λ on positive integers: 1, if n is even, λ(n) = (70) 0, otherwise. Lemma 7. For the Chebyshev polynomials Tn , the following relations can be proven:
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
1
Tn ω dx =
−1
1
n3 πλ(n), 2
xTn ω dx =
−1
1
(71)
n(n2 − 1) πλ(n − 1), 2
(72)
n(n2 − 2) πλ(n). 2
x 2 Tn ω dx =
−1
97
(73)
Proof. First of all, we show (72). For n = 2k + 1, the identities (35) and (36) imply that 1
xTn ω dx
1 =n
−1
xU2k ω dx
1
1 k k
= 2n x T2i ω dx = 2n 2i xU2i−1 ω dx = 2nk(k + 1)π.
−1
−1
i=1
i=1
(74)
−1
= (xU2k−1 ) − U2k−1 , the same arguments Since k = 12 (n − 1), we have obtained (72). Using the identity xU2k−1 show (73). Finally, the identity (71) follows directly from (73) and:
1
1 − x 2 Tn ω dx = n
−1
1 (xUn−1 − nTn )ω dx = nπλ(n),
(75)
−1
2
where Eq. (31) is used.
The asymptotic behavior of the two weights ω0 (= ωn ) and ω˜ 0 (= −ω˜ n ) is established by the following lemma. Lemma 8. For n large, up to multiplicative constants, one has: ω0 ∼ n−1
and ω˜ 0 ∼ n−3 .
(76)
Proof. Eq. (48) ensures that ω˜ 0 ∼ n−3 . Now the estimation of the behavior of ω0 is required. Using the identity (32) and the expression of the polynomial n in (24): n =
(x + 1)χn+1 (1) 2χn+1
(77)
,
it is possible to rewrite lˆn as: lˆn = − n (1) + 1/2 l˜n +
n 2 (1) (x + 1) (Tn−1 − Tn+1 ). 8χn+1
(78)
Hence, the weight ω0 is given by: ω0 = −
n (1) +
1 2 1 n ω˜ n + − Tn+1 )ω dx. x + 2x + 1 (Tn−1 2 8χn+1 (1)
(79)
−1
On the other hand, using the expression for χn+1 in (10), the properties of Chebyshev polynomials given in (45), and Lemma 7, for n large, up to multiplicative constants, one has:
n (1)ω˜ n
∼n
−1
,
n (1) 8χn+1
1
2 x + 2x + 1 (Tn−1 − Tn+1 )ω dx ∼ n−1 .
(80)
−1
It turns out that the asymptotic behavior of ω0 is n−1 . This completes the proof.
2
98
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
3. Convergence analysis The convergence of the discrete solution vn of (3) (obtained by collocation at the superconsistent Chebyshev grid) to the solution u of problem (1) can now be proven. This result was not available in [8], although there it was shown experimentally that the error |u − vn | decays faster than the error |u − un |, and this last quantity is known to converge to zero exponentially. The positivity of the weights of the quadrature formula (17), resulting from Lemma 5, plays an important role in the theorem that follows. Theorem 9. The discrete superconsistent Chebyshev solution vn of problem (3) converges, for n → +∞, to the exact solution u of problem (1) and the rate of convergence is of spectral type. {z}
Proof. Let In−2 f be the interpolant of f at the nodes z, i.e.:
{z} In−2 f (zi ) = f (zi ),
i = 1, . . . , n − 1.
(81)
From (3) we can write: {z}
−vn = In−2 f
in ]−1, 1[,
(82)
{z}
since vn and In−2 have both degree n − 2 and coincide in n − 1 points. { x}
Similarly, by denoting by In−2 f the interpolant of f at the nodes x , i.e.:
{ x} In−2 f (xi ) = f (xi ),
i = 1, . . . , n − 1,
(83)
we deduce that: { x}
−un = In−2 f
in ]−1, 1[.
(84)
1 = H 1 (−1, 1) in the usual way, by the After introducing the weighted Sobolev spaces L2ω = L2ω (−1, 1) and H0,ω 0,ω triangle inequality, one gets:
u − vn H 1 u − un H 1 + un − vn H 1 . 0,ω
0,ω
0,ω
(85)
It is already known that the first term on the right-hand side of (85) converges spectrally to zero, so that it is now necessary to estimate the second one. From (82) and (84) we get: { x}
{z}
−(un − vn ) = In−2 f − In−2 f
in ]−1, 1[.
(86)
By using an inequality proven in [1], Eq. (86) gives: 1 un − vn 2H 1 C3 0,ω
(un − vn ) (un − vn )ω dx
−1
1 = −C3
(un − vn ) (un − vn )ω dx
−1
1 = C3
{x } {z} In−2 f − In−2 f (un − vn )ω dx,
(87)
−1
where C3 is a constant not depending on n. The integration by parts takes into account that (un − vn )(±1) = 0. At this point, it is possible to continue with two alternative proofs. The first one only works for n even. Thanks to the Schwarz inequality, we can write:
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
1 un − vn 2H 1 0,ω
C3
99
{x } {z} In−2 f − In−2 f (un − vn )ω dx
−1
2 {x } {z} 1 − x C3 In−2 f − In−2 f √ x
√ x (un − vn ) 1 − x2 2
(88)
.
L2ω
Lω
By applying the Poincaré inequality (see [1]) we have that: √ x (un − vn ) C4 un − vn 1 H0,ω 1 − x 2 L2ω
(89)
since |x| 1. The positive constant C4 does not depend on n. Substituting in (88), one obtains: 2 {x } {z} 1 − x , I un − vn H 1 C5 f − I f √ n−2 n−2 0,ω x L2ω
(90)
where the positive constant C5 does not depend on n. { x} {z} For n even the point x = 0 belongs both to the grid x and to the grid z. Therefore, the difference In−2 f − In−2 f { x}
{z}
2 2
) ∈ P02n−1 , the quadrature formula (17) is exact. vanishes at x = 0. Then for the polynomial (In−2 f − In−2 f )2 (1−x x This brings us to: n−1
{x } {x } − x 2 2 2 (1 − zi2 )2 I f − I {z} f 1 √ = I f − f (z ) ωi , (91) i n−2 n−2 n−2 zi x L2ω n i=1,i = 2
{z}
where we noticed that In−2 f = f at the grid z, and that the boundary terms in the quadrature formula are equal to zero. Afterwards, taking into account (65), it is possible to write: n−1
1 2 {x } 2 {x } (1 − zi2 )2 ωi π max In−2 f − f max In−2 f − f (zi ) n x∈[−1,1] zi i = 2 zi n i=1,i = 2 2 1 {x } π max I {x } f − f 2 1 . = π max In−2 f − f (92) x∈[−1,1] x∈[−1,1] n−2 z n2 −1 n The last inequality in (92) holds because, as it was observed in the previous section, in the central part of the interval [−1, 1], the quantity zi−1 behaves like n−1 . Putting all together, we obtain the estimate: 1 { x} (93) u − vn H 1 u − un H 1 + C6 √ f − In−2 f L∞ , 0,ω 0,ω n for some positive constant C6 not depending on n. In conclusion, the spectral convergence, for n → +∞, of vn to the exact solution u comes from the convergence { x} of un to u and of In−2 f to f . This is a well-established result (see, for example, [2]). The error estimate we just proven is not optimal, due to the term n−1/2 appearing on the right-hand side of (93). We believe that such an estimate could be ameliorated, although, at the moment, we have no idea how to do that. An alternative proof, valid independently of the parity of n, can be given by using the first inequality in (62). In √ { x} this case, we can avoid putting x in the denominator of Eq. (88) and estimate directly the L2ω norm of (In−2 f − {z}
In−2 f )(1 − x 2 ) with the help of the quadrature formula. Unfortunately, due to the term the optimality, arriving to the same estimate as in (93). 2
C n
in (62), we lose once again
References [1] C. Canuto, A. Quarteroni, Spectral and pseudo-spectral methods for parabolic problems with nonperiodic boundary conditions, Calcolo 18 (1981) 197–218.
100
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
[2] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods in Fluid Dynamics, Springer, New York, 1988. [3] D. Funaro, Spectral Elements for Transport-Dominated Equations, Lecture Notes in Computational Science and Engineering, vol. 1, Springer, New York, 1997. [4] D. Funaro, A superconsistent Chebyshev collocation method for second-order differential operators, Numerical Algorithms 28 (2001) 151– 157. [5] D. Funaro, Superconsistent discretizations, Journal of Scientific Computing 17 (1–4) (2002) 67–80. [6] D. Funaro, Superconsistent discretizations of integral type equations, Applied Numerical Mathematics 48 (2004) 1–11. [7] D. Funaro, G. Pontrelli, A general class of finite-difference methods for the linear transport equation, Communications in Mathematical Sciences 3 (3) (2005) 403–423. [8] L. Fatone, D. Funaro, V. Scannavini, Finite-difference preconditioners for superconsistent pseudospectral approximations, submitted for publication. [9] F. Marcellán, B.P. Osilenker, I.A. Rocha, Sobolev-type orthogonal polynomials and their zeros, Rendiconti di Matematica VII (17) (1997) 423–444. [10] T.J. Rivlin, The Chebyshev Polynomials, John Wiley & Sons, New York, 1974. [11] G. Szegö, Orthogonal Polynomials, American Mathematical Society, New York, 1939.
A convergence analysis for the superconsistent Chebyshev method Lorella Fatone a,∗ , Daniele Funaro a , Gang Joon Yoon b a Dipartimento di Matematica Pura e Applicata, Università di Modena e Reggio Emilia, Via Campi 213/b, 41100 Modena, Italy b School of Mathematics, Korea Institute for Advanced Study, 207-43 Cheongnyangni 2-dong,
Dongdaemun-gu, Seoul 130-722, Republic of Korea Available online 12 December 2006
Abstract The superconsistent collocation method is based on collocation nodes which are different from those used to represent the solution. The two grids are chosen in such a way that the continuous and the discrete operators coincide on a space as larger as possible (superconsistency). There are many documented situations in which this technique provides excellent numerical results. Unfortunately very little theory has been developed. Here, a theoretical convergence analysis for the superconsistent discretization of the second derivative operator, when the representation grid is the set of Chebyshev Gauss–Lobatto nodes is carried out. To this end, a suitable quadrature formula is introduced and studied. © 2006 IMACS. Published by Elsevier B.V. All rights reserved. PACS: 02.60.Lj; 02.70.Jn Keywords: Convergence; Chebyshev polynomials; Consistency; Collocation grids
1. Introduction An incredible variety of numerical schemes, for the solution of a wide range of problems, can be generated by using two grids at the same time. The first grid is used to describe the approximated solution. As a matter of fact, the computed values are directly related to it. Through a sort of duality process, the second grid is used instead to enforce the equations to be satisfied. There, the residual is required to vanish. This approach allows for the creation of new schemes or the improvement of existing ones. Actually, we may interpret most of the well-known schemes as methods in which the two grids coincide. In order to appreciate the performances of this technique, the two grids must be somehow related. The first one is in general chosen for convenience (because we want the solution to be expressed in that way). The second one should be taken with the aim of obtaining the best from the approximation scheme. A very effective way of choosing the second grid is to build it by asking the discrete operator to be as much “accurate” as possible. For example, when the first grid is fixed, the superconsistent method indicates how to obtain the second grid, by requiring the discretized operator to coincide with the exact operator in a finite dimensional space as larger as possible.
* Corresponding author.
E-mail addresses: [email protected] (L. Fatone), [email protected] (D. Funaro), [email protected] (G.J. Yoon). 0168-9274/$30.00 © 2006 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2006.11.001
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
89
The applications of this idea are unlimited. For the Legendre spectral collocation method in two dimensions we address the reader to [3], where the superiority of the superconsistent method is shown, mainly experimentally, especially for transport-dominated advection–diffusion equations. A survey of possible other applications is provided in [5]. Superconsistent collocation techniques are also considered for integral type equations in [6]. Usually, the extra costs of computing the residual grid and implementing the corresponding discrete operator, are well compensated by the advantages (in term of accuracy and computer time) that one can achieve. So far, except for the results in [7] (dealing with finite-differences for pure hyperbolic problems), very little theory has been developed, despite the large amount of successful experimental material available. This is also true for the simplest examples. This paper, which is quite technical, is aimed to finally provide some theoretical proof of convergence. Some of these results were anticipated in [8]. In that paper Legendre and Chebyshev collocation methods on special (superconsistent) grids were considered for a simple one-dimensional second-order equation. The main interest was to study the properties of suitable preconditioners for the spectral discretization matrices. The technique was successively applied to one-dimensional transport-dominated equations, showing also in that situation the reliability of the method and the effectiveness of the new preconditioners. In order to study the symmetric Chebyshev case, a suitable quadrature formula was introduced. It was noted that such a formula is related to an interesting set of orthogonal polynomials of Sobolev type and conjectures were made about its properties. As a matter of fact, the knowledge of the behavior of nodes and weights of the formula is of crucial relevance for the development of a convergence theory for the approximation of the differential equations examined. The purpose of this paper is to prove these conjectures and add further insight. This will be used to carry out a rigorous convergence analysis for the Chebyshev superconsistent collocation method. Unfortunately, the theory here presented is able to cover only a very simple case, which is marginal with respect to the potentiality of the superconsistent method. In this circumstance, the improvements obtainable by adopting the two-grids approach are irrelevant if compared with the standard Chebyshev collocation method (which can be seen as another two-grids method, having coincident grids). However, together with the results developed for the quadrature formula (that certainly have a general interest for other types of applications), there is now the hope that this preliminary analysis, carried out on a very basic problem, could open the path to more serious generalizations. d2 In order to start, let us briefly recall some general assumptions. Let L be the second derivative operator: L = − dx 2. Then consider the following Dirichlet boundary-value problem: Lu = −u = f
in ]−1, 1[,
with u(−1) = u(1) = 0,
(1)
where f is a given function. The standard Chebyshev collocation method to approximate u is constructed by taking the points of the grid x to be the zeros of Tn . Then, one looks for un ∈ P0n such that: −un (xi ) = f (xi )
for i = 1, . . . , n − 1,
(2)
is the space of polynomials of degree less than or equal to n, vanishing at x = ±1. where In the superconsistent method (see [4]), one introduces another collocation grid z and looks for a new solution vn ∈ P0n in such a way that: P0n
−vn (zi ) = f (zi )
for i = 1, . . . , n − 1.
(3)
Each node zi , i = 1, . . . , n − 1, of the new grid z is required to be solution to the equation: (zi ) = 0, (Lχn+1 )(zi ) = −χn+1
i = 1, . . . , n − 1,
(4)
where
χn+1 = 1 − x 2 Tn ∈ P0n+1 .
(5)
Actually, it is easy to realize that Eq. (4) has exactly n − 1 distinct zeros. This choice is made for the following reasons. If one expands vn in the usual Lagrange basis related to the grid x ∪ {±1}, a (n − 1) × (n − 1) discretization matrix An−1 corresponding to (3) is obtained. Thus, we observe that both the exact operator L and the discrete operator, when applied to χn+1 , are vanishing at the nodes z. The first property is true due to (4), the second one is a consequence of the fact that the matrix An−1 only acts on the values attained at the nodes x , that in the case of the function χn+1 are zero. This means that the kernel of the difference between
90
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
the exact and discrete operators, not only contains the space P0n , but also the enlarged space P0n+1 . This justifies the name superconsistent. Note that we only enlarged the approximation space of one dimension, but this seems to be sufficient for achieving good results in applications. As documented by experiments in [8], in terms of accuracy, the superconsistent method turns out to be superior to the standard collocation method, so that the error |u − vn | is generally less than the error |u − un |. Actually, the Lebesgue constant associated to the distribution of nodes z ∪ {±1} is less than the one corresponding to the nodes x ∪ {±1}. The convergence, for n → +∞, of the discrete solution un , obtained by collocation at the Chebyshev nodes, to the exact solution u of problem (1), is proven in [1], where estimates for the spectral decay of the error are given. Such a proof, based on variational techniques, makes use of the exactness of Gaussian formulas up to degree 2n − 1. In order to analyze the convergence also for the discrete solution vn , it is necessary to introduce and study other suitable high-order integration formulas based on the grid z. This is what will be done in the coming sections. 2. Study of a quadrature formula of Sobolev type are even (odd) functions. Let us start by observing that, if n is odd (even), then χn+1 and χn+1 It should be noted that the grid z is different from the grid x . As a matter of fact, in [8] the following estimates regarding the location of the nodes are given: ⎧ ⎨ xj < zj < ξj +1 , j = 1, . . . , [n/2], j = [n/2] + 1, . . . , n − 1, ξj < zj < xj , (6) ⎩ z = 0, n even and j = n/2, j −1)π , j = 1, . . . , n. However, it is where [∗] is the integer part of ∗ and ξ are the zeros of Tn , i.e.: ξj = − cos (2j 2n interesting to point out that there may be cases in which the grid x and its corresponding superconsistent grid coincide. In fact, in [8], the following theorem was proven. 2
d Theorem 1. Concerning the differential operator L = − dx 2 , with homogeneous Dirichlet boundary conditions, one has that z = x if and only if the points xi , i = 1, . . . , n − 1, are the zeros of the derivative Pn of the Legendre polynomial Pn .
With the help of some trigonometry, another characterization of the grid z can be given. Assuming that, for i = 1, . . . , n − 1, each point of z is written as zi = cos θi , for a certain angle θi , then the following lemma holds. Lemma 2. For each node zi = cos θi , i = 1, . . . , n − 1, the following equation holds: tan nθi =
n cos θi sin θi n2 sin2 θi + 1
,
i = 1, . . . , n − 1.
(7)
Proof. The Chebyshev polynomial Tn (x) = cos nθ , with x = cos θ , satisfies the Sturm–Liouville problem: 1 − x 2 Tn − xTn + n2 Tn = 0. The differentiation of Eq. (8) gives the two relations: 1 − x 2 Tn = − n2 Tn + xTn , 1 − x 2 Tn = − n2 + 1 Tn + xTn . Therefore, from (4) we have: −χn+1 (zi ) = 0 ⇐⇒ n2 + 1 Tn (zi ) + zi Tn (zi ) = 0,
(8)
(9) (10) i = 1, . . . , n − 1.
(11)
Since: n sin nθ , x = cos θ, sin θ n sin nθ cos θ n2 cos nθ Tn (x) = − , 3 sin θ sin2 θ Tn (x) =
(12) x = cos θ,
(13)
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
each θi , i = 1, . . . , n − 1, is solution of: n n2 + 1 sin nθi sin2 θi + n sin nθi cos2 θi − n2 cos nθi cos θi sin θi = 0,
91
(14)
that is:
sin nθi n2 sin2 θi + 1 = n cos nθi cos θi sin θi ,
(15)
2
from which we arrive to (7).
Due to the above considerations, one can see that the difference between two consecutive nodes behaves like n−1 in the central part of the interval [−1, 1], and like n−2 near the endpoints ±1. Then, let us recall some other lemmas, previously obtained in [8], characterizing the superconsistent Chebyshev nodes. Lemma 3. Let χn+1 = (1 − x 2 )Tn . Then, the following orthogonality relation holds: 1
χn+1 qn ω dx = 0,
−1
1 with ω(x) = √ , 1 − x2
(16)
for all polynomials qn of degree n with qn (±1) = 0 and qn (±1) = 0. This lemma provides a link between the Chebyshev superconsistent grid z ∪ {±1} and the so-called Sobolevtype orthogonal polynomials (see, for example, [9]). In fact, given a function f , the set of nodes z is related to the integration formula: 1 f ω dx ≈
n
f (zi )ωi + f (−1)ω˜ 0 + f (1)ω˜ n ,
(17)
i=0
−1
where z0 = −1, zn = 1. The weights ωi , i = 0, . . . , n, and ω˜ 0 , ω˜ n , are defined as follows: 1 ωi =
lˆi ω dx,
i = 0, . . . , n,
(18)
−1
1 ω˜ 0 =
l˜0 ω dx,
1 ω˜ n =
−1
l˜n ω dx,
(19)
−1
where the Lagrange polynomial basis used in (18), (19) is defined as: lˆi ∈ Pn+2 with lˆi (zj ) = δij l˜0 ∈ Pn+2 with l˜0 (zi ) = 0, l˜n ∈ Pn+2 with l˜n (zi ) = 0,
and lˆi (±1) = 0, i = 0, . . . , n, l˜0 (−1) = 1, i = 0, . . . , n,
l˜n (−1) = 0,
(20) l˜0 (1) = 0, l˜n (1) = 1.
(21) (22)
The n + 3 polynomials satisfying conditions (20)–(22) are unique. It is possible to rewrite them with the help of the polynomial of degree n + 1: Ω(x) = (x − z0 ) · · · (x − zn ) = cn x 2 − 1 χn+1 (x),
with cn = −
1 , 2n−1 n2 (n + 1)
(23)
and defining: i (x) =
Ω(x) , − zi )
Ω (zi )(x
i = 0, . . . , n.
(24)
92
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
Thus, it is not hard to show that: ˆl0 (x) = (x − 1) 1 0 (−1)(x + 1) − 1 (x + 3) 0 (x), 2 4 (1 − x 2 ) i (x), i = 1, . . . , n − 1, lˆi (x) = (1 − zi2 ) ˆln (x) = −(x + 1) 1 n (1)(x − 1) + 1 (x − 3) n (x), 2 4 1 l˜0 (x) = 1 − x 2 0 (x), 2 1 l˜n (x) = − 1 − x 2 n (x). 2 Furthermore, since the zeros z0 , z1 , . . . , zn are symmetric with respect to x = 0, that is: zi = −zn−i ,
i = 0, . . . , n,
we come to: lˆi (−x) = lˆn−i (x),
i = 0, . . . , n,
(25)
(26) and l˜0 (−x) = −l˜n (x).
(27)
Consequently, as expected, one has ωi = ωn−i , i = 0, . . . , n, and ω˜ 0 = −ω˜ n . In [8] the following result has been proven. Lemma 4. The quadrature formula (17) is exact if f is a polynomial of degree less than or equal to 2n − 1. Some other properties of Chebyshev polynomials will be useful in the following. In particular for Tn , it is not difficult to get the relation (see [10] or [11]): n (28) 1 − x 2 Tn = {Tn−1 − Tn+1 }. 2 Differentiating the above equation yields the relation: n n + 2xTn − Tn+1 . (29) 1 − x 2 Tn = Tn−1 2 2 On the other hand, by recalling (12), we can write Tn (x) = nUn−1 (x), with x = cos θ , where Un−1 is the Chebyshev polynomial of the second kind, defined by: 1 sin nθ , x = cos θ. (30) Un−1 (x) = Tn (x) = n sin θ Moreover, from (8) one has: n Tn = (xUn−1 − nTn ). (31) 1 − x2 By taking the second derivative in (28), we arrive to: n = {Tn−1 − Tn+1 } χn+1 2
n (n − 1)xUn−2 − (n − 1)2 Tn−1 − (n + 1)xUn + (n + 1)2 Tn+1 . (32) = 2(1 − x 2 ) In addition, the following relations will be also useful (see again [10]): Tn = Un − xUn−1 ,
(33)
U2n+1 = 2Un Tn+1 ,
(34)
xU2n−1 = 1 + 2T2 + · · · + 2T2n−2 + T2n ,
(35)
U2n = 1 + 2T2 + · · · + 2T2n .
(36)
Now the following theorem can be proven.
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
93
Theorem 5. The weights ωi , i = 0, . . . , n in (18) and ω˜ 0 in (19) are positive (as a consequence, ω˜ n is negative). , the corresponding weight ω is positive. Let Q be Proof. First of all, let us show that for any zi = 0, zero of χn+1 i i the polynomial of degree n − 3 defined by
Qi (x) =
(x) χn+1
x 2 − zi2
(37)
.
Then, let fi be the polynomial of degree 2n − 2 defined by: 2 fi (x) = 1 − x 2 Q2i (x).
(38)
It is easy to realize that fi is a non-negative even function such that fi (±1) = 0, fi (zi ) = fi (zn−i ) > 0 and fi (zk ) = 0 for k = i, n − i. The quadrature formula ensures: 1 fi ω dx = ωi fi (zi ) + ωn−i fi (zn−i ) = 2ωi fi (zi ) > 0,
(39)
−1
showing that the corresponding weight ωi must be positive. The next step is to show the positiveness of ω˜ 0 . By definition (see (19)), one has: 1 ω˜ 0 =
1
l˜0 ω dx =
(−1) 4χn+1
−1
1
1 − x 2 (1 − x)ω dx. χn+1
(40)
−1
From (10) we can write: 1 ω˜ 0 = 4χn+1 (−1)
1
2 n + 1 Tn + xTn 1 − x 2 (x − 1)ω dx.
(41)
−1
Since the derivatives of Chebyshev polynomials form a set of orthogonal functions with respect to (1 − x 2 )ω, the following orthogonality relation holds: 1
1 − x 2 ω dx = 0, Tn rn−1
∀rn−1 ∈ Pn−1 ,
(42)
−1
which leads to: ω˜ 0 =
1 (−1) 4χn+1
1
xTn 1 − x 2 (x − 1)ω dx.
(43)
−1
The Sturm–Liouville equation (8) and the orthogonality relations of Chebyshev polynomials (including (42)) bring for n > 2 to: 1 ω˜ 0 = (−1) 4χn+1 1 = 4χn+1 (−1) 1 = 4χn+1 (−1)
1
x xTn − n2 Tn (x − 1)ω dx
−1
1
x 2 Tn (x − 1)ω dx
−1
1 −1
Tn (1 − x)
1 1 − x 2 dx
+ −1
Tn (x
− 1)ω dx
94
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
1 = 4χn+1 (−1)
1
Tn (x − 1)ω dx.
(44)
−1
In order to compute ω˜ 0 explicitly, we first observe that (see [10]): Tn (1) = n2 ,
Tn (1) =
n4 − n2 , 3
and Tn (1) =
n2 (n2 − 1)(n2 − 4) . 15
(45)
Then, using (10), we get:
χn+1 (−1) = (−1)n χn+1 (1) 2(−1)n 2 2 = (−1)n n2 + 1 Tn (1) + Tn (1) = n 2n + 1 . 3 On the other hand, (34)–(36) ensure that: 1
Tn (x − 1)ω dx = (−1)n nπ.
(46)
(47)
−1
Consequently, one recovers the following value of ω˜ 0 : ω˜ 0 =
3π = −ω˜ n . 8(2n3 + n)
(48)
)2 the It is now the time to check that the remaining weights ω0 = ωn are positive. To this end, we apply to f = (χn+1 quadrature formula (17) obtaining:
1
f ω dx = 2ω0 f (1) + ω˜ 0 f (−1) − ω˜ 0 f (1).
(49)
−1
Due to the fact that f (1) > 0 and f (−1) = −f (1) < 0, considering that the integral above is positive and that ω˜ 0 > 0, we deduce the positivity of ω0 . At the same time we entail that ωn is also positive. is allowed to have zm = 0 as a root. Therefore, There is one more case left, corresponding to n = 2m, where χn+1 we have to deal with the weight ωm . When m = 1, it is easy to see that: 1 ω1 =
2 1 − x 2 ω(x) dx > 0.
(50)
−1
Now let us consider the case when m 2. Let us introduce a polynomial f of degree n + 3 = 2m + 3 as: 2 f (x) = 1 − x 2 Q(x),
(51)
where now Q is the polynomial given by: Q(x) =
(x) χ2m+1
. x By applying to f , defined in (51), the quadrature formula (17), we obtain the equation: 1 −1
2 χ2m+1 ω dx = Q(0)ωm . 1 − x2 x
(52)
(53)
Then, considering (23), it is easy to see that: Q(0) = lim Q(x) = x→0
m−1 (−1)m−1 2 zk , cn i=1
(54)
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
where the coefficient cn is provided in (23). Therefore one gets the sign of Q(0): sgn Q(0) = (−1)m .
95
(55)
On the other hand, Eq. (32) and the orthogonality of Ui with respect to (1 − x 2 )ω give: 1 −1
2 χ2m+1 ω dx = m 1 − x2 x
1
−1
1 (2m + 1)2 T2m+1 − (2m − 1)2 T2m−1 1 − x 2 2 dx x 1
= m 4m2 + 1
−1
1 T2m+1 − T2m−1 1 − x 2 2 dx, x
(56)
where we used the relation: T2m+1 + T2m−1 = 2xT2m
(57)
and that the weighted integral of this last quantity is zero. Going ahead, the relations (33), (34) and the orthogonality condition (42) bring to: 1 −1
1 T2m+1 − T2m−1 1 − x 2 2 dx = x
1
−1
1 U2m+1 − U2m−1 1 − x 2 2 dx x
1
1 Um Tm+1 − Um−1 Tm 1 − x 2 2 dx. x
=2 −1
(58)
When m = 2i is even, we have: 1 −1
1 Um Tm+1 − Um−1 Tm 1 − x 2 2 dx = x
1
−1
1 =
1 T2i+1 U2i 1 − x 2 2 dx − x
1
−1
1 2 1 − x 2 2 dx + U2i
−1
U2i−1 1 − x2 T2i ω dx x
1 T2i2 ω dx > 0.
(59)
−1
When m = 2i + 1 is odd, the same arguments bring to: 1 −1
1 Um Tm+1 − Um−1 Tm 1 − x 2 2 dx = x
1
−1
U2i+1 T2i+2 ω dx − 1 − x2 x
1
−1
U2i
1 T2i+1 1 − x 2 2 dx < 0. x
Hence we have shown that 1 χ 2 2 2m+1 ω dx = (−1)m , 1−x sgn x
(60)
(61)
−1
which, from (53), implies ωm > 0. This completes the proof.
2
In addition, the following results can be proven. Lemma 6. For each n 2, one has the two inequalities: n 1 1
2 2 C 2 2 φ ω dx φ (zi )ωi + φ (−1)ω˜ 0 + φ (1)ω˜ n 2 φ 2 ω dx, n −1
i=0
−1
for all φ polynomials of degree n, where the positive constant C does not depend on n.
(62)
96
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
Proof. First, let us show the second inequality. Since φ can be written as a multiple of Tn plus lower degree terms, due to the exactness of the quadrature formula up to degree 2n − 1, it is enough to prove the inequality for φ = Tn . Therefore, it is necessary to show that: n
Tn2 (zi )ωi
+ 2Tn (−1)Tn (−1)ω˜ 0
+ 2Tn (1)Tn (1)ω˜ n
i=0
1 B
Tn2 ω dx,
(63)
−1
where B does not depend on n. Since Tn (x)Tn (x) is always an odd function and ω˜ n = −ω˜ 0 , the part containing the boundary terms on the left-hand side of (63) is less than zero: 2Tn (−1)Tn (−1)ω˜ 0 + 2Tn (1)Tn (1)ω˜ n = 4Tn (−1)Tn (−1)ω˜ 0 < 0.
(64)
Thus, we just have to deal with the summation. Let us note that, by taking f = 1 in (17) , we have: 1 ω dx = π =
n
(65)
ωi .
i=0
−1
Moreover, since |Tn (x)| 1, ∀x ∈ [−1, 1] and ωi > 0, i = 0, . . . , n, the sum in (63) can be easily bounded. As a matter of fact: n
Tn2 (zi )ωi
i=0
n
i=0
π ωi = π = 2 · = 2 2
1 Tn2 ω dx.
(66)
−1
Therefore (63) is true by choosing B = 2. In order to show the first inequality in (62), due to the exactness of the quadrature formula (17) up to degree 2n − 1, )2 (1 − x 2 ) ∈ P . Therefore, it is necessary to have: we now choose φ 2 = (χn+1 2n C n
1
2 (χn+1 )2 1 − x 2 ω dx 4 χn+1 (1) ω˜ 0 .
(67)
−1
Such an inequality is correct since Eqs. (46)–(48) ensure that the right-hand side of (67) behaves like n5 . On the other hand, the left-hand side of (67) behaves like n6 . In fact Eq. (10) entails: 1
(χn+1 )2
1 − x 2 ω dx =
−1
1
2 2 n + 1 (Tn )2 + (xTn )2 + 2x n2 + 1 Tn Tn 1 − x 2 ω dx,
(68)
−1
and by taking only the first integral in (68) one has: 1
2 2 2 n + 1 (Tn )2 1 − x 2 ω dx = n2 n2 + 1
−1
that grows as n6 .
π
2 π , sin2 (nθ ) dθ = n2 n2 + 1 2
(69)
0
2
We continue by collecting other results about our quadrature formula. In particular the asymptotic behavior of the weights ω0 (= ωn ) and ω˜ 0 (= −ω˜ n ) is investigated. Let us start by defining a function λ on positive integers: 1, if n is even, λ(n) = (70) 0, otherwise. Lemma 7. For the Chebyshev polynomials Tn , the following relations can be proven:
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
1
Tn ω dx =
−1
1
n3 πλ(n), 2
xTn ω dx =
−1
1
(71)
n(n2 − 1) πλ(n − 1), 2
(72)
n(n2 − 2) πλ(n). 2
x 2 Tn ω dx =
−1
97
(73)
Proof. First of all, we show (72). For n = 2k + 1, the identities (35) and (36) imply that 1
xTn ω dx
1 =n
−1
xU2k ω dx
1
1 k k
= 2n x T2i ω dx = 2n 2i xU2i−1 ω dx = 2nk(k + 1)π.
−1
−1
i=1
i=1
(74)
−1
= (xU2k−1 ) − U2k−1 , the same arguments Since k = 12 (n − 1), we have obtained (72). Using the identity xU2k−1 show (73). Finally, the identity (71) follows directly from (73) and:
1
1 − x 2 Tn ω dx = n
−1
1 (xUn−1 − nTn )ω dx = nπλ(n),
(75)
−1
2
where Eq. (31) is used.
The asymptotic behavior of the two weights ω0 (= ωn ) and ω˜ 0 (= −ω˜ n ) is established by the following lemma. Lemma 8. For n large, up to multiplicative constants, one has: ω0 ∼ n−1
and ω˜ 0 ∼ n−3 .
(76)
Proof. Eq. (48) ensures that ω˜ 0 ∼ n−3 . Now the estimation of the behavior of ω0 is required. Using the identity (32) and the expression of the polynomial n in (24): n =
(x + 1)χn+1 (1) 2χn+1
(77)
,
it is possible to rewrite lˆn as: lˆn = − n (1) + 1/2 l˜n +
n 2 (1) (x + 1) (Tn−1 − Tn+1 ). 8χn+1
(78)
Hence, the weight ω0 is given by: ω0 = −
n (1) +
1 2 1 n ω˜ n + − Tn+1 )ω dx. x + 2x + 1 (Tn−1 2 8χn+1 (1)
(79)
−1
On the other hand, using the expression for χn+1 in (10), the properties of Chebyshev polynomials given in (45), and Lemma 7, for n large, up to multiplicative constants, one has:
n (1)ω˜ n
∼n
−1
,
n (1) 8χn+1
1
2 x + 2x + 1 (Tn−1 − Tn+1 )ω dx ∼ n−1 .
(80)
−1
It turns out that the asymptotic behavior of ω0 is n−1 . This completes the proof.
2
98
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
3. Convergence analysis The convergence of the discrete solution vn of (3) (obtained by collocation at the superconsistent Chebyshev grid) to the solution u of problem (1) can now be proven. This result was not available in [8], although there it was shown experimentally that the error |u − vn | decays faster than the error |u − un |, and this last quantity is known to converge to zero exponentially. The positivity of the weights of the quadrature formula (17), resulting from Lemma 5, plays an important role in the theorem that follows. Theorem 9. The discrete superconsistent Chebyshev solution vn of problem (3) converges, for n → +∞, to the exact solution u of problem (1) and the rate of convergence is of spectral type. {z}
Proof. Let In−2 f be the interpolant of f at the nodes z, i.e.:
{z} In−2 f (zi ) = f (zi ),
i = 1, . . . , n − 1.
(81)
From (3) we can write: {z}
−vn = In−2 f
in ]−1, 1[,
(82)
{z}
since vn and In−2 have both degree n − 2 and coincide in n − 1 points. { x}
Similarly, by denoting by In−2 f the interpolant of f at the nodes x , i.e.:
{ x} In−2 f (xi ) = f (xi ),
i = 1, . . . , n − 1,
(83)
we deduce that: { x}
−un = In−2 f
in ]−1, 1[.
(84)
1 = H 1 (−1, 1) in the usual way, by the After introducing the weighted Sobolev spaces L2ω = L2ω (−1, 1) and H0,ω 0,ω triangle inequality, one gets:
u − vn H 1 u − un H 1 + un − vn H 1 . 0,ω
0,ω
0,ω
(85)
It is already known that the first term on the right-hand side of (85) converges spectrally to zero, so that it is now necessary to estimate the second one. From (82) and (84) we get: { x}
{z}
−(un − vn ) = In−2 f − In−2 f
in ]−1, 1[.
(86)
By using an inequality proven in [1], Eq. (86) gives: 1 un − vn 2H 1 C3 0,ω
(un − vn ) (un − vn )ω dx
−1
1 = −C3
(un − vn ) (un − vn )ω dx
−1
1 = C3
{x } {z} In−2 f − In−2 f (un − vn )ω dx,
(87)
−1
where C3 is a constant not depending on n. The integration by parts takes into account that (un − vn )(±1) = 0. At this point, it is possible to continue with two alternative proofs. The first one only works for n even. Thanks to the Schwarz inequality, we can write:
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
1 un − vn 2H 1 0,ω
C3
99
{x } {z} In−2 f − In−2 f (un − vn )ω dx
−1
2 {x } {z} 1 − x C3 In−2 f − In−2 f √ x
√ x (un − vn ) 1 − x2 2
(88)
.
L2ω
Lω
By applying the Poincaré inequality (see [1]) we have that: √ x (un − vn ) C4 un − vn 1 H0,ω 1 − x 2 L2ω
(89)
since |x| 1. The positive constant C4 does not depend on n. Substituting in (88), one obtains: 2 {x } {z} 1 − x , I un − vn H 1 C5 f − I f √ n−2 n−2 0,ω x L2ω
(90)
where the positive constant C5 does not depend on n. { x} {z} For n even the point x = 0 belongs both to the grid x and to the grid z. Therefore, the difference In−2 f − In−2 f { x}
{z}
2 2
) ∈ P02n−1 , the quadrature formula (17) is exact. vanishes at x = 0. Then for the polynomial (In−2 f − In−2 f )2 (1−x x This brings us to: n−1
{x } {x } − x 2 2 2 (1 − zi2 )2 I f − I {z} f 1 √ = I f − f (z ) ωi , (91) i n−2 n−2 n−2 zi x L2ω n i=1,i = 2
{z}
where we noticed that In−2 f = f at the grid z, and that the boundary terms in the quadrature formula are equal to zero. Afterwards, taking into account (65), it is possible to write: n−1
1 2 {x } 2 {x } (1 − zi2 )2 ωi π max In−2 f − f max In−2 f − f (zi ) n x∈[−1,1] zi i = 2 zi n i=1,i = 2 2 1 {x } π max I {x } f − f 2 1 . = π max In−2 f − f (92) x∈[−1,1] x∈[−1,1] n−2 z n2 −1 n The last inequality in (92) holds because, as it was observed in the previous section, in the central part of the interval [−1, 1], the quantity zi−1 behaves like n−1 . Putting all together, we obtain the estimate: 1 { x} (93) u − vn H 1 u − un H 1 + C6 √ f − In−2 f L∞ , 0,ω 0,ω n for some positive constant C6 not depending on n. In conclusion, the spectral convergence, for n → +∞, of vn to the exact solution u comes from the convergence { x} of un to u and of In−2 f to f . This is a well-established result (see, for example, [2]). The error estimate we just proven is not optimal, due to the term n−1/2 appearing on the right-hand side of (93). We believe that such an estimate could be ameliorated, although, at the moment, we have no idea how to do that. An alternative proof, valid independently of the parity of n, can be given by using the first inequality in (62). In √ { x} this case, we can avoid putting x in the denominator of Eq. (88) and estimate directly the L2ω norm of (In−2 f − {z}
In−2 f )(1 − x 2 ) with the help of the quadrature formula. Unfortunately, due to the term the optimality, arriving to the same estimate as in (93). 2
C n
in (62), we lose once again
References [1] C. Canuto, A. Quarteroni, Spectral and pseudo-spectral methods for parabolic problems with nonperiodic boundary conditions, Calcolo 18 (1981) 197–218.
100
L. Fatone et al. / Applied Numerical Mathematics 58 (2008) 88–100
[2] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods in Fluid Dynamics, Springer, New York, 1988. [3] D. Funaro, Spectral Elements for Transport-Dominated Equations, Lecture Notes in Computational Science and Engineering, vol. 1, Springer, New York, 1997. [4] D. Funaro, A superconsistent Chebyshev collocation method for second-order differential operators, Numerical Algorithms 28 (2001) 151– 157. [5] D. Funaro, Superconsistent discretizations, Journal of Scientific Computing 17 (1–4) (2002) 67–80. [6] D. Funaro, Superconsistent discretizations of integral type equations, Applied Numerical Mathematics 48 (2004) 1–11. [7] D. Funaro, G. Pontrelli, A general class of finite-difference methods for the linear transport equation, Communications in Mathematical Sciences 3 (3) (2005) 403–423. [8] L. Fatone, D. Funaro, V. Scannavini, Finite-difference preconditioners for superconsistent pseudospectral approximations, submitted for publication. [9] F. Marcellán, B.P. Osilenker, I.A. Rocha, Sobolev-type orthogonal polynomials and their zeros, Rendiconti di Matematica VII (17) (1997) 423–444. [10] T.J. Rivlin, The Chebyshev Polynomials, John Wiley & Sons, New York, 1974. [11] G. Szegö, Orthogonal Polynomials, American Mathematical Society, New York, 1939.