Collocation Method via Jacobi Polynomials for Solving Nonlinear Ordinary Differential Equations

Abstract

We extend a collocation method for solving a nonlinear ordinary differential equation (ODE) via Jacobi polynomials. To date, researchers usually use Chebyshev or Legendre collocation method for solving problems in chemistry, physics, and so forth, see the works of (Doha and Bhrawy 2006, Guo 2000, and Guo et al. 2002). Choosing the optimal polynomial for solving every ODEs problem depends on many factors, for example, smoothing continuously and other properties of the solutions. In this paper, we show intuitionally that in some problems choosing other members of Jacobi polynomials gives better result compared to Chebyshev or Legendre polynomials.

1. Introduction

The Jacobi polynomials with respect to parameters 𝛼>−1, 𝛽>−1 (see, e.g., [1, 2]) are sequences of polynomials 𝑃_𝑛^(𝛼,𝛽)(𝑥);(𝑛=0,1,2,…) satisfying the following relation ¹₋₁(1−𝑥)^𝛼(1+𝑥)^𝛽𝑃_𝑛^(𝛼,𝛽)(𝑥)𝑃_𝑚^(𝛼,𝛽)(𝑥)ℎ𝑑𝑥=0,𝑚≠𝑛,_𝑛,𝑚=𝑛, (1.1) where ℎ_𝑛=2^𝛼+𝛽+1Γ(𝛼+𝑛+1)Γ𝛽+𝑛+12𝑛+𝛼+𝛽+1𝑛!Γ𝑛+𝛼+𝛽+1. (1.2)

These polynomials are eigenfunctions of the following singular Sturm-Liouville equation: 1−𝑥²Φ^′′(𝑥)+𝛽−𝛼−𝛽+𝛼+2𝑥Φ^′(𝑥)+𝑛𝑛+𝛽+𝛼+1Φ(𝑥)=0. (1.3)

A consequence of this is that spectral accuracy can be achieved for expansions in Jacobi polynomials so that𝑃_𝑛^(𝛼,𝛽)(1)=Γ(𝑛+𝛼+1)(),𝑃𝑛!Γ𝛼+1_{𝑛(𝛼,𝛽)}()=()−1−1^𝑛Γ𝑛+𝛽+1.𝑛!Γ𝛽+1 (1.4)

The Jacobi polynomials can be obtained from Rodrigue's formula as (1−𝑥)^𝛼(1+𝑥)^𝛽𝑃_𝑛^(𝛼,𝛽)(𝑥)=(−1)^𝑛𝑑^𝑛2^𝑛𝑛!𝑑𝑥^𝑛(1−𝑥)^𝛼+𝑛(1+𝑥)^𝛽+𝑛. (1.5) Furthermore, we have that 𝑑𝑃𝑑𝑥^{𝑛(𝛼,𝛽)}(𝑥)=12𝑛+𝛼+𝛽+1𝑃^{(𝛼+1,𝛽+1)}_𝑛−1(𝑥). (1.6)

The Jacobi polynomials are normalized such that 𝑃_𝑛^(𝛼,𝛽)(1)=𝑛=𝑛+𝛼Γ(𝑛+𝛼+1)𝑛!Γ(𝛼+1). (1.7)

An important consequence of the symmetry of weight function 𝑤(𝑥) and the orthogonality of Jacobi polynomial is the symmetric relation 𝑃_𝑛^(𝛼,𝛽)(1)=(−1)^𝑛𝑃_𝑛^(𝛽,𝛼)(−1), (1.8) that is, the Jacobi polynomials are even or odd depending on the order of the polynomial. In this form the polynomials may be generated using the starting form 𝑃₀^(𝛼,𝛽)(𝑥)𝑃=1,_1(𝛼,𝛽)(𝑥)=12𝛼−𝛽+(𝜆+1)𝑥, (1.9) such that𝜆=𝛼+𝛽+1, (1.10) which is obtained from Rodrigue's formula as follows:𝑃_𝑛^(𝛼,𝛽)(𝑥)=12^𝑛𝑛_𝑖=0𝑖()𝑛+𝛼𝑛+𝛽𝑛−𝑘1−𝑥^{(𝑛−𝑘)}(1+𝑥)^𝑘. (1.11)

