# Chebyshev-Fourier spectral methods for nonperiodic boundary value problems.

1. Introduction and FormulationA 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 T > 1, let [G.sub.n] be the space of 2T-periodic functions of the form

g [member of] [G.sub.n]: g(x) = [a.sub.0]/2 + [n.summation over (k=1)] ([a.sub.k] cos [pi]kx/T + [b.sub.k] sin [pi]kx/T). (1)

The Fourier extension of f [member of] [L.sup.2](-1,1), defined on the interval [-1,1] to the interval [-T, T], is the solution to the optimization problem

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (2)

In his paper, Huybrechs considered a square-integrable function f [member of] [L.sup.2](-1,1) that is not necessarily smooth or periodic. The idea to obtain a spectrally accurate Fourier series is to extend the given function f to a function g that is periodic on a wider interval [-T,T], T > 1. The Fourier series of the constructed function is obviously point-wise convergent to f on the interval [-1,1].

Huybrechs analyzed this problem for the choice T = 2 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

Ly(x) = f(x), x [member of] [-1,1], (3)

with boundary conditions

By(x) = 0, x [member of] {-1,1}, (4)

where L is a linear differential operator:

L = a[alpha](x) [d.sup.2]/d[x.sup.2] + [beta](x) d/dx + [gamma](x)I, (5)

I is the identity operator, and B 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 [-1, 1] is only a matter of simplification, since it is well known how to map an arbitrary interval [a, b] 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 [-1, 1] to an interval [-T, T] 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 T = 2 and proposed to use the set of basis functions

[D.sub.n] := [C.sub.n] [union] [S.sub.n], (6)

where

[C.sub.n] = {1/[square root of 2]} [union] [{cos [pi]kx/2}.sup.n.sub.k=1], [S.sub.n] = [{sin [pi]kx/2}.sup.2.sub.k=1]. (7)

Note that the set [C.sub.n] consists of even and the set [S.sub.n] of odd functions. The function space is the span of [D.sub.n].

Huybrechs showed that [D.sub.[infinity]] is not a basis of [L.sup.2](-1, 1) but a tight frame with frame bound 2, that is, the frame obeys a generalized Parseval's identity, and that the set Dm consists of all eigenfunctions of the Laplace operator on [-1,1] subject to either homogeneous Dirichlet or Neumann boundary conditions. That is due to the fact that the functions in [D.sub.[infinity]] are linearly dependent. On the contrary, the set [D.sub.n] is a basis for a finite dimensional subspace of [L.sup.2](-1,1), for any finite n. All the functions in [D.sub.n] are linearly independent.

It makes perfect sense to look for an orthonormal basis on the interval [-1,1]. Since the even functions in [C.sub.n] and the odd functions in [S.sub.n] are mutually orthogonal, the orthonormalization problem divides into two problems. Let us define two spaces, denoted by

[C.sub.n] := span [C.sub.n], [S.sub.n] := span [S.sub.n]. (8)

The following two theorems are stated and proved by Huybrechs in [6].

Theorem 3. Let [T.sup.h.sub.k](y) be the unique normalized sequence of orthogonal polynomials satisfying

4/[pi] [[integral].sup.1.sub.0] [T.sup.h.sub.k](y)[y.sup.e]/[square root of 1 - [y.sup.2]] dy = [[delta].sub.k,l], l = 0, ..., k - 1. (9)

Then the set [{[T.sup.k.sub.0](cos([pi]x/2))}.sup.n.sub.k=0] is an orthonormal basis for [C.sub.n] on [-1, 1].

Theorem 4. Let [U.sup.h.sub.k](y) be the unique normalized sequence of orthogonal polynomials satisfying

4/[pi] [[integral].sup.1.sub.0] [U.sup.k.sub.k](y)[y.sup.e] [square root of 1 - [y.sup.2]]dy = [[delta].sub.k,l], l = 0, ..., k - 1. (10)

Then the set [{[U.sup.h.sub.k](cos([pi]x/2)) sin([pi]x/2)}.sup.n-1.sub.k=0] is an orthonormal basis for [S.sub.n] on [-1,1],

