Abstract
A new class of spectral methods for solving two-point boundary value problems for linear ordinary differential equations is presented in the paper. Although these methods are based on trigonometric functions, they can be used for solving periodic as well as nonperiodic problems. Instead of using basis functions periodic on a given interval , we use functions periodic on a wider interval. The numerical solution of the given problem is sought in terms of the half-range Chebyshev-Fourier (HCF) series, a reorganization of the classical Fourier series using half-range Chebyshev polynomials of the first and second kind which were first introduced by Huybrechs (2010) and further analyzed by Orel and Perne (2012). The numerical solution is constructed as a HCF series via differentiation and multiplication matrices. Moreover, the construction of the method, error analysis, convergence results, and some numerical examples are presented in the paper. The decay of the maximal absolute error according to the truncation number for the new class of Chebyshev-Fourier-collocation (CFC) methods is compared to the decay of the error for the standard class of Chebyshev-collocation (CC) methods.
1. Introduction and Formulation
A standard problem in approximation theory is to compute the coefficients of a Fourier series to approximate smooth and periodic functions. This can be efficiently done by using the FFT, which is a stable and well-understood method that yields spectral convergence. Things look very different when dealing with nonperiodic or nonsmooth functions. This is due to the so called Gibbs phenomenon which causes oscillations near the points of discontinuity and/or near the boundary as well as slow decay of Fourier coefficients. There are several possibilities to overcome these difficulties and successfully deal with such functions. Some of them were analyzed by Gottlieb and Shu in [1] and by Tadmor in [2]. One possibility is to use some periodizing transformation and compute the Fourier series of the transformed function. There is a widely used transformation which yields Chebyshev polynomials of the first kind. Much about this can be found in the literature by Boyd in [3], by Fornberg and Sloan in [4], or by Trefethen in [5]. These polynomials are arranged as a Chebyshev series.
Another approach was recently presented by Huybrechs in [6] where he analyzed the problem which was stated by Boyd in [7] and by Bruno et al. in [8].
Problem 1. For , let be the space of -periodic functions of the form The Fourier extension of , defined on the interval to the interval , is the solution to the optimization problem
In his paper, Huybrechs considered a square-integrable function that is not necessarily smooth or periodic. The idea to obtain a spectrally accurate Fourier series is to extend the given function to a function that is periodic on a wider interval , . The Fourier series of the constructed function is obviously point-wise convergent to on the interval .
Huybrechs analyzed this problem for the choice and developed three numerical methods for solving it. Besides proving the existence and uniqueness of the solution, he characterized the solution with two nonclassical families of orthogonal polynomials related to Chebyshev polynomials of the first and second kind: the half-range Chebyshev polynomials of the first and second kind. The second and third methods are based on projection and collocation. The convergence is in most cases exponential rather than superalgebraic.
The half-range Chebyshev polynomials of the first and second kind and the corresponding half-range Chebyshev-Fourier (HCF) series were further explored by Orel and Perne in [9]. The authors presented an efficient method for the construction of these polynomials using the modified Chebyshev algorithm for the computation of the recursion coefficients via the three-term recurrence relation formula. Moreover, the authors developed some necessary tools for the construction of spectral methods using these polynomials.
In this paper, instead of approximation problems, like Problem 1, which was successfully solved by Huybrechs in [6], we are interested in solving BVPs in ODEs via spectral methods, that is, solving problems as defined below.
Problem 2. Find the numerical solution of a linear two-point boundary value problem of the form with boundary conditions where is a linear differential operator: is the identity operator, and is a set of linear boundary differential operators.
There are several different classes of numerical methods to solve boundary value problems. In the company of such methods as finite differences (FDM) and finite elements (FEM) we are interested in spectral methods (SM). It is well known that spectral methods approximate the solution in a finite dimensional subspace of a Hilbert space. The basis functions used are defined globally (on the whole interval). On the contrary, the basis functions for finite element methods are defined locally (only on a small interval).
Different approaches are used if the underlying problem is periodic or nonperiodic. In the case of periodic problems, the natural basis functions are trigonometric functions. In other words, the solution is approximated with a truncated Fourier series. Equidistant mesh points are used to discretize the interval on which the problem is defined. In the case of nonperiodic problems, orthogonal polynomials are used, especially the Chebyshev polynomials of the first kind. In this case, the solution is approximated with a truncated Chebyshev series. There are two types of Chebyshev mesh points. Chebyshev points of the second kind, where the Chebyshev polynomials of the first kind reach their extreme values, are generally used in spectral methods to discretize the underlying interval. Chebyshev points of the first kind are the zeros of the Chebyshev polynomials of the first kind. Usually, Chebyshev points in nonperiodic problems outweigh equidistant points because they are denser near the boundary of the interval. Besides, such a distribution of points overwhelms problems caused by Gibbs and/or Runge phenomenon. There are different types of spectral methods depending on the method used for computing the expansion coefficients, for example, Galerkin method, Tau method, or collocation. Spectral methods based on collocation are usually called pseudospectral methods. Classical references on spectral methods include textbooks by Boyd, [3], Canuto et al. [10], Fornberg [11], Gottlieb and Orszag [12], and Trefethen [5], and more recent ones are Canuto et al. [13] and Shen et al. [14].
There were several attempts to solve nonperiodic problems using a trigonometric basis. Adcock in [15, 16] solved the problem with modified Fourier series, using Galerkin method to compute expansion coefficients. Huybrechs in [6] proposed a new set of trigonometric functions, which includes sines and cosines as well as half-sines and half-cosines.
Our intention is to provide a new class of spectral methods using approaches described in [6, 9], combined with collocation to compute expansion coefficients. This approach yields pseudospectral method for solving nonperiodic problems with tools used for solving periodic problems. We focus on Dirichlet boundary conditions, although it is not difficult to extend this approach to Neumann or mixed boundary conditions. The generalization to higher order linear boundary problems is also possible. The restriction to the interval is only a matter of simplification, since it is well known how to map an arbitrary interval to the unit interval.
The paper is organized as follows. In Section 2 we define the form of the exact solution to the discrete problem by introducing the set of basis functions proposed by Huybrechs in [6] and two nonclassical families of orthogonal polynomials. In Section 3 follows the developement of a spectral method for the solution of the boundary value problem, based on a pseudospectral (collocation) approach. Error analysis and convergence results based on error analysis approach in [13, 14] are addressed in Section 4. Finally, we present some numerical examples in Section 5. Section 6 concludes the paper.
2. Approximation with Half-Range Chebyshev Polynomials
Let us first review the exact solution of the discrete Problem 1 using the half-range Chebyshev polynomials. The idea is to extend a nonperiodic function on the interval to an interval on which the given function is periodic and to use the set of trigonometric functions, proposed by Huybrechs in [6]. He analyzed the extension for and proposed to use the set of basis functions where Note that the set consists of even and the set of odd functions. The function space is the span of .
Huybrechs showed that is not a basis of but a tight frame with frame bound , that is, the frame obeys a generalized Parseval’s identity, and that the set consists of all eigenfunctions of the Laplace operator on subject to either homogeneous Dirichlet or Neumann boundary conditions. That is due to the fact that the functions in are linearly dependent. On the contrary, the set is a basis for a finite dimensional subspace of , for any finite . All the functions in are linearly independent.
It makes perfect sense to look for an orthonormal basis on the interval . Since the even functions in and the odd functions in are mutually orthogonal, the orthonormalization problem divides into two problems. Let us define two spaces, denoted by The following two theorems are stated and proved by Huybrechs in [6].
Theorem 3. Let be the unique normalized sequence of orthogonal polynomials satisfying Then the set is an orthonormal basis for on .
Theorem 4. Let be the unique normalized sequence of orthogonal polynomials satisfying Then the set is an orthonormal basis for on .
The polynomials and are called half-range Chebyshev polynomials of the first and second kind, respectively. They have the same weight functions as the Chebyshev polynomials of the first and second kind but are orthogonal on the interval rather than on the interval . The orthogonal polynomials are guaranteed to exist, because the weight functions are positive and integrable. The construction and additional properties of these polynomials were studied by Orel and Perne in [9]. An arbitrary function can be then expanded as a half-range Chebyshev-Fourier series:
Huybrechs in [6] proved the existence and uniqueness of the exact solution to Problem 1 as a truncated half-range Chebyshev-Fourier (HCF) series.
Theorem 5. For a given , the solution to Problem 1 is where
Convergence of the HCF series (11) was extensively studied by Huybrechs in [6], where the following theorem and its corollary were stated and proved.
Theorem 6. Let , where , be analytic in the region bounded by an ellipse with major semiaxis length and with foci and . The corresponding domain of analyticity of is denoted by . If is analytic in the domain , with , then the solution to the Problem 1 satisfies with unless is analytic and periodic on .
Corollary 7. Under the conditions of Theorem 6, the coefficients and of in the form of (12) and (13) satisfy .
3. Construction of Chebyshev-Fourier-Collocation Methods
In this section we construct and analyze a new class of spectral methods for the solution of Problem 2, which we will then call Chebyshev-Fourier-collocation (CFC) methods. The numerical solution is sought as a half-range Chebyshev-Fourier series defined in (12). The series is expanded in terms of trigonometric functions and rearranged in terms of half-range Chebyshev polynomials defined in Theorems 3 and 4. Collocation is used for the computation of expansion coefficients defined in (13). Let us turn back to Problem 2 and rewrite it as a two-point boundary value problem (3) with Dirichlet boundary conditions (4): The numerical solution of the above problem is sought as a truncated HCF series, introduced in Theorem 5: In other words, we seek for coefficients and , so that the numerical solution (17) is as good as possible. Here, is the truncation number and for a given value of we use orthogonal polynomials: half-range Chebyshev polynomials of the first and half-range Chebyshev polynomials of the second kind.
We commence, by dividing the interval with collocation points, the Chebyshev points of the second kind into subintervals: where
Besides, we compute the first derivative of the truncated series (17), where and denote the coefficients of the first derivative of the truncated HCF series:
We obtain a similar form (20) for the first derivative as for the approximation of the solution (17). The coefficients and are linearly dependent on the coefficients and and are computed via the differentiation matrix : where and . Since the coefficients depend only on and only on , the differentiation matrix is block-antidiagonal: where and . This matrix was constructed and studied in detail by Orel and Perne in [9]. Furthermore, we proceed by computing the second derivative to obtain a similar representation: where the coefficients of the second derivative of the truncated series are denoted by and . Again, we compute these coefficients using differentiation matrix : where .
Solving BVPs via collocation assumes that the numerical solution of the BVP exactly solves the BVP in the interior collocation points (19) , . After inserting the truncated series (17), (20), and (23) into the differential equation (3) we obtain a system of linear equations: Besides, we have two additional equations, originating in boundary value conditions:
Let us denote by the collocation matrix. The entries of this matrix are the values of the basis functions computed in the collocation points (19). If the set of all basis functions is denoted by , then the entries of the matrix are In our case the set of basis function is composed by two sets, (half-)sines and (half-)cosines, rearranged as half-range Chebyshev polynomials of the first and second kind: and , respectively.
We have already denoted by the set of sought expansion coefficients and , arranged as a vector. Besides, we denote by the vector of values of the right-hand side function at interior Chebyshev collocation points . The first and last elements of vector are Dirichlet boundary conditions at the endpoints of the interval .
Since the coefficients of the differential equation are not necessarily constant, but functions of the independent variable , it is in general necessary to expand these coefficient functions into truncated HCF series; for example, for , we have Now we need to multiply the truncated HCF series , , and with the truncated HCF series (17) of the numerical solution and its derivatives (20) and (23). In order to perform these operations, we introduce the multiplication matrices , , and , which were constructed and studied in detail by Orel and Perne in [9]. As before, for , the multiplication matrix is a transformation matrix between the coefficients and of the truncated HCF series (17) and the coefficients and of the truncated multiplication of the truncated HCF series (28) and (17) using the coefficients and of the truncated HCF series (28): where denotes the vector of the expansion coefficients of the truncated HCF series of the multiplication . The multiplication matrix is a block matrix: where , , , and . In the case of a differential equation with constant coefficients, the matrix is scalar; otherwise the matrix is dense.
Let us now denote by the differential operator matrix where is the differentiation matrix (21) and , , and are multiplication matrices (29). We are now ready to rewrite the linear system of (25) and (26) into a matrix form where the matrix is obtained by replacing the first and the last row of the matrix with the first and the last row of the collocation matrix to satisfy the boundary conditions (26). The solution of the system (32) gives the spectral coefficients of the truncated HCF series for the numerical solution of the two-point boundary value problem (16). Multiplication with collocation matrix yields the solution values at collocation (Chebyshev) points (19).
4. Error Analysis
In this section we focus on error analysis and the rate of convergence of the new class of Chebyshev-Fourier-collocation spectral methods introduced in Section 3. In order to prove convergence of the method and to estimate the error of the numerical solution, we follow error estimation techniques presented by Canuto et al. in [13] and by Shen et al. in [14]. Since the HCF series is a generalized trigonometric series reorganized in terms of half-range Chebyshev polynomials of the first and second kind, it is convenient to follow the steps for error estimating of Fourier spectral methods.
Let us restrict to the problem of solving linear BVPs of second order with homogeneous Dirichlet boundary conditions on the interval : where is a linear differential operator defined in (5). Let us assume that the coefficient function on the interval and let us further assume that the coefficient function is differentiable and both functions and are bounded and strictly positive on the interval . Finally, let us assume that the condition holds for every .
The required Hilbert space is , the space of all square integrable functions on the interval . In this space the operator is unbounded. We denote by the appropriate inner product in and by the associated norm. Moreover, let denote the space of all trigonometric polynomials of degree that satisfy the boundary conditions (34) and let us further denote by the appropriate discrete inner product with the associated norm . The collocation solution of (33) and (34) satisfies the equations Here, the nodes , , are the interior collocation points defined in (19) and the operator is an approximation to the operator , obtained by replacing exact derivatives by interpolation derivatives (21) and (24) and expanding coefficient functions in HCF series.
Equations (33) and (34) can be equivalently written in a weak form as a bilinear form: where satisfies the boundary conditions (34). The collocation method (36) can be then rewritten as where . Analysis of convergence properties requires the existence of a dense Hilbert subspace of . An appropriate choice in our analysis is the Sobolev space where the derivative in the sense of distributions belongs to . In the following we denote the first derivative by and the second derivative by . This Sobolev space is equipped with the Sobolev norm such that for all . Note that for every , the space is contained in . Moreover, our analysis requires that the operator , or more exactly the bilinear form , satisfies the coercivity condition and the continuity condition Since where we use integration by parts in the second row, the coercivity condition (41) for the bilinear form (37) is satisfied with Similarly, since where we use integration by parts in the second row, the Cauchy-Schwartz inequality in the third row, and inequalities and in the fourth row, the continuity condition (42) is satisfied with Both constants and are independent of . Let us note that the coercivity and continuity conditions are sufficient but not necessary as is shown in the forthcoming examples.
Let us further denote , where is the exact, is the numerical solution of (33) and (34), and is a projection operator from to . Under the assumptions of the first Strang lemma (see Theorem 1.3, page 14, [14]), that is, the coercivity (41) and the continuity (42) condition are satisfied with constants and defined in (44) and (46), the problem (38) admits a unique numerical solution , satisfying Moreover, the first Strang lemma states that if the coercivity and continuity conditions are fulfilled, the error estimate of the numerical solution reads as follows:
Here, is a projection operator from to according to the discrete inner product and is then a trigonometric polynomial of degree , matching at the interior collocation points (19) and vanishing at the boundary points. The method is convergent, if all three parts of the inequality (48) converge to with . Since (see (), page 294, [13]) where is the order of smoothness of the right-hand side function and the solution , and the above is true for the first two parts of inequality (48). For the third part of that inequality we use the fact that for every , where is the Sobolev space equipped with the norm Moreover, we observe that the differential operator is bounded in this space; see Leoni [17] or Ziemer [18]. Finally, we obtain the estimate All constants , , and are independent of and . We have proved the following theorem which states the error estimate and the rate of convergence for solutions and right-hand side functions that are continuously differentiable to a certain order .
Theorem 8. Let be the exact solution of the problem (33) with boundary conditions (34) and let be the numerical solution obtained by the class of Chebyshev-Fourier-collocation (CFC) methods, constructed in Section 3. Let one assume that the solution and the right-hand side function are -times continuously differentiable. Then the estimated error of the approximation of the solution for the constructed class of CFC methods is where the constants and are defined in (44) and (46).
If the solution and the right-hand side function are smooth or analytic functions in some domain, containing the interval , then under the conditions of Theorem 6, the Chebyshev-Fourier-collocation method has a spectral rate of convergence.
Theorem 9. Let be the exact solution of the problem (33) with boundary conditions (34) and let be the numerical solution obtained by the class of Chebyshev-Fourier-collocation (CFC) methods, constructed in Section 3. Let one assume that the solution and the right-hand side function are analytic functions in the domain , defined in Theorem 6. Then the estimated error of the approximation of the solution for the constructed class of CFC methods is where is defined in (15).
5. Numerical Examples
In the following examples we compare numerical solutions obtained with two different classes of collocation spectral methods. One is a standard Chebyshev-collocation approach, where we approximate the solution using Chebyshev series; the other one is the new class of methods constructed in Section 3, where approximation of the solution with half-range Chebyshev-Fourier series is used. Much about classical Chebyshev polynomials can be found, for example, in [19]. Both methods were implemented in MATLAB by the authors.
We observe that the performance of the Chebyshev-Fourier-collocation (CFC) method is comparable with the standard Chebyshev-collocation (CC) method. However, since the absence of a fast technique for the computation of the expansion coefficients, comparable with the FFT, the performance in terms of computational costs is worse for the new class of methods. Throughout this section we use the abbreviations and .
Example 1. As a first example we consider a second order linear differential equation with nonconstant coefficients: The exact solution is
In Figure 1(a), the comparison of the decay of the maximal absolute error of the two numerical solutions, obtained with CC and CFC methods with respect to the truncation number of terms in the Chebyshev and Chebyshev-Fourier expansions, is depicted. For both methods we have to compute coefficients. The exact solution of this problem is a smooth and analytic function and both methods reach machine accuracy. However, as seen from the figure, the CC numerical solution converges more rapidly, since the maximal absolute error reaches machine accuracy at , instead of for the CFC numerical solution. Note the exponential decay of the maximal absolute error for both methods. In this case we obtain spectral accuracy as predicted by Theorem 9 in (55).
| (a) Maximal error dependent on |
| (b) Maximal error dependent on |
| (c) Maximal error dependent on |
| (d) Maximal error dependent on |
Example 2. As a second example we consider the Airy differential equation The exact solution is a highly oscillatory function where and are Airy functions; see, for example, [19].
Figure 1(b) shows the comparison of the decay of the maximal absolute error of the two numerical solutions obtained with CC and CFC methods with respect to the truncation number . The exact solution is a smooth but oscillatory function. For this reason the maximal absolute error begins to decay not until is big enough to overcome problems with the resolution. However, both methods again yield exponential decay of the maximal absolute error with respect to , which gives spectral accuracy. We note that the convergence of the Chebyshev-Fourier-collocation method is faster in comparison with the Chebyshev-collocation method.
Example 3. As a third example we consider a second order linear differential equation with nonconstant and nonsmooth coefficients: The exact solution is
Figure 1(c) shows the comparison of the decay of the maximal absolute error for the two numerical solutions with CC and CFC methods with respect to the truncation number . In this example, the exact solution is nonsmooth, actually it is only six-times continuously differentiable. For this reason the decay of the maximal absolute error is slow for both methods; however, it is faster for the CFC method. The error decays as . This is in accordance with the result (54) of Theorem 8. We note that in this case the Chebyshev-Fourier-collocation method outweighs the Chebyshev-collocation method.
Example 4. As a last example we again consider a second order linear differential equation with nonconstant and nonsmooth coefficients: The exact solution is
Figure 1(d) shows the comparison of the decay of the maximal absolute error for the two numerical solutions with CC and CFC methods with respect to the truncation number . In this example, the exact solution is again nonsmooth; actually it is only eight times continuously differentiable. For this reason the decay of the maximal absolute error is slow for both methods; however, it is faster than in the previous example; actually the error decays as . This is again in accordance with the result (54) of Theorem 8. We note that in this case the Chebyshev-Fourier-collocation method again outweighs the Chebyshev-collocation method.
6. Conclusions
Based on the papers by Huybrechs [6] and by Orel and Perne [9], we first discussed two nonclassical families of orthogonal polynomials: the half-range Chebyshev polynomials of the first and second kind and the associated half-range Chebyshev-Fourier series. We have briefly recalled the solution of the approximation problem stated in Problem 1, using the HCF series, and analyzed by Huybrechs in [6]. The main focus of the paper is on analyzing Problem 2 via collocation spectral methods using truncated HCF series. These series are generalized Fourier series, where the basic trigonometric functions are rearranged in terms of the half-range Chebyshev polynomials of the first and second kind. The usage of such reorganization which yields orthogonality of the expansion basis functions allows solving nonperiodic problems with tools, otherwise reserved for periodic problems.
In the paper we have constructed a new class of spectral methods, based on collocation, the Chebyshev-Fourier-collocation methods. The idea is similar to the widely used Chebyshev spectral methods, where instead of using classical Chebyshev series for the approximation of the solution, we expand the numerical solution as a half-range Chebyshev-Fourier series. We deal with linear two-point boundary value problems of the second order. Generalization to higher order and/or different intervals is also possible. Error analysis and convergence theory, presented in Section 4, show that for problems that are smooth enough, that is, analytic functions, we obtain spectral accuracy, that is, the maximal absolute error depending on the truncation number decays exponentially. Otherwise, we obtain algebraic convergence. These results are comparable with the ones for the standard Chebyshev-collocation methods. Examples shown at the end of the paper demonstrate exponential convergence for smooth and analytic functions and comparability with standard classes of spectral methods. Furthermore, for some problems the CFC method outweighs the CC method.
As far as we are interested in convergence theory, things are more or less comparable with standard classes of spectral methods. This is regrettably not the case if we are concerned with the computational costs. For computing spectral coefficients for the HCF series we do not yet have such a marvelous tool as it is the FFT for Fourier or Chebyshev series. This is an open problem that needs further investigation. Moreover, as a part of future work with HCF series, we are interested in evolutive time-dependent partial differential equations, for example, heat or wave equations. More of this theme will be discussed in a subsequent paper.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.