# Collocation Method via Jacobi Polynomials for Solving Nonlinear Ordinary Differential Equations

Authors: Ahmad Imani; Azim Aminataei; Ali Imani
Journal: International Journal of Mathematics and Mathematical Sciences
 (2011)

---

## 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.


---

## Body

## 1. Introduction

The Jacobi polynomials with respect to parametersα>-1, β>-1 (see, e.g., [1, 2]) are sequences of polynomials Pn(α,β)(x);(n=0,1,2,…) satisfying the following relation
(1.1)∫-11(1-x)α(1+x)βPn(α,β)(x)Pm(α,β)(x)dx={0,m≠n,hn,m=n,
where
(1.2)hn=2α+β+1Γ(α+n+1)Γ(β+n+1)(2n+α+β+1)n!Γ(n+α+β+1).These polynomials are eigenfunctions of the following singular Sturm-Liouville equation:(1.3)(1-x2)Φ′′(x)+[β-α-(β+α+2)x]Φ′(x)+n(n+β+α+1)Φ(x)=0.A consequence of this is that spectral accuracy can be achieved for expansions in Jacobi polynomials so that(1.4)Pn(α,β)(1)=Γ(n+α+1)n!Γ(α+1),Pn(α,β)(-1)=(-1)nΓ(n+β+1)n!Γ(β+1).The Jacobi polynomials can be obtained from Rodrigue's formula as(1.5)(1-x)α(1+x)βPn(α,β)(x)=(-1)ndn2nn!dxn[(1-x)α+n(1+x)β+n].
Furthermore, we have that (1.6)ddx(Pn(α,β)(x))=12(n+α+β+1)Pn-1(α+1,β+1)(x).The Jacobi polynomials are normalized such that(1.7)Pn(α,β)(1)=(n+αn)=Γ(n+α+1)n!Γ(α+1).An important consequence of the symmetry of weight functionw(x) and the orthogonality of Jacobi polynomial is the symmetric relation (1.8)Pn(α,β)(1)=(-1)nPn(β,α)(-1),
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.9)P0(α,β)(x)=1,P1(α,β)(x)=12(α-β+(λ+1)x),
such that(1.10)λ=α+β+1,
which is obtained from Rodrigue's formula as follows:(1.11)Pn(α,β)(x)=12n∑i=0n(n+αi)(n+βn-k)(1-x)(n-k)(1+x)k.Following the two seminal papers of Doha [3, 4] let f(x) be an infinitely differentiable function defined on [−1, 1]; then we can write (1.12)f(x)=∑n=0∞anPnα,β(x),
and, for the qth derivative of f(x), (1.13)f(q)(x)=∑n=0∞an(q)Pn(α,β)(x),an(0)=an.Then,(1.14)an(q)=2-q∑i=0∞(n+i+q+λ-1)qCn+i,n(α+q,β+q,α,β)an+i+q,∀n≥0,q≥1,
where (1.15)Cn+i,n(α+q,β+q,α,β)=(n+i+2q+λ-1)n(n+α+q+1)iΓ(n+λ)i!Γ(2n+λ)×F23(-i2n+i+2q+λn+α+1;1n+q+α+12n+λ+1),(m)n=m(m+1)⋯(m+n-1)=(m+n-1)!(m-1)!.
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 (1.16)xmPj(α,β)(x)=∑n=02mamn(j)Pj+m-n(α,β)(x),∀m,j≥0,
with (1.17)P-r(α,β)(x)=0,r≥1,
where
(1.18)amn(j)=(-1)n2j+m-nm!(2j+2m-2n+λ)Γ(j+m-n+λ)Γ(j+α+1)Γ(j+β+1)Γ(j+m-n+α+1)Γ(j+m-n+β+1)Γ(j+λ)×∑max⁡(0,j-n)min⁡(j+m-n,j)(j+m-nk)Γ(j+k+λ)2k(n+k+j)!Γ(3j+2m-2n-k+λ+1)×∑i=0j-k(-1)lΓ(2j+m-n-k-l+α+1)Γ(j+m+l-n+β+1)l!(j-k-l)!Γ(j-l+α+1)Γ(k+l+β+1)×F12(j-k-n,j+m+l-n+β+1;3j+2m-2n-k+λ+1;2).
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 p-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 Section2, 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(2.1)∑i=0sPi(x)(y′′)i+∑i=0lQi(x)(y′)i+∑i=0mRi(x)(y)i=F(x),
where s,l,m≥0, Pi(x), Qi(x), Ri(x), and F(x) are defined on the interval a≤x≤b. The most general boundary conditions are(2.2)γ0y(a)+γ1y′(a)=δ,α0y(b)+α1y′(b)=μ,
where the real coefficients γi,αi, δ, and μ are appropriate constants and the solution is expressed in the form(2.3)y=∑r=0ncrPr(α,β)(x),
under the certain conditions (r=0,1,2,…).Any finite rangea≤w≤b can be transformed to the basic range -1≤x≤1 with the change of variable(2.4)w=(b-a)x+(b+a)2.Hence, there is no loss of generality in takinga=-1 and b=1. To obtain such a solution, the Jacobi collocation method can be used. Collocation points can be taken as the roots of Pn(α,β)(x) [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:(2.5)∑i=0sPi(xj)(∑r=0ncr(Pr(α,β)(xj))′′)i+∑i=0lQi(xj)(∑r=0ncr(Pr(α,β)(xj))′)i+∑i=0mRi(xj)(∑r=0ncrPr(α,β)(xj))i=F(xj),j=0,1,2,3,…,n.The above system of equations constructsn+1 nonlinear equations. We impose the boundary conditions to these equations as follows:(2.6)γ0(∑r=0ncrPr(α,β)(a))+γ1(∑r=0ncrPr(α,β)(a))′=δ,α0(∑r=0ncrPr(α,β)(b))+α1(∑r=0ncrPr(α,β)(b))′=μ.
Equations (2.5) and (2.6) construct the n+3 nonlinear equations and so we eliminate two equations from (2.1). Thus we have n+1 equations and solve this system with Newton iteration method.


## 3. Illustrative Numerical Example

Consider the following nonlinear boundary value problem:(3.1)cos⁡(x)(y′′)2-sin2(x)y′=0,
subject to boundary conditions(3.2)y(0)=0,y(π2)=1.The exact solution isy=sin(x). 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.Numerical solution of our example for different values ofα and β. (a1)α=0, β=0;(a2) is error. (b1)α=-0.4, β=-0.4;(b2) is error. (c1)α=-0.5, β=-0.5;(c2) is error. (d1)α=-0.4, β=0.6;(d2) is error. (e1)α=0.6, β=0.6;(e2) is error. (f1)α=2, β=2;(f2) is error.
(a)(b)(c)(d)(e)(f)At first we obtain the roots ofPn(α,β)(x). For this goal we must determine the α and β parameters, and Table 1 shows the roots.Table 1
The roots ofPn(α,β)(x) for different values of α and β.
α=0,α=-0.4,α=-0.5,α=-0.4,α=0.6,α=2,β=0β=-0.4β=-0.5β=0.6α=0.6β=2−0.8611−0.9103−0.69470.9284−0.7996−0.69470.86110.91030.6947−0.75620.79960.6947−0.3400−0.3729−0.25050.4879−0.3038−0.25060.34000.37290.2505−0.48790.30380.2506To apply the method we assume that the solution is in the following form:(3.3)y=∑i=03ciPi(α,β)(x).
By substituting y 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:(3.4)cos⁡(xj)(∑i=03ci(Pi(α,β)(xj))′′)2-sin2(xj)(∑i=03ciPi(α,β)(xj))′=0,j=0,1,2,3.
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 β.
α=0,
β=0y=-0.00000000000000001119663700028745952445299264438P0(0,0)(x)+0.91522505937265707954150015923236P1(0,0)(x)-0.000000000000000022393274000574919048905985288767P2(0,0)(x)-0.059677652961529090927417051547721P3(0,0)(x).α=-0.4,
β=-0.4y=-0.000000000000000009996746174510901192865699495P0(-0.4,-0.4)(x)+0.89426538703542076709570198303656P1(-0.4,-0.4)(x)-0.000000000000000011996095223297443607007612548232P2(-0.4,-0.4)(x)-0.04198983897503790669842860235471P3(-0.4,-0.4)(x).α=-0.5,
β=-0.5y=+0.0000000000000000311183604968765645740547197244P0(-0.5,-0.5)(x)+0.88795403952304732584548430952028P1(-0.5,-0.5)(x)+0.000000000000000031118360496876564574054719724411P2(-0.5,-0.5)(x)-0.036586425139067737724980961913629P3(-0.5,-0.5)(x).α=-0.4,
β=0.6y=+0.40613910584685335757685249174382P0(-0.4,0.6)(x)+0.48556046784167101822183413097838P1(-0.4,0.6)(x)-0.018735918702094560665272214650799P2(-0.4,0.6)(x)-0.021911869735090742971371812015998P3(-0.4,0.6)(x).α=0.6,
β=0.6y=+0.00000000000000000185396073409571348071696317294P0(0.6,0.6)(x)+0.93795482243803032015459053856673P1(0.6,0.6)(x)+0.0000000000000000059326745972604669042025846148469P2(0.6,0.6)(x)-0.078409581597779066879404197800853P3(0.6,0.6)(x).α=2,
β=2y=+0.000000000000000002334262949979475861874127912143P0(2,2)(x)+0.97023755744312467875938569010755P1(2,2)(x)+0.00000000000000001400557769987685517124476747286P2(2,2)(x)-0.10421967857002249036265594452531P3(2,2)(x).Table 3
Representation of the error by increasing the number of terms in our polynomial solution.
(a)N=10N=12N=14(α,β)=(0,0)1.02×10-92.13×10-102.65×10-12(α,β)=(-0.4,-0.4)2.26×10-103.23×10-112.12×10-13(α,β)=(-0.5,-0.5)2.33×10-101.18×10-101.05×10-13(α,β)=(-0.4,0.6)7.16×10-101.01×10-103.55×10-12(α,β)=(0.6,0.6)6.37×10-82.24×10-93.70×10-11(α,β)=(2,2)4.25×10-61.66×10-72.39×10-9(b)N=16N=18N=20(α,β)=(0,0)1.31×10-151.21×10-182.09×10-22(α,β)=(-0.4,-0.4)4.23×10-167.45×10-199.76×10-22(α,β)=(-0.5,-0.5)1.62×10-153.28×10-201.12×10-23(α,β)=(-0.4,0.6)5.60×10-144.08×10-195.35×10-22(α,β)=(0.6,0.6)2.36×10-131.85×10-173.19×10-20(α,β)=(2,2)4.17×10-113.06×10-137.11×10-16


## 4. Conclusion

In this paper, we show that this method is very accurate, even for small value ofn (i.e., n=4). We see that Chebyshev polynomial (the cases c1 and c2) and Legendre polynomial (the cases a1 and a2) have not obtained the best approximation (see the case α=-0.4,β=0.6). We observe that the error is symmetric when α=β (see the cases a2,b2,c2,e2, and f2). When we increase the number of terms in our computation, we monitor that the error decreases and our computed results are the best.



---

*Source: 673085-2011-05-26.xml*