The polynomials [T.sup.h.sub.n](x) and [U.sup.h.sub.n](x) 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 [0,1] rather than on the interval [-1,1]. 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 f [member of] [L.sup.2](-1,1) can be then expanded as a half-range Chebyshev-Fourier series:

f(x) = [[infinity].summation over (k=0)] [a.sub.k][T.sup.h.sub.k](cos [pi]x/2) + [[infinity].summation over (k=0)][b.sub.k][U.sup.h.sub.k](cos [pi]x/2) sin [pi]x/2. (11)

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 f [member of] [L.sup.2](-1,1), the solution to Problem 1 is

[g.sub.n](x) = [n.summation over (k=0)][a.sub.k][T.sup.h.sub.k](cos [pi]x/2) + [n-1.summation over (k=0)][b.sub.k][U.sup.h.sub.k](cos [pi]x/2) sin [pi]x/2, (12)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (13)

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 [??](y) = f(x), where x = (2/[pi]) [cos.sup.-1] y, be analytic in the region bounded by an ellipse with major semiaxis length R and with foci 0 and 1, The corresponding domain of analyticity of f is denoted by D(R), If f is analytic in the domain D(R), with R > 1/2, then the solution [g.sub.n] to the Problem 1 satisfies

[parallel]f - [g.sub.n][parallel] ~ [[rho].sup.-n], (14)

with

[rho] = min (3 + 2[square root of 2], 2R + [square root of 4[R.sup.2] - 1]), (15)

unless f is analytic and periodic on [-2,2],

Corollary 7. Under the conditions of Theorem 6, the coefficients [a.sub.k] and [b.sub.k] of [g.sub.n] in the form of (12) and (13) satisfy [a.sub.k], [b.sub.k] ~ [[rho].sup.-n],

3. Construction of Chebyshev-FourierCollocation 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):

[alpha](x) [d.sup.2]/d[x.sup.2+] + [beta](x) dy/dx + [gamma](x)y = f(x), x [member of][-1, 1], y(-1) = A, y(1) = B. (16)

The numerical solution of the above problem is sought as a truncated HCF series, introduced in Theorem 5:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (17)

In other words, we seek for coefficients [a.sub.k] and [b.sub.k], so that the numerical solution (17) is as good as possible. Here, N is the truncation number and for a given value of N we use 2N + 1 orthogonal polynomials: N + 1 half-range Chebyshev polynomials of the first and N half-range Chebyshev polynomials of the second kind.

We commence, by dividing the interval [-1,1] with 2N+1 collocation points, the Chebyshev points of the second kind into 2 N subintervals:

-1 = [x.sub.0] < [x.sub.1] < ... < [x.sub.2N] = 1, (18)

where

[x.sub.i] = -cos ([pi]i/2N), i = 0,1, ..., 2N. (19)

Besides, we compute the first derivative of the truncated series (17), where [a'.sub.k] and [b'.sub.k] denote the coefficients of the first derivative of the truncated HCF series:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (20)

We obtain a similar form (20) for the first derivative as for the approximation of the solution (17). The coefficients [a'.sub.k] and [b'.sub.k] are linearly dependent on the coefficients [a.sub.k] and [b.sub.k] and are computed via the differentiation matrix D:

u' = Du, (21)