Following the two seminal papers of Doha [3, 4] let 𝑓(𝑥) be an infinitely differentiable function defined on [−1, 1]; then we can write 𝑓(𝑥)=^∞_𝑛=0𝑎_𝑛𝑃_𝑛^𝛼,𝛽(𝑥), (1.12) and, for the 𝑞th derivative of 𝑓(𝑥), 𝑓^(𝑞)(𝑥)=^∞_𝑛=0𝑎_𝑛^(𝑞)𝑃_𝑛^(𝛼,𝛽)(𝑥),𝑎_𝑛⁽⁰⁾=𝑎_𝑛. (1.13)

Then, 𝑎_𝑛^(𝑞)=2^∞−𝑞_𝑖=0𝑛+𝑖+𝑞+𝜆−1_𝑞𝐶_{𝑛+𝑖,𝑛}𝛼+𝑞,𝛽+𝑞,𝛼,𝛽𝑎_{𝑛+𝑖+𝑞},∀𝑛≥0,𝑞≥1, (1.14) where 𝐶_{𝑛+𝑖,𝑛}𝛼+𝑞,𝛽+𝑞,𝛼,𝛽=𝑛+𝑖+2𝑞+𝜆−1_𝑛𝑛+𝛼+𝑞+1_𝑖Γ(𝑛+𝜆)()×𝑖!Γ2𝑛+𝜆₃𝐹₂,(𝑚)−𝑖2𝑛+𝑖+2𝑞+𝜆𝑛+𝛼+1;1𝑛+𝑞+𝛼+12𝑛+𝜆+1_𝑛()⋯()=()!=𝑚𝑚+1𝑚+𝑛−1𝑚+𝑛−1()!.𝑚−1 (1.15) For the proof of the above, see [3]. The formula for the expansion coefficients of a general-order derivative of an infinitely differentiable function in terms of those of the function is available for expansions in ultraspherical and Jacobi polynomials in Doha [5]. Another interesting formula is 𝑥^𝑚𝑃_𝑗^(𝛼,𝛽)(𝑥)=^2𝑚_𝑛=0𝑎_𝑚𝑛𝑗𝑃^(𝛼,𝛽)_{𝑗+𝑚−𝑛}(𝑥),∀𝑚,𝑗≥0, (1.16) with 𝑃^(𝛼,𝛽)_−𝑟(𝑥)=0,𝑟≥1, (1.17) where 𝑎_𝑚𝑛𝑗=(−1)^𝑛2^{𝑗+𝑚−𝑛}𝑚!2𝑗+2𝑚−2𝑛+𝜆Γ𝑗+𝑚−𝑛+𝜆Γ𝑗+𝛼+1Γ𝑗+𝛽+1ΓΓΓ×𝑗+𝑚−𝑛+𝛼+1𝑗+𝑚−𝑛+𝛽+1𝑗+𝜆_{min𝑗+𝑚−𝑛,𝑗max0,𝑗−𝑛}𝑘𝑗+𝑚−𝑛Γ𝑗+𝑘+𝜆2^𝑘×𝑛+𝑘+𝑗!Γ3𝑗+2𝑚−2𝑛−𝑘+𝜆+1_{𝑗−𝑘𝑖=0}(−1)^𝑙Γ2𝑗+𝑚−𝑛−𝑘−𝑙+𝛼+1Γ𝑗+𝑚+𝑙−𝑛+𝛽+1Γ×𝑙!𝑗−𝑘−𝑙!Γ𝑗−𝑙+𝛼+1𝑘+𝑙+𝛽+1₂𝐹₁𝑗−𝑘−𝑛,𝑗+𝑚+𝑙−𝑛+𝛽+1;3𝑗+2𝑚−2𝑛−𝑘+𝜆+1;2. (1.18) For the proof see [4]. Chebyshev, Legendre, and ultraspherical polynomials are particular cases of the Jacobi polynomials. These polynomials have been used both in the solution of boundary value problems [6] and in computational fluid dynamics [7, 8]. For the ultraspherical coefficients of the moments of a general-order derivative of an infinitely differentiable function, see [5]. Collocation method is a kind of spectral method that uses the delta function as a test function. Test functions have an important role because these functions are applied to obtain the minimum value for residual by using inner product. Attention to this point is very important since Tau and collocation methods do not give the best approximation solution of ODEs. Due to adding the boundary conditions, we are forced to eliminate several conditions of orthogonality properties, and this decreases the accuracy of the ODEs solution. Previous explanations do not mean that the accuracy of the collocation method is less than that of the Galerkin method. Many studies on the collocation method have recently appeared in the literature on the numerical solution of boundary value problems [9–12]. The collocation solution is a piecewise continuous polynomial function. The error has been analyzed for different conditions by different authors. Frank and Ueberhuber [13] showed equivalence between solutions of collocation method and fixed points of Iterated Defect Correction (IDeC) method. They proved that IDeC method can be regarded as an efficient scheme for solving collocation equations. Çelik [14, 15] investigated the corrected collocation method for the approximate computation of eigenvalues of Sturm-Liouville and periodic Sturm-Liouville problems by a Truncated Chebyshev Series (TCS). On the other hand, some results on Jacobi approximation were used for analyzing the 𝑝-version of finite element method, boundary element method, spectral method for axisymmetrical domain, and various rational spectral methods [16–21].

