Tchebychev Polynomial Approximations for $m^{th}$ Order Boundary Value Problems

Higher order boundary value problems (BVPs) play an important role modeling various scientific and engineering problems. In this article we develop an efficient numerical scheme for linear $m^{th}$ order BVPs. First we convert the higher order BVP to a first order BVP. Then we use Tchebychev orthogonal polynomials to approximate the solution of the BVP as a weighted sum of polynomials. We collocate at Tchebychev clustered grid points to generate a system of equations to approximate the weights for the polynomials. The excellency of the numerical scheme is illustrated through some examples.


Introduction
Boundary value problems are an interesting topic in the field of applied mathematics, physics and engineering [2,3,4,6,7,13]. Higher order BVPs arise in modeling various physical, chemical and biological realities. As for example, the Schrdinger equation can be presented by a second order BVP [21]; the free vibration analysis of beam structures is governed by a 4 th order differential equation; and that of ring structures by a 6 th order differential equation; an ordinary convection yields a 10 th order BVP; an ordinary overstability yields a 12 th order BVP [20]. Also an historically important example of a 4 th order BVP is the Orr Sommerfeld equation from the field of hydrodynamic stability [17]. Numerical approximations for these problems are challenging because of the higher order derivatives and boundary conditions involving higher order derivatives of the unknown function.
In this study we consider the m th order non-homogeneous linear differential equation y (m) (t)+a 1 (t)y (m−1) (t)+a 2 (t)y (m−2) (t)+· · ·+a m−1 (t)y (1) (t)+a m (t)y(t) = f (t), t ∈ (a, b), (1) with a set of m boundary conditions. Here a i (t), i = 1, 2, · · · , m and f (t) are real valued functions. Here we consider (1) with two exemplary set of boundary conditions though the scheme is not limited to them.
The existence and uniqueness of analytic solutions of the higher order BVP (1) have been well studied [3,11,12] and many references therein. The numerical approximation on the solution of the m th -order BVP (1) is sparse. So we are not interested in discussing the existence and uniqueness issues rather we focus on an efficient numerical scheme for (1).
The authors concentrate on numerical methods for 4 th order BVPs in [5,18]. A low-order numerical method is outlined in [18]. Taher et. al [16] study 4 th order Sturm-Liouville BVPs. They use Tchebychev differentiation matrices to approximate the model. They present some examples to demonstrate the scheme.
A detailed study about a 6 th order BVPs can be found in [10]. They update various existing MATLAB codes to solve the 6 th order differential equation with various boundary conditions using finite difference based schemes. A good discussion about the existence and uniqueness of solutions of such BVPs can be found in this article.
A locally supported Lagrange polynomial approximation of higher order BVPs can be found in [20]. They approximate a 6 th order differential equation and an 8 th order differential equation, respectively to illustrate the proposed method with two different set of boundary conditions.
Recently Adomian decomposition methods for 4 th order BVPs are used by [1]. Shi and Cao [22] approximates a higher order eigenvalue problem using Haar wavelets. Whereas [19] applies a differential quadrature method for a 4 th order BVP.
The famous and popular texts on higher order BVPs using spectral method are [17,21]. In both of the study the authors consider Tchebychev nodes to approximate the higher order derivatives in the MATLAB programming environment. They present numerical solutions of some 2 nd and 4 th order BVPs and eigenvalue problems. Both the authors use some differentiation matrices to approximate the higher order derivatives. Then they discuss issues regarding boundary conditions. They discuss some difficulties of handling boundary conditions, specially boundary conditions related to the higher order derivatives. The limitation of the scheme is that nonuniqueness of approximate solutions may occur [21] because of the boundary conditions. From all the discussions above and the references therein we experience that the main difficulty of using these schemes is to implement/handle boundary conditions. Studying all these articles we notice that most techniques are developed focusing some particular boundary conditions. Most of the articles are restricted to 4 th order problems. It is also noticed that some articles discuss schemes for even order BVPs and some for odd order BVPs. Thus a complete study about a general higher order BVPs is missing. Also, in most cases, the authors compute Tchebychef differentiation matrices (TDM) to approximate all higher order derivatives. The implementation of the boundary conditions in TDM is complicated [17,21]. Here we aim to develop a scheme where we do not need to compute TDMs, and handling boundary conditions is easy.
Here we propose a simple and efficient scheme for a general higher order BVP (1). This scheme works well for both the even and the odd order differential equations with any type of boundary conditions. The efficient recipe is based on the Tchebychev collocation method. In this scheme one can handle boundary conditions without any complication. Actually we convert the higher order BVP to a first order BVP. We approximate the resulting first order system of differential equations using orthogonal polynomials defined in a bounded interval [a, b]. Here all boundary conditions involving higher order derivatives have been converted to boundary conditions for the respective transformed functions. Then we present the unknown functions and the boundary values by weighted sums of truncated Tchebychev polynomials. As a result the scheme gives us a system of linear equations for the Tchebychev weights/coefficients which can be solved by using a standard linear algebra solver. We do not have to compute Tchebychev differentiation matrices for higher order derivatives here. Thus the scheme becomes simple and efficient. The proposed scheme preserves spectral accuracy which has many applicabilities in physical and engineering models. Here along with the solution of (1) the proposed scheme generates an approximation of the first m − 1 derivatives of the unknown function which is important in the study of some physical realities.
The rest of the article is organised as follows: In Section 2 we convert the higher order differential equation (1) to a set of first order differential equations. We discuss some basic results on Tchebychev polynomials and recurrence relations to support our study in Section 3. A simple scheme for the resulting system of first order differential equations has been developed in Section 4. We finish with some numerical examples and discussion in Section 5.
2 Reduction of order for the BVP (1) Any higher order differential equation can be converted into a system of first order differential equations. Thus a higher order BVP can be converted into a first order system of differential equations with the same boundary conditions. Here in this study we aim to develop a scheme to approximate higher order BVPs after converting them to first order system of differential equations. To that end using the transformations we rewrite the m th order linear BVPs (1) as a system of ordinary differential equations (3) The system of equations (3) can be written as where The boundary condition for the k th derivative of y(t) is now transformed into the boundary condition of the unknown function y k (t) for all k = 0, 1, 2, · · · , m − 1. So in the similar fashion the transformed boundary conditions can be written as or for m even y k+1 (a) = α k , k = 0, 1, 2, · · · , m/2 − 1, From now on we consider the system of first order differential equations (4) with the set of transforms boundary conditions (5) and (6).

