Barycentric prolate interpolation and pseudospectral differentiation

In this paper, we provide further illustrations of prolate interpolation and pseudospectral differentiation based on the barycentric perspectives. The convergence rates of the barycentric prolate interpolation and pseudospectral differentiation are derived. Furthermore, we propose the new preconditioner, which leads to the well-conditioned prolate collocation scheme. Numerical examples are included to show the high accuracy of the new method. We apply this approach to solve the second-order boundary value problem and Helmholtz problem.


Introduction
In many cases, we are confronted with wave phenomena, such as wave scattering, signal processing, and antenna theory, which are characterized by bandlimited functions (whose Fourier transforms are compactly supported) [3,14]. It is well known that the natural tool for effectively representing bandlimited functions on an interval is prolate spheroidal wave functions (PSWFs) [7,28,32]. Hence, there has been a growing interest in developing prolate spheroidal wave functions, which also offer an alternative to Chebyshev and other orthogonal polynomials for pseudospectral/collocation and spectral-element algorithms [5].
It is well known that a simple way of approximating a function f (x) is to choose a sequence of points {x j } n j =0 and find the function P (x) from the values of f (x) at these interpolation nodes, i.e., set P (x j ) = f (x j ), 0 ≤ j ≤ n. The standard Yan Tian yantian@csrc.ac.cn tool for interpolation and approximation algorithms was investigated in [15,32]. We highlight that a very popular alternative nowadays is to use barycentric interpolation formula, and the favourable numerical aspects of this way are summarized by Berrut and Trefethen [1,24]. So, the related issues are worthy of investigation.
The purpose of this paper is to have new insights into prolate interpolation and pseudospectral differentiation based on the Prolate-Gauss-Lobatto points. The inspiration behind the proposed numerical method is the remarkable advantage of barycentric interpolation formula [1]. The main contributions reside in the following aspects.
• We give the barycentric prolate interpolation and differentiation formula, which enjoys a more stable approximability and efficiency than the formulas given early [32,34]. • We give the error analysis of the barycentric prolate interpolation and differentiation based on the error analysis of the standard prolate interpolation and differentiation [22]. • We offer a preconditioning matrix that nearly inverts the second-order prolate pseudospectral differentiation matrix, leading to a well-conditioned collocation approach for second-order boundary value problems.
The structure of the paper is as follows. In Section 2, we review some results of PSWFs and prolate interpolation and pseudospectral differentiation, while Section 3 combines the barycentric form with the prolate interpolation, which yields the new barycentric prolate interpolation and pseudospectral differentiation scheme. Furthermore, the convergence analysis of barycentric prolate interpolation and pseudospectral differentiation is given. In Section 4, we introduce the preconditioning matrix that nearly inverts the second-order barycentric prolate differentiation matrix. Sections 5 and 6 demonstrate the analysis via several numerical experiments and apply to the second-order boundary value problem and Helmholtz problem.

Preliminaries
The PSWFs were introduced in the 1960s by Slepian et al. in a series of papers [20,21]. Firstly, we briefly recall some preliminary properties of the PSWFs. All of these can be found in [2, 6, 13-16, 18, 20, 21, 30, 31, 34]. Prolate spheroidal wave functions of order zero are the eigenfunctions ψ j (x) of the Helmholtz equation in prolate spheroidal coordinates: for x ∈ (−1, 1) and c ≥ 0. A series of papers [18,32,34] have shown that the ψ j (x) are also the eigenfunctions of the following integral eigenvalue problem: Here, {χ j := χ j (c)} ∞ j =0 and {λ j := λ j (c)} ∞ j =0 are the associated eigenvalues corresponding to the differential operator and integral operator, and the constant c is known as the "bandwidth parameter." The eigenvalues {λ j } ∞ j =0 satisfy |λ 0 | > |λ 1 | > |λ 2 | > . . . > 0, which decay exponentially to nearly 0. Specifically, based on the work of Wang and Zhang [31], we have: where the notation A B means that for B = 0, the ratio A B → 1 in the sense of some limiting process.
An important issue related to the PSWFs is the choice of bandlimit parameter c. For general functions, we do not have a simple optimal c. This is due to the fact that an arbitrary function has many different modes and each mode has a distinct optimal c [7]. Regardless of whether the function being represented is bandlimited or not, all the useful choices of c must satisfy [5]: As recommended in [30,32], in practice, a quite safe choice is c = N 2 . With this in mind, we choose c = N 2 . Guidelines on the suitable choice of c can be found in [3]. The practical rule for pairing up (c, N) has been given in [13,32].
Denoting the zeros of (1 − x 2 )ψ n (x) by the Prolate-Gauss-Lobatto points (PGL points). For computation, Boyd [6] described Newton's iteration method with some care in selecting initial guesses. Generally, [12] gives the efficient algorithm for computing zeros of special functions, such as PSWFs. With the Prolate-Gauss-Lobatto points at our disposal, we will introduce the prolate interpolation and pseudospectal differentiation.