The rest of this paper is organized as follows. In Section 2, the method of solution is developed. In Section 3, we report our numerical findings and demonstrate the accuracy of the proposed scheme by considering a numerical example. Finally, a brief conclusion is presented in Section 4.

2. Jacobi Collocation Method

We consider the class of ODEs with the following form^𝑠_𝑖=0𝑃_𝑖(𝑥)𝑦^′′𝑖+^𝑙_𝑖=0𝑄_𝑖(𝑥)𝑦′^𝑖+^𝑚_𝑖=0𝑅_𝑖(𝑥)𝑦^𝑖=𝐹(𝑥), (2.1) where 𝑠,𝑙,𝑚≥0, 𝑃_𝑖(𝑥), 𝑄_𝑖(𝑥), 𝑅_𝑖(𝑥), and 𝐹(𝑥) are defined on the interval 𝑎≤𝑥≤𝑏. The most general boundary conditions are𝛾₀𝑦(𝑎)+𝛾₁𝑦^′(𝑎)𝛼=𝛿,₀𝑦(𝑏)+𝛼₁𝑦^′(𝑏)=𝜇, (2.2) where the real coefficients 𝛾_𝑖,𝛼_𝑖, 𝛿, and 𝜇 are appropriate constants and the solution is expressed in the form𝑦=^𝑛_𝑟=0𝑐_𝑟𝑃_𝑟^(𝛼,𝛽)(𝑥), (2.3) under the certain conditions (𝑟=0,1,2,…).

Any finite range 𝑎≤𝑤≤𝑏 can be transformed to the basic range −1≤𝑥≤1 with the change of variable()()𝑤=𝑏−𝑎𝑥+𝑏+𝑎2. (2.4)

Hence, there is no loss of generality in taking 𝑎=−1 and 𝑏=1. To obtain such a solution, the Jacobi collocation method can be used. Collocation points can be taken as the roots of 𝑃_𝑛^(𝛼,𝛽)(𝑥) [22]. Until now we have not the explicit formulas that determine these roots. Therefore, we use the numerical method to determine these roots, for example, the Newton method. Substituting the Jacobi collocation points in (2.1), the following expression can be obtained:^𝑠_𝑖=0𝑃_𝑖𝑥_𝑗^𝑛_𝑟=0𝑐_𝑟𝑃_𝑟^(𝛼,𝛽)𝑥_𝑗^′′𝑖+^𝑙_𝑖=0𝑄_𝑖𝑥_𝑗^𝑛_𝑟=0𝑐_𝑟𝑃_𝑟^(𝛼,𝛽)𝑥_𝑗^′𝑖+^𝑚_𝑖=0𝑅_𝑖𝑥_𝑗^𝑛_𝑟=0𝑐_𝑟𝑃_𝑟^(𝛼,𝛽)𝑥_𝑗^𝑖=𝐹𝑥_𝑗,𝑗=0,1,2,3,…,𝑛. (2.5)