Tchebychev polynomial approximation
Before moving onto our main study we discuss some basic properties to support our scheme. In this article we are concerned with approximating solutions of (4) using Tchebychev polynomials as trials functions and Tchebychev nodes for collocation. Polynomial interpolations based on Tchebychev nodes are often used to approximate smooth function. The pseudospectral methods perform well, in cases where both solutions and coefficients are not smooth [8,9] as well. In general, Spectral methods [17] are a class of spatial discretizations for differential equations. They can be categorised as Galerkin, tau and collocation(or pseudospectral) spectral methods. Galerkin and tau work with the coefficients of a global expansion whereas pseudospectral work with the values at collocation points. There are two key components for the formulation of Spectral methods: • Trial function, which is also called the expansion or approximation functions.
• Test function, which is also known as weight functions.
Let L n u(t) be the n th degree polynomial approximation of u(t). We state an upper bound for the solutions generated by pseudospectral approximation.
where Ω is the domain for u(t), and where φ j (t) are orthogonal polynomials defined in Ω.
Theorem 3.1. [15] Let u(t) be finite at every point of the finite interval (a, b) and such that b a u 2 (t)dt exists. Then if {φ n (t)} represents any set of orthogonal polynomials corresponding to (a, b), and ψ is a weight function.
Proof. For exact details see [15]. Now a days, a lot of attention has been grown to the study of the Tchebychev orthogonal polynomial approximation focusing various real life models. The efficiency of the method is very important, and have been well studied by several authors [17,21]. The goal of this section is to recall some properties of the Tchebychev polynomials, state some known results, and derive useful formulas that are important for this study.
The Tchebychev points are unequally spaced points over [−1, 1]. These points are the horizontal coordinates of a unit circle with center (0, 0). It is to note that they are numbered from right to left. These points are the extreme points of the nth degree Tchebychev polynomial of the first kind. These points are defined by t i = cos iπ n , i = 0, 1, 2, · · · , n.
This polynomials satisfy the following properties: Any analytic function u(t) can be approximated by a truncated Tchebychev series as which can be written as where T (t) = [T 0 (t), T 1 (t), · · · , T n (t)], and c = [c 0 , c 1 , · · · , c n ] ′ .
Here, following the above convergence result we see that if the function u belongs to C ∞ class, the produced error of approximation as n tends to infinity, approaches zero with exponential rate (O(e −c + n ), for some c + > 0). This phenomenon is usually referred to as "spectral accuracy". Assuming the function is differentiable the derivatives of u(t) can be written as [14] It is well known that the coefficients c of u (k) (t) and u (k+1) (t) can be written by the following three term recurrence relation can be simplified into the following relation [14] c (k+1) where c (0) r = c r . Accordingly, the truncated system c (k) can be computed as where when n is odd, and when n is even,