where u = [([a.sub.0], ..., [a.sub.N], [b.sub.0], ..., [b.sub.N-1]).sup.T] and u' = ([a'.sub.0], ..., [a'.sub.N], [b'.sub.0], ..., [b'.sub.N-1])T. Since the coefficients [a'.sub.j] depend only on [b.sub.k] and [b'.sub.j] only on [a.sub.k], the differentiation matrix D [member of] [R.sup.(2N+1)x(2N+1)] is blo ck-antidiagonal:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (22)

where [H.sub.1] [member of] [R.sup.(N+1)xN] and [H.sub.2] [member of] [R.sup.Nx(N+1)]. 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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (23)

where the coefficients of the second derivative of the truncated series are denoted by [a".sub.k] and [b".sub.k]. Again, we compute these coefficients using differentiation matrix D:

u" = [D.sup.2]u, (24)

where u" = [([a".sub.0], ..., [a".sub.N], [b".sub.0], ..., [b".sub.N-1]).sup.T].

Solving BVPs via collocation assumes that the numerical solution of the BVP exactly solves the BVP in the interior collocation points (19) [x.sub.i], i = 1,2, ..., 2N- 1. After inserting the truncated series (17), (20), and (23) into the differential equation (3) we obtain a system of linear equations:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (25)

Besides, we have two additional equations, originating in boundary value conditions:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (26)

Let us denote by C [member of] [R.sup.(2N+1)x(2N+1)] 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 [{[[phi].sub.j]}.sup.2N.sub.j=0], then the entries of the matrix C = [[c.sub.ij]] are

[c.sub.ij] = [[phi].sub.j]([x.sub.i]). (27)

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: [T.sup.h.sub.k](cos([pi]x/2)) and [U.sup.h.sub.k](cos([pi]x/2)) sin([pi]x/2), respectively.

We have already denoted by u the set of sought expansion coefficients [a.sub.k] and [b.sub.k], arranged as a vector. Besides, we denote by v = [(A, f([x.sub.1]), ..., f([x.sub.2N-1]), B).sup.T] the vector of values of the right-hand side function f at interior Chebyshev collocation points [x.sub.1], [x.sub.2], ..., [x.sub.2N-1]. The first and last elements of vector v are Dirichlet boundary conditions at the endpoints of the interval [-1,1].

Since the coefficients of the differential equation are not necessarily constant, but functions of the independent variable x, it is in general necessary to expand these coefficient functions into truncated HCF series; for example, for [eta] [member of] [[alpha], [beta], [gamma]}, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (28)

Now we need to multiply the truncated HCF series [[alpha].sub.N], [[beta].sub.N], and [[gamma].sub.N] with the truncated HCF series [[gamma].sub.N] (17) of the numerical solution and its derivatives [y'.sub.N] (20) and [y".sub.N] (23). In order to perform these operations, we introduce the multiplication matrices [F.sub.[alpha]], [F.sub.[beta]], and [F.sub.[gamma]], which were constructed and studied in detail by Orel and Perne in [9]. As before, for [eta] [member of] [[alpha], [beta], [gamma]}, the multiplication matrix [F.sub.[eta]] [member of] [R.sup.(2N+1)x(2N+1)] is a transformation matrix between the coefficients [a.sub.k] and [b.sub.k] of the truncated HCF series (17) and the coefficients [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] of the truncated multiplication of the truncated HCF series (28) and (17) using the coefficients [[mu].sub.k] and [v.sub.k] of the truncated HCF series (28):

[??] = [F.sub.[eta]]u, (29)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] denotes the vector of the expansion coefficients of the truncated HCF series of the multiplication [[eta].sub.N](x) x [y.sub.N](x). The multiplication matrix [F.sub.[eta]] [member of] [R.sup.(2n+1)x(2n+1)] is a block matrix:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (30)

where [G.sub.1] [member of] [R.sup.(N+1)x(N+1)], [G.sub.2] [member of] [R.sup.(N+1)xN], [G.sub.3] [member of] [R.sup.Nx(N+1)], and [G.sub.4] [member of] [R.sup.NxN]. In the case of a differential equation with constant coefficients, the matrix [F.sub.[eta]] is scalar; otherwise the matrix is dense.

Let us now denote by L [member of] [R.sup.(2N+1)x(2N+1)] the differential operator matrix

L = [F.sub.[alpha]][D.sup.2] + [F.sub.[beta]]D + [F.sub.[gamma]], (31)

where D is the differentiation matrix (21) and [F.sub.[alpha]], [F.sub.[beta]], and [F.sub.[gamma]] are multiplication matrices (29). We are now ready to rewrite the linear system of (25) and 26) into a matrix form

Uu = v, (32)