The above system of equations constructs 𝑛+1 nonlinear equations. We impose the boundary conditions to these equations as follows:𝛾₀^𝑛_𝑟=0𝑐_𝑟𝑃_𝑟^(𝛼,𝛽)(𝑎)+𝛾₁^𝑛_𝑟=0𝑐_𝑟𝑃_𝑟^(𝛼,𝛽)(𝑎)^′𝛼=𝛿,_0𝑛𝑟=0𝑐_𝑟𝑃_𝑟^(𝛼,𝛽)(𝑏)+𝛼₁^𝑛_𝑟=0𝑐_𝑟𝑃_𝑟^(𝛼,𝛽)(𝑏)^′=𝜇. (2.6) Equations (2.5) and (2.6) construct the 𝑛+3 nonlinear equations and so we eliminate two equations from (2.1). Thus we have 𝑛+1 equations and solve this system with Newton iteration method.

3. Illustrative Numerical Example

Consider the following nonlinear boundary value problem:cos(𝑥)𝑦^′′2−sin²(𝑥)𝑦′=0, (3.1) subject to boundary conditions𝑦(0)𝑦𝜋=0,2=1. (3.2)

The exact solution is 𝑦=sin(𝑥). We solve it with the above-mentioned method (for different values of 𝛼 and 𝛽), and we show the results. We compute the numerical solution for 4 terms in Figure 1.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1 

Numerical solution of our example for different values of 𝛼 and 𝛽. (𝑎1)𝛼=0, 𝛽=0;(𝑎2) is error. (𝑏1)𝛼=−0.4, 𝛽=−0.4;(𝑏2) is error. (𝑐1)𝛼=−0.5, 𝛽=−0.5;(𝑐2) is error. (𝑑1)𝛼=−0.4, 𝛽=0.6;(𝑑2) is error. (𝑒1)𝛼=0.6, 𝛽=0.6;(𝑒2) is error. (𝑓1)𝛼=2, 𝛽=2;(𝑓2) is error.

At first we obtain the roots of 𝑃_𝑛^(𝛼,𝛽)(𝑥). For this goal we must determine the 𝛼 and 𝛽 parameters, and Table 1 shows the roots.

Table 1 

The roots of 𝑃_𝑛^(𝛼,𝛽)(𝑥) for different values of 𝛼 and 𝛽.

To apply the method we assume that the solution is in the following form:𝑦=³_𝑖=0𝑐_𝑖𝑃_𝑖^(𝛼,𝛽)(𝑥). (3.3) By substituting 𝑦 in (3.1), we obtain a system of equations, and, with attention to Gaussian nodes, we must solve the set of nonlinear equations. The boundary conditions are imposed on the system of equations with elimination of two equations, and we add these boundary conditions to the system as follows:cos𝑥_𝑗³_𝑖=0𝑐_𝑖𝑃_𝑖^(𝛼,𝛽)𝑥_𝑗^′′2−sin²𝑥_𝑗³_𝑖=0𝑐_𝑖𝑃_𝑖^(𝛼,𝛽)𝑥_𝑗^′=0,𝑗=0,1,2,3. (3.4) The above system can be solved with Newton iteration method. We show our results for every parameter as in Table 2. Representation of the error by increasing the number of terms in our polynomial solution is shown in Table 3. It is observable that as the number of terms increases, the absolute error decreases.

Table 2 

The solution of our example for different values of 𝛼 and 𝛽.

Table 3 

Representation of the error by increasing the number of terms in our polynomial solution.

4. Conclusion

In this paper, we show that this method is very accurate, even for small value of 𝑛 (i.e., 𝑛=4). We see that Chebyshev polynomial (the cases 𝑐₁ and 𝑐₂) and Legendre polynomial (the cases 𝑎₁ and 𝑎₂) have not obtained the best approximation (see the case 𝛼=−0.4,𝛽=0.6). We observe that the error is symmetric when 𝛼=𝛽 (see the cases 𝑎₂,𝑏₂,𝑐₂,𝑒₂, and 𝑓₂). When we increase the number of terms in our computation, we monitor that the error decreases and our computed results are the best.