Prolate interpolation and pseudospectal differentiation
In this subsection, firstly, we review some facts about prolate interpolation and pseudospectal differentiation.
The key idea for interpolation is to search the prolate cardinal functions i (x) := i (x; c), which are designed to satisfy the interpolation property: Then, the function f (x) is approximated by The standard route to get the derivatives is by directly differentiating the prolate cardinal basis i (x).
Generally, we define the prolate cardinal functions i (x) as where {x k } N k=0 are the Prolate points, which are zeros of s p (x). It follows that the standard interpolation is [15,32]: The standard differentiation generated from the cardinal basis (2.7) can be computed by: In the following, let us consider the prolate interpolation and pseudospectal differentiation scheme through the barycentric form.

Barycentric prolate interpolation and differentiation
In this section, we start with the barycentric interpolation [1,9,11], which are important pieces of the puzzle for our new approach, and then give the new insights into prolate interpolation, which are called the barycentric prolate interpolation. The differentiation matrices are derived through the barycentric interpolation formula. Then, the convergence analysis of barycentric prolate interpolation and differentiation is given.

Barycentric interpolation formula
Let {x j } N j =0 be a set of distinctive nodes in [−1, 1], which are arranged in ascending order, together with corresponding numbers f (x j ). We assume that the nodes are real and they are zeros of the function s(x), i.e., s(x j ) = 0 for 0 ≤ j ≤ N. Thus, the Lagrange interpolating basis is defined by: Accordingly, the interpolation in Lagrange form for the function f (x) is The barycentric formula is the alternative Lagrange form, and for computations, it is generally recommended that one should use barycentric interpolation formula [1,11], which has stability or robustness property that proves advantageous in some application. The barycentric interpolation is defined as are the barycentric weights. To this end, it suffices to note that the barycentric weights {w j } N j =0 can be written as different quantity. As with the polynomial interpolation [1,24], s(x) can be written as: such that barycentric weights become .
For certain special sets of nodes {x j } N j =0 , the explicit expressions of the barycentric weights w j were available in [1,26,27]. For general point sets, the barycentric weights w j can be evaluated by the fast multipole method [23]. These observations lead to an efficient method for computing prolate interpolants based on the Prolate-Gauss-Lobatto points through a new definition of non-zero barycentric weights.