where the matrix U [member of] [R.sup.(2N+1)x(2N+1)] is obtained by replacing the first and the last row of the matrix CL with the first and the last row of the collocation matrix C 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 C 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 [-1,1]:

Lu = f, (33)

u (-1) = u (1) = 0, (34)

where L is a linear differential operator defined in (5). Let us assume that the coefficient function [alpha](v) [equivalent to] -1 on the interval [-1, 1] and let us further assume that the coefficient function ft is differentiable and both functions [beta] and [gamma] are bounded and strictly positive on the interval [-1,1]. Finally, let us assume that the condition [gamma](v) - [beta]'(x)/2 > 0 holds for every x [member of] [-1,1].

The required Hilbert space is [L.sup.2](-1,1), the space of all square integrable functions on the interval [-1,1]. In this space the operator L is unbounded. We denote by

(u, v) = [[integral].sup.1.sub.-1]u(x)v(x)dx (35)

the appropriate inner product in [L.sup.2](-1,1) and by [parallel]u[parallel] = [(u,u).sup.1/2] the associated norm. Moreover, let [T.sub.N] [subset] [L.sup.2](-1,1) denote the space of all trigonometric polynomials of degree [less than or equal to] N that satisfy the boundary conditions (34) and let us further denote by [(u, v).sub.N] the appropriate discrete inner product with the associated norm [parallel]u[[parallel].sub.N] = [(u, u).sup.1/2.sub.N]. The collocation solution [u.sup.N] [member of] [T.sub.N] of (33) and 34) satisfies the equations

[L.sub.N] [u.sup.N] ([x.sub.k]) = f ([x.sub.k]), (36)

[u.sup.N]([x.sub.0]) = [u.sup.N]([x.sub.2N]) = 0.

Here, the nodes [x.sub.k], k = 1,2, ..., 2N - 1, are the interior collocation points defined in (19) and the operator [L.sub.N] is an approximation to the operator L, 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:

(Lu, v) = (f, v), [for all]v [member of] [L.sup.2](-1,1), (37)

where u satisfies the boundary conditions (34). The collocation method (36) can be then rewritten as

[([L.sub.N][u.sup.N], v).sub.N] = [(f, v).sub.N], v [member of] [T.sub.N], (38)

where [u.sup.N] [member of] [T.sub.N]. Analysis of convergence properties requires the existence of a dense Hilbert subspace of [L.sup.2](-1,1). An appropriate choice in our analysis is the Sobolev space

[H.sup.1](-1,1) = {v [member of] [L.sup.2](-1,1); dv/dx [member of] [L.sup.2](-1,1)}, (39)

where the derivative dv/dv in the sense of distributions belongs to [L.sup.2](-1,1). In the following we denote the first derivative by v' [equivalent to] dv/dv and the second derivative by v" [equivalent to] [d.sup.2]v/d[x.sup.2]. This Sobolev space is equipped with the Sobolev norm

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (40)

such that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Note that for every N > 0, the space [T.sub.N] is contained in [H.sup.1](-1,1). Moreover, our analysis requires that the operator L, or more exactly the bilinear form (Lu, v), satisfies the coercivity condition

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (41)

and the continuity condition

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (42)

Since

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (43)

where we use integration by parts in the second row, the coercivity condition (41) for the bilinear form (37) is satisfied with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (44)

Similarly, since

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (45)

where we use integration by parts in the second row, the Cauchy-Schwartz inequality in the third row, and inequalities [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] in the fourth row, the continuity condition (42) is satisfied with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (46)

Both constants [[alpha].sup.*] and A are independent of N. 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 e = [u.sup.N] - [R.sub.N]u, where u [member of] [L.sup.2](-1,1) is the exact, [u.sup.N] [member of] [T.sub.N] is the numerical solution of (33) and (34), and [R.sub.N] is a projection operator from [L.sup.2](-1,1) to [T.sub.N]. 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 [[alpha].sup.*] and A defined in (44) and (46), the problem (38) admits a unique numerical solution [u.sup.N] [member of] [T.sub.N], satisfying

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (47)

