Solving Higher-Order Boundary and Initial Value Problems via Chebyshev–Spectral Method: Application in Elastic Foundation

In this work, we introduce an efficient scheme for the numerical solution of some Boundary and Initial Value Problems (BVPs-IVPs). By using an operational matrix, which was obtained from the first kind of Chebyshev polynomials, we construct the algebraic equivalent representation of the problem. We will show that this representation of BVPs and IVPs can be represented by a sparse matrix with sufficient precision. Sparse matrices that store data containing a large number of zero-valued elements have several advantages, such as saving a significant amount of memory and speeding up the processing of that data. In addition, we provide the convergence analysis and the error estimation of the suggested scheme. Finally, some numerical results are utilized to demonstrate the validity and applicability of the proposed technique, and also the presented algorithm is applied to solve an engineering problem which is used in a beam on elastic foundation.


Introduction
Initial value problems, boundary value problems, and other related problems have many applications in physics, chemistry, and different research areas (see, e.g., in [1][2][3][4][5][6][7][8][9][10][11] and the references given therein). Most higher order of these types of equations and problems do not have exact analytical solutions, thus approximation techniques must be applied.
In this work, we shall discuss the following BVPs.
with β k0 u(−1) + n−1 ∑ j=1 β kj u (j) (1) = l k , (k = 0, 1, 2, . . . , n − 1), (2) and also the following IVPs, with n−1 ∑ j=0 β kj u (j) (0) = l k , (k = 0, 1, 2, . . . , n − 1), where a i , β kj ∈ R are constants, (i = 0, 1, ..., n), (j, k = 0, 1, ..., n − 1), and R(x) and u(x) are given and unknown functions, respectively. Among numerical approaches, spectral methods are very successful and powerful tools for the numerical solution of differential and integral equations involving derivatives of integer and non-integer order. Effective properties that have encouraged many experts in numerical analysis to use spectral methods for different kinds of mathematical problems are spectral accuracy and the ease of applying these methods. One of the most important advantages of these methods is their high accuracy, the so-called convergence of "infinite order". Therefore, when the exact solution is infinitely differentiable, the numerical approximation converges faster than M −n , where M is the order of approximation and n > 0 is constant [12,13].
Pseudospectral, Galerkin, and Tau methods are several kinds of spectral methods which can be obtained from a weighted residual method [12]. The Tau method, which can be explained as a special case of spectral methods, was invented by Lanczos [14]. It has become increasingly popular with applications in many disciplines such as porous media, viscoelasticity, electrochemistry, and other problems where high accuracy is desired. Further detailed information of this procedure can be found in [15][16][17].
In the literature, there are several approaches to obtain numerical methods to solve boundary value problems and other related problems such as Chebyshev-type, Runge-Kutta type, Wavelet-Galerkin method, and Shannon approximation [5,6,[18][19][20][21][22][23]. Wolf [24] has been concerned with the numerical solution of non-singular integral and integro-differential equations by iteration with Chebyshev series. In [25], a Chebyshev series were used for nonlinear differential equations. Various works, such as the works in [26,27], introduced and discussed a Chebyshev technique for solving nonlinear optimal control problems. Mezzadri et al. [28] have been concerned with nonlinear programming methods for the solution of optimal control problems via a Chebyshev technique.
The proposed approach in this work is referred to as operational in the sense that it makes it possible to transform a given differential equation into an algebraic equation. An efficient algorithm of the Tau approximation based on Chebyshev polynomials is presented to numerically solve the boundary and initial value problems. One of the main advantages of this methodology is that the matrices generated of BVPs and IVPs are sparse. The proposed procedure can be very attractive to users who are concerned with memory conservation, much less computational cost and much more computational speed (see Remark 2). Furthermore, the simplicity of the scheme and low algorithm run time are among its other advantages. Ultimately, to demonstrate the validity and applicability of the method, an engineering problem of the fourth-order ODE which is arising in elastic foundation is solved.
The organization of this paper is as follows. In Section 2, we present the basic properties of Chebyshev polynomials. Converting boundary value problem (1) and initial value problem (3) to a matrix form based on the new proposed algorithm is shown and the convergence analysis for the first kind of Chebyshev expansion is given in Section 3. The last part of this work is concerned with four numerical examples to illustrate the method.
The Chebyshev polynomials T n (x) are defined on the interval [−1, 1]. In order to use these polynomials on the interval x ∈ [a, b], we introduce x = b−a 2 t + b+a 2 .

Outline of the Method for Boundary and Initial Value Problems
Here, Chebyshev-Tau method is expressed to solve boundary value problem (1) and initial value problem (3). The first part is concerned with the matrix structure produced by proposed method, the convergence, and the error analysis for the mentioned polynomials. The second part will be discussed with the numerical solvability of the system of algebraic equations arising in the presented method.

Operational Matrix of Derivative
Let functions u(x) and Du(x) be expanded by the Chebyshev polynomials in [−1, 1] as where V = (ν 0 , ν 1 , ν 2 , . . .), W 1 = (ω 1,0 , ω 1,1 , ω 1,2 , . . .) and T x = (T 0 (x), T 1 (x), T 2 (x), . . .) . Moreover, for any integrable function u(x) on [−1, 1], we define ν n as follows, Taking the derivative of (10), we have Due to (11) and (12), we obtain by using Equation (7), we can rewrite the last equation as follows, or equivalently where Let us introduce matrix E as a coefficient of matrix V, then from (14) we have W T 1 = EV T . The upper triangular matrix E, via easy calculations, is as follows, then for any k ∈ N, we obtain Proof. As we pointed out previously, Equation (14) can be written to the following matrix structure, Furthermore, due to Equation (11), we have From Equation (7), we can write Given the assumption, it follows that Using the process of obtaining Equation (15), we conclude Due to (15) and the last equation, we get W T 2 = E 2 V T . Therefore, by repeating this scheme, it follows that W T k = E k V T , and we can write The following theorem is concerned with the convergence and the error analysis for the Chebyshev expansion. Theorem 3. Suppose u(x) is a function that satisfies 1 −1 |u (2) (x)| 2 dx < ∞ and |u (2) (x)| ≤ K for some constant K. Then, we have u(x) can be expanded as ∑ ∞ n=0 ν n T n (x) and the series converges to u(x) uniformly, in other words Moreover, for an error estimation the following inequality holds, Proof. By showing the absolute convergence of the series ∑ ∞ n=0 |ν n |, it can be concluded that the series ∑ ∞ n=0 ν n T n (x) converges to u(x) uniformly. According to definition of ν n we have u (cos θ) sin(n + 1)θ sin θdθ , thus for n > 1, (for n = 0, c n = 1), we have Furthermore, due to assumption and by simple computation, we can obtain Finally, the last inequality shows that the series ∑ ∞ n=0 |ν n | ia absolutely convergent. Let us set due to orthogonality of T n (x) and value of |ν n |, we have:

Description of the Proposed Method for BVPs and IVPs
The main object of this section is applying the obtained consequence for constructing the Tau approximate solution with Chebyshev polynomials of the boundary value problem (1).
We suppose that where Using the above relations, (1) can be written as According to the linear independence of Chebyshev polynomials, we have (a n E n N + a n−1 E n−1 N + . . .
Let us set X = a n E n N + a n−1 E n−1 N + . . . + a 1 E N + a 0 I N+1 , similarly, for the Equation (2) Due to (22) and (23), the desired matrix can be written as where X n is obtained by eliminating the last n row of matrix X and By solving a sparse system of above algebraic equation, we can obtained unknown v i , (i = 0, 1, . . . , N) and finally, u N (x), the desired approximation, can be computed from the following relation, The following Algorithm 1 summarizes our proposed method.

X n based on X from (22).
Step 5. Set: Step 6. Obtain the solution V N from the system of (24) and set:

Construction of the Proposed Method for IVPs
Here, we consider the initial value problems (3) as It is clear that the given scheme for BVPs (1) and Algorithm 1, with small changes, can be easily applied to these problems, thus we refrain from going into details. Only Step 5 in Algorithm 1 will be changed, and we have therefore, we can obtain V N from the system of (24) and finally the desired solution can be achieved from u N(x) = T N,x V T N .

Remark 2.
As a consequence of the described method, a significant advantage for increasing the computation speed is to the sparsity of matrices E N , E 2 N ,...,E n N . If we consider matrix E N with an odd dimension, then the number of non-zero matrix elements is N(N + 2) 4 . Moreover, the number of non-zero matrix elements for matrix E N with an even dimension is [ N + 1 2 ] 2 . Furthermore, the number of non-zero matrix elements for matrix E 2 N with an even and odd dimension are , respectively and so it goes. Figures 1-3 indicate that as the derivative order rise, the sparsity of matrices and computational speed also increase.

Numerical Results
In this section, the numerical results of four test problems are reported. We solve the examples using proposed Tau approximation method based on Chebyshev basis functions. In the tables presented here, the maximal differences between exact and approximation solution has been shown by "Maximal Errors". We compare our obtained results with in the some existing numerical methods. The numerical results of the given method have been shown in Figures 4-7.
with the conditions u(−1) = 1 e 1 + We take N = 7, for computational details of the described technique in Section 3. By applying mentioned algorithm, the following matrices will be obtained.
The left and right hand-side of (24) are, respectively, By solving the linear system (24), we obtain vector V 7 as follows, By putting these values in Equation (18), an approximate solution will be obtained. The numerical results are given in Table 1. Example 2. Consider the first-order initial value problem: The exact solution of this problem under the initial condition u(0) = 0 is given by This problem was also solved in [4,21] by a method based on Legendre and Chebyshev polynomials. The best reported maximum error is O(10 −4 ). Table 2 shows that we can attain good numerical results compare to the numerical results in [4,21] for n ≥ 10. Consider an engineering problem of the fourth-order initial value problem which is arising in elastic foundation [20] as follows, ).
To obtain the approximate solution of Examples 3 and 4, we apply the proposed Tau method. The obtained numerical results in Table 3, show that the desired accuracy is obtained.

Conclusions
In this paper, an efficient and accurate numerical algorithm has been applied to solve boundary and initial value problems by using the Chebyshev-Tau method. Due to the wide application of boundary value problems in many different fields, an engineering problem which is arising in elastic foundation is chosen as test example. By utilizing proposed method, the sparsity of the obtained derivative operational matrix makes us able to solve the linear system with much less computational cost and much more computational speed. The reported results indicated that the proposed scheme can obtain appropriate approximate solutions in comparing with the exact solution. Moreover, by increasing the order of Chebyshev polynomial, a better approximation was achieved.