Barycentric prolate interpolation and pseudospectral differentiation
Using the barycentric form, this subsection give a new definition of barycentric prolate weights, which leads to remarkably simple and efficient schemes for the construction of rational barycentric interpolation, which is denoted by barycentric prolate interpolation.
The fist question is how to choose barycentric weights. In the similar manner as deriving the barycentric formula (3.13) form Lagrange interpolation (3.11), since the Prolate-Gauss-Lobatto points are the roots of (1 − . Correspondingly, we define the prolate barycentric weights to be . (3.14) According to the foregoing observations, it is desirable to define a new interpolation which is called barycentric prolate interpolation. Moreover, the interpolation property is stable with respect to the nonzero weights, as noticed in [33].

Definition 3.1 The barycentric prolate interpolation can be expressed as
are the Prolate-Gauss-Lobatto points and w j = .
The error analysis will be derived in the next subsection. In fact, from the numerical evidences in Figs

Remark 2
The barycentric formula has natural advantages for applications to fast multipole method [1], which is an useful and efficient tool to improve the complexity of centain sums (3.15) The idea of using FMM to accelerate the interpolation and pseudospectral differentiation can be traced back to [4,8] and we see from [15] that the FMM was used to accelerate the standard prolate interpolation and differentiation. It is noteworthy to point out that the new scheme (3.15) can also be accelerated by the FMM through a very similar process in [17].
is normalized Legendre polynomial and α k is the eigenvector of a matrix, which is complex and time-consuming. However, it is obvious that the factor s p (x) has dropped out in the (3.15), and this feature has practical consequences.
Furthermore, defining the cardinal basis function of the barycentric prolate interpolation (3.15) as It leads to the differentiation matrices which have the explicit formulas [23]: are the Prolate-Gauss-Lobatto points and w j = .

Remark 4
It is obvious that the standard differentiation method (2.10) involves the first-order differentiation value j (x i ) (2.9), which causes error propaganda for large number N. The barycentric prolate differentiation only involves the barycentric weights value. Hence, the barycentric prolate differentiation form is stable even for large N, which has been shown in Fig. 6.

Convergence properties of barycentric prolate interpolation and differentiation
Results can also be obtained for the convergence properties of barycentric prolate interpolation and differentiation.
Lemma 3.1 [22] Let f be the entire function, R be the boundary of the square is the interpolant of f (x) at the Prolate-Gauss-Lobatto points (2.8), then it follows for χ n > c 2 and −1 < x < 1 that Remark 5 We remark that the condition "Let f be the entire function" in Lemma 3.1 can be refined as "Let f be analytic in a region bounded by the square [18,35].

Theorem 3.1 Let f be analytic in a region bounded by the square
is the interpolant of f (x) at the Prolate-Gauss-Lobatto points by fomula (2.8) and (3.15), then it follows for χ n > c 2 and −1 < x < 1 that Proof Due to Lemma 3.1, when P n [1](x) interpolates the constant function f (x) = 1, let P n [1](x) = 1 + E n (x), we provide the error that (3.25) Then, we have: Combining with (3.25), we obtain: where ε n and ε n are defined in (3.24). The proof is completed. Remark 7 A function f may be less smooth than the case we have considered; numerical results illustrate that it might be also suited to this fast convergence. However, it appears open to know about exactly how the convergence rates of barycentric prolate interpolation depend on the degree of smoothness of f .

A well-conditioned prolate-collocation method
As everyone knows, the second-order prolate differentiation matrix is apparently unstable even for slightly large N [19]. Fortunately, Wang et al. [32] offered a new basis leading to well-conditioned collocation linear systems. In this subsection, we give a different way to evaluate the Birkhoff interpolation basis, which generates the preconditioner P in , such that the eigenvalues of P in D (2) in are nearly concentrated around one.
Consider the second-order BVPs with Dirichlet boundary conditions: Following the work of Wang [32], the Birkhoff interpolation p(x) of f (x) can be uniquely determined by: are the Birkhoff interpolation basis and satisfy: and We omit the proof, since it is very similar to that in [29]. In order to avoid the instability and low-efficiency of the Lagrange interpolation, the barycentric form is used which is recommended by [1].
To construct the Birkhoff interpolation basis, we give the numerical scheme for integral (4.32) Introducing the change of variable allows us to rewrite the definite integrals (4.34) further as Since the integrands in (4.36) can be computed exactly using an Gauss quadrature at Legendre points. Based on fast O(N) operations for the computation of Gaussian quadrature due to Hale and Townsend [10], we get the fast scheme for the Birkhoff j (x i ), and define the matrices (4.37) Due to (4.27), h k (x) in (3.16) can be approximated by According to the fact that h k (x) satisfying h j (x i ) = δ ij , it follows that where I M is an M × M identity matrix, and the matrix D (2) in is the same as in (3.17). We depict in Fig. 1 the distribution of the largest and smallest eigenvalues of B in D (2) in at the Prolate-Gauss-Lobatto points. This agrees with (4.39).

Remark 8
Obviously, the new system (4.42) does not involve the direct multiplication of the preconditioner, and the round-off errors in forming differentiation matrices can be alleviated.

Remark 9
The use of Birkhoff interpolation as basis functions for deriving precondition is mimic to the preconditioning technique in [32]. However, [32] search for the Birkhoff interpolation basis {B j (x)} through expansion in a different finite dimensional space, and then solving the coefficients by the interpolation conditions. This process involves inverting a matrix of PSWF values, which is time-consuming. My idea of constructing the basis {B j (x)} in (4.31)-(4.32) is actually inspired by polynomial-based algorithms in [29] and the new insights reside in two aspects. First, in order to avoid the instability of the Lagrange interpolation, the barycentric form was used. Second, through changing the variable, the integrals in (4.32) were computed by the fast Gauss quadrature proposed by Hale and Townsend.

Numerical tests
In this section, we illustrate the numerical results in this paper. All the numerical results in this paper are carried out by using Matlab R2014a on a desktop (4.0 GB RAM, 2 Core2 (64 bit) processors at 3.17 GHz) with Windows 7 operating system. Figure 2 illustrates the convergence of the barycentric prolate interpolation formulas for the two analytic functions:

Example 1
f (x) = e sin(6x) and f (x) = 2 sin(10x). For each n, the error is defined by which is measured at 1000 random points in [−1, 1]. As we can see, the convergence is exponential and is almost indistinguishable for different c. Moreover, it is shown that the optimal c depends on the function being approximated [7]. In the following, we will take c = n/2 for general functions, which is recommended in [30,32].  Figs. 3 and 4. It is seen that the errors for these approaches decrease very fast. Furthermore, the barycentric prolate interpolation has better stability than that of the Lagrange formulation for a large number of points.

Example 3
For the wave functions f (x) = sin(25x) 2−x 2 and f (x) = (cos(25x) + sin(x))/(1 + 4x 2 ), we compare the barycentric prolate interpolation (c = n/2) with the barycentric interpolation in the polynomial case, whose nodes and barycentric weights are computed in the chebfun system by the command legpts [24]. Figure 5 illustrates the barycentric prolate interpolation yields spectrally accurate results using even fewer points than barycentric interpolation in the polynomial case. at Prolate-Gauss-Lobatto points by the barycentric prolate differentiation (3.18)- (3.19) and standard method (2.9)-(2.10). Results of these calculations are shown in Fig. 6. As can be seen, since the standard method (2.10) involves the first-order differentiation value, it causes error propaganda for a large number n. There is a good performance of prolate barycentric differentiation, which gives us the motivations for the application.

Application
Different from the usual collocation scheme using the standard Lagrange differentiation, barycentric prolate differentiation (3.18)- (3.19), combining with the usual spectral collocation method and GMRES, has been implemented and tested on the highly oscillatory problem and two-dimensional Helmholtz problem. The comparison with CC points-based method is reminiscent when the solution is highly oscillatory.

Example 5
The second example is one where the solution is very oscillatory The behavior of the prolate barycentric differentiation matrix is demonstrated in Fig. 7. It is clear that this method is rapidly convergent and stable, which is better than the usual collocation method based on CC points.  Example 6 We extend the barycentric prolate pseudospectral method to 2D Helmholtz problem [25], which arises in the analysis of wave propagation: where u = 0 on the boundary and k is a real parameter. For such a problem, we set up a grid based on Prolate-Gauss-Lobatto points independently in each direction called a tensor product grid. To solve such a problem for the particular choices k = 9, f (x, y) = exp(−10[(y − 1) 2 + (x − 1/2) 2 ]). The solution appears as a mesh plot in Fig. 8. Compared with the value u(0, 0) is accurate to nine digits at Chebyshev grid [25] when N = 24, the new barycentric prolate differentiation scheme (3.19) achieves the accuracy to eleven digits at the same number of points. On the right side of Fig. 8, the absolute error at u(0, 0) is illustrated when N = 4 : 2 : 38, which show the fast convergence rate at Prolate-Gauss-Lobatto points.  Table 1 compares the condition number and errors of the spectral collocation (SC) scheme (4.40), direct preconditioned (M-PC) scheme (4.41), and the new basis preconditioned collocation (B-PC) scheme (4.42), respectively. We also show the iteration number for solving the systems by GMRES. Table 1 clearly indicates that the two preconditioned schemes are well-conditioned and the new basis preconditioned collocation (B-PC) scheme has desired performance.

Conclusion
In this paper, we have developed a new scheme for the prolate interpolation and prolate spectral differentiation. The solver is based on the barycentric interpolation, which allows for stable approximation and the error analysis of barycentric prolate interpolation and differentiation are given. What's more, the new preconditioning skill is proposed for the usual prolate-collocation scheme. The numerical examples demonstrate the performance of the proposed algorithms.