Tchebychev polynomials on [a, b]
Now the Tchebychev polynomials can be defined over an interval [a, b] by The collocation points in [a, b] can be defined by Also using the transformed Tchebychev polynomials a unknown function v * (t) and the derivatives can be approximated by

Polynomial approximation of the first order system of BVPs
In this section we motivate ourselves to solve the system of differential equations (4) with a given set of boundary conditions. To that end we recall (4) with solutions as a truncated series of Tchebychef polynomials given by which can be expressed as , and C i = [C i,0 C i,1 C i,2 · · · C i,n ] ′ ; i = 1, 2, · · · , m, and T * j (t) are given by (10) in [a b]. It is easy to see that y i (t) approximates y (i−1) (t) for i = 1, 2, 3, · · · , m. Here we aim to compute C i,j values so that (13) approximates (4) satisfying all the boundary conditions. We write the derivatives of the unknown function y i (t) as Combining (13), (14) and (12) the derivatives involved in (4) can be expressed as Thus y(t) and all the derivatives in (4) can be approximated by Using the above notations the system of ODEs (4) gives us a system of algebraic equations for the Tchebychef weights C as where To solve the system for C i,j we collocate (15) at translated Tchebychev nodes (11). Collocating at the prescribed grid points (15) yields where U (i) , i = 0, 1, are m(n + 1) column vectors, · · · f (t n )] ′ , a m(n + 1) column vector, , a m(n + 1) × m(n + 1) matrix, and , a m(n + 1) × m(n + 1) matrix. Thus (16) gives a full discrete system of equation of the form WC = F .
(17) needs to be solved for m(n + 1) unknowns C after imposing boundary conditions.
So the matrix W is nonsingular even if B is singular, since Now we discuss boundary conditions. We have m boundary conditions. Depending on the physical conditions modelled the k th left and right boundary conditions (5) or (6) (or BCs of any type) can be presented by n j=0 C k,j T * j (t n ) = α k , and n j=0 C k,j T * j (t 0 ) = β k , which can be written in vector form as So we need to solve the system of equations (17) and (18) for C i,j . Now (17) and (18) involve the Tchebychev weights C i,j for the unknown functions y i at t = t 0 and t = t n . We need to apply boundary conditions (18) in (17) before we attempt to solve them for C i,j . Here we combine (17) and (18) to solve the system for the weights.

Right Boundary Conditions
The first m rows of W, F and Ψ have been computed at t = t 0 . Let us assign them as submatrices W 1 , F 1 , and Ψ 1 . Now the k th row of W 1 , and F 1 are assigned for the unknown function y k (t), k = 1, 2, · · · , m at t = t 0 . Thus we replace the k th row of W 1 by the k th row of Ψ 1 (W 1 (k, :) = Ψ 1 (k, :)), and respective k th element of F 1 by the boundary value β k (F 1 (k, 1) = β k ) for all boundary conditions at t 0 = b.