Moreover, the first Strang lemma states that if the coercivity and continuity conditions are fulfilled, the error estimate of the numerical solution [u.sup.N] reads as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (48)

Here, [Q.sub.N] is a projection operator from [L.sup.2](-1,1) to TN according to the discrete inner product and [Q.sub.N]v is then a trigonometric polynomial of degree N, matching v 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 0 with N [right arrow] [infinity]. Since (see (5.5.15), page 294, [13])

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (49)

where m is the order of smoothness of the right-hand side function f and the solution w, and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (50)

the above is true for the first two parts of inequality (48). For the third part of that inequality we use the fact that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for every v [member of] [H.sup.2](-1,1), where

[H.sup.2](-1,1) = {v [member of] [L.sup.2](-1,1); v', v" [member of] [L.sup.2](-1,1)}, (51)

is the Sobolev space equipped with the norm

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (52)

Moreover, we observe that the differential operator L is bounded in this space; see Leoni [17] or Ziemer 18]. Finally, we obtain the estimate

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (53)

All constants [C.sub.1], [C.sub.2], and D are independent of N and m. 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 m.

Theorem 8. Let u [member of] [L.sup.2](-1,1) be the exact solution of the problem (33) with boundary conditions (34) and let [u.sup.N] [member of] [T.sub.N] 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 m-times continuously differentiable. Then the estimated error of the approximation of the solution for the constructed class of CFC methods is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (54)

where the constants [[alpha].sup.*] and A are defined in (44) and (46).

If the solution u and the right-hand side function f are smooth or analytic functions in some domain, containing the interval [-1,1], then under the conditions of Theorem 6, the Chebyshev-Fourier-collocation method has a spectral rate of convergence.

Theorem 9. Let u [member of] [L.sup.2](-1,1) be the exact solution of the problem (33) with boundary conditions (34) and let [u.sup.N] [member of] [T.sub.N] 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 D(R), defined in Theorem 6. Then the estimated error of the approximation of the solution for the constructed class of CFC methods is

[parallel]u - [u.sup.N][parallel] ~ [[rho].sup.-N], (55)

where [rho] 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 y' [equivalent to] dy/dv and y" [equivalent to] [d.sup.2]y/d[x.sup.2].

Example 1. As a first example we consider a second order linear differential equation with nonconstant coefficients:

y" + xy' = (2 + [x.sup.2]) cos x, (56)

y (-1) = sin (1), y (1) = sin (1).

The exact solution is

y(x) = x sin x. (57)

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 N of terms in the Chebyshev and Chebyshev-Fourier expansions, is depicted. For both methods we have to compute N + 1 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 N = 14, instead of N = 34 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).

Example 2. As a second example we consider the Airy differential equation

y" - (x - 1000) y = 0, (58)

y(-1) = 1, y(i) = i

The exact solution is a highly oscillatory function

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (59)

where Ai and Bi 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 N. The exact solution is a smooth but oscillatory function. For this reason the maximal absolute error begins to decay not until N is big enough to overcome problems with the resolution. However, both methods again yield exponential decay of the maximal absolute error with respect to N, 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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (60)

The exact solution is

y(x) = [x.sup.5][absolute value of x]. (61)

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 N. 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 O(1/[N.sup.5]). 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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (62)

The exact solution is

y (x) = [x.sup.7] [absolute value of x]. (63)

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 N. 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 O(1/[N.sup.7]). 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 N 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.

http://dx.doi.org/10.1155/2014/572694

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

[1] D. Gottlieb and C.-W. Shu, "On the Gibbs phenomenon and its resolution," SIAM Review, vol. 39, no. 4, pp. 644-668, 1997.

[2] E. Tadmor, "Filters, mollifiers and the computation of the Gibbs phenomenon," Acta Numerica, vol. 16, pp. 305-378, 2007

[3] J. P. Boyd, Chebyshev and Fourier Spectral Methods, Dover Publications, Mineola, NY, USA, 2nd edition, 2001.

[4] B. Fornberg and D. M. Sloan, "A review of pseudospectral methods for solving partial differential equations," in Acta Numerica, 1994, Acta Numerica, pp. 203-267, Cambridge University Press, Cambridge, UK, 1994.