Left Boundary Conditions
The last m rows of W, F and Ψ have been computed at t = t n . Let us assign them as submatrices W n , F n , and Ψ n . Now the k th row of W n , and F n are assigned for the unknown function y k (t), k = 1, 2, · · · , m at t = t n . Thus we replace the k th row of W n by the k th row of Ψ n (W n (k, :) = Ψ n (k, :)), and respective k th element of F n by the boundary value α k (F n (k, 1) = α k ) for all boundary conditions at t n = a.
Thus replacing the first m rows of W by the rearranged W 1 , the last m rows of W by the rearranged W n , first m elements of F by the rearranged F 1 , and the last m elements of F by the rearranged F n we get the system of equations of the formW wherẽ The m(n + 1) × m(n + 1) system of equations (19) can be solved for C by using any standard linear system solver. Here • C 1,j , j = 0, 1, 2, · · · , n are the Tchebychev weights for the solution y(t) of (1), • C 2,j , j = 0, 1, 2, · · · , n are the Tchebychev weights for y ′ (t) of (1), • C 3,j , j = 0, 1, 2, · · · , n are the Tchebychev weights for y ′′ (t) of (1), . . .
Thus we get a simple algorithm to solve the higher order linear differential equation (1) with any set of boundary conditions. In the next section we inspect some examples to show the efficiency of the scheme.

Numerical experiments and discussions
In this section we present numerical solutions of some higher order BVPs using the method outlined in the previous section. First define the set of clustered nodes (11). We use the transformation (2) to convert the higher order BVPs to a system of first order BVPs. Here we consider some examples (BVPs) with exact solutions so that we can present accuracy of our scheme umerically.
The exact solution is y(t) = 1 − cos t + cot 1 sin t.
The following Figure 1 demonstrates the numerical solution of the 2 nd order BVP. Here from the error computation we observe that the accuracy of the scheme is of an exponential order and we need 8 nodes for a 10 digit accurate solution.
The following Figure 2 demonstrates the numerical solution of the 4 th order BVP. Here from the error computation we observe that the accuracy of the scheme is of an exponential order. Here we notice that we need 9 nodes for a 10 digit accurate solution.
Example 3. We consider the following BVP The exact solution is The following Figure 3 demonstrates the numerical solution of the 5 th order BVP. Here from the error computation we observe that the accuracy of the scheme is of an exponential order. Here we need 11 nodes for a 10 digit accurate solution.   with y(0) = 1, y ′ (0) = 0, y (2) (0) = −1, y(1) = 0, y ′ (1) = −e, y (2) (1) = −2e. The exact solution is The following Figure 4 demonstrates the numerical solution of the 6 th order BVP. Here from the error computation we observe that the accuracy of the scheme is of an exponential order and 11 nodes are enough for a 10 digit accurate solution.
The following Figure 5 demonstrates the numerical solution of the 9 th order BVP. Here from the error computation we observe that the accuracy of the scheme is of an exponential order. Here we see that we need 13 nodes for a 8 digit accurate solution. We discuss in detail the formulations of a spectral collocation method for a general m th order BVPs using Tchebychev polynomials as basis functions. From these computations we notice that the result is minimum 8 − 15 digit accurate when n > 10 which agree with other published works available in the literature. Here the advantage of using this scheme is that we can handle boundary conditions simply and efficiently. It preserves spectral accuracy as one expects from polynomial approximations at the Tchebychev notes. Therefore, we conclude that the Tchebychev polynomial approach for the transformed system of first order BVPs produce much accurate results than all other previously published lower order works. From these computations we see that our method compete very well with other methods cited in this article.
There are some benefits of using the transformed scheme for (1). Here we do not need to compute higher order derivatives of the unknown function. In this presented scheme, by converting the higher order differential equation (1) to the first order system of differential equations (4) we avoid computing Tchebychev differentiation matrices (what one needs while using Tchebychev polynomials to approximate higher order derivatives for a direct scheme for (1) [17]). All boundary conditions y (i) (a) = α i and y (i) (b) = β i have been converted to boundary conditions y i (a) = α i , and y i (b) = β i respectively, which are simple to handle in our proposed set up. Whereas applying y (i) (a) and y (i) (b) in the Tchebychev differentiation matrices is really complicated (see [17,21] for the exact detail). Thus the scheme we present here is simply efficient for solving the higher order BVP (1), and preserves spectral accuracy.