[5] L. N. Trefethen, Spectral Methods in MATLAB, vol. 10 of Software, Environments, and Tools, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, Pa, USA, 2000.

[6] D. Huybrechs, "On the Fourier extension of nonperiodic functions," SIAM Journal on Numerical Analysis, vol. 47, no. 6, pp. 4326-4355, 2010.

[7] J. P. Boyd, "A comparison of numerical algorithms for Fourier extension of the first, second, and third kinds," Journal of Computational Physics, vol. 178, no. 1, pp. 118-160, 2002.

[8] O. P. Bruno, Y. Han, and M. M. Pohlman, "Accurate, high-order representation of complex three-dimensional surfaces via Fourier continuation analysis," Journal of Computational Physics, vol. 227, no. 2, pp. 1094-1125, 2007.

[9] B. Orel and A. Perne, "Computations with half-range Chebyshev polynomials," Journal of Computational and Applied Mathematics, vol. 236, no. 7, pp. 1753-1765, 2012.

[10] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral Methods in Fluid Dynamics, Springer Series in Computational Physics, Springer, New York, NY, USA, 1988.

[11] B. Fornberg, A Practical Guide to Pseudospectral Methods, vol. 1 of Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, Cambridge, UK, 1996.

[12] D. Gottlieb and S. A. Orszag, Numerical Analysis of Spectral Methods: Theory and Applications, CBMS-NSF Regional Conference Series in Applied Mathematics, no. 26, Society for Industrial and Applied Mathematics, Philadelphia, Pa, USA, 1977.

[13] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral Methods: Fundamentals in Single Domains, Scientific Computation, Springer, Berlin, Germany, 2006.

[14] J. Shen, T. Tang, and L.-L. Wang, Spectral Methods: Algorithms, Analysis and Applications, vol. 41 of Springer Series in Computational Mathematics, Springer, Heidelberg, Germany, 2011.

[15] B. Adcock, "Spectral methods and modified Fourier series," Tech. Rep. NA2007/08, DAMTP, University of Cambridge, 2007

[16] B. Adcock, "Univariate modified Fourier methods for second order boundary value problems," BIT Numerical Mathematics, vol. 49, no. 2, pp. 249-280, 2009.

[17] G. Leoni, A First Course in Sobolev Spaces, vol. 105 of Graduate Studies in Mathematics, American Mathematical Society, Providence, RI, USA, 2009.

[18] W. P. Ziemer, Weakly Differentiable Functions: Sobolev Spaces and Functions of Bounded Variation, vol. 120 of Graduate Texts in Mathematics, Springer, New York, NY, USA, 1989.

[19] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, vol. 55 of National Bureau of Standards Applied Mathematics Series, For Sale by the Superintendent of Documents, U.S. Government Printing Office, Washington, DC, USA, 1964.

Bojan Orel (1) and Andrej Perne (2)

(1) Faculty of Computer and Information Science, University of Ljubljana, Trzaska Cesta 25,1000 Ljubljana, Slovenia

(2) Faculty of Electrical Engineering, University of Ljubljana, Trzaska Cesta 25,1000 Ljubljana, Slovenia

Correspondence should be addressed to Andrej Perne; andrej.perne@fe.uni-lj.si

Received 15 January 2014; Revised 9 May 2014; Accepted 13 May 2014; Published 1 June 2014

Academic Editor: Miroslaw Lachowicz

Printer friendly Cite/link Email Feedback | |

Title Annotation: | Research Article |
---|---|

Author: | Orel, Bojan; Perne, Andrej |

Publication: | Journal of Applied Mathematics |

Article Type: | Report |

Date: | Jan 1, 2014 |

Words: | 6680 |

Previous Article: | Effects of controller and nonuniform temperature profile on the onset of Rayleigh-Benard-Marangoni electroconvection in a micropolar fluid. |

Next Article: | Applications of Schauder's fixed point theorem to semipositone singular differential equations. |

Topics: |