# A Jacobi Dual-Petrov-Galerkin Method for Solving Some Odd-Order Ordinary Differential Equations

Authors: E. H. Doha; A. H. Bhrawy; R. M. Hafez
Journal: Abstract and Applied Analysis
 (2011)

---

## Abstract

A Jacobi dual-Petrov-Galerkin (JDPG) method is introduced and used for solving fully integrated reformulations of third- and fifth-order ordinary differential equations (ODEs) with constant coefficients. The reformulated equation for theJth order ODE involves n-fold indefinite integrals for
n=1,…,J. Extension of the JDPG for ODEs with polynomial coefficients is treated using the Jacobi-Gauss-Lobatto quadrature. Numerical results with comparisons are given to confirm the reliability of the proposed method for some constant and polynomial coefficients ODEs.


---

## Body

## 1. Introduction

A well-known advantage of spectral methods is high accuracy with relatively fewer unknowns when compared with low-order finite-difference methods [1, 2]. On the other hand, spectral methods typically give rise to full matrices, partially negating the gain in efficiency due to the fewer degrees of freedom. In general, the use of the Jacobi polynomials (Pn(α,β) with α,β∈(-1,∞) and n is the polynomial degree) has the advantage of obtaining solutions of ordinary differential equations (ODEs) in terms of the Jacobi indices (see for instance, [3–5]). Several such pairs (α,β) have been used for approximate solutions of ODEs (see [6–10]). We avoid developing approximation results for each particular pair of indices and instead carry out a study with general indices. With this motivation, we introduce in this paper a family of the Jacobi polynomials with general indices.Third-order differential equations have applications in many engineering models, see for instance [11–14]). Fifth-order differential equations generally arise in the mathematical modeling of viscoelastic flows and other branches of mathematical physics and engineering sciences, see [15–17]. Existence and uniqueness of solutions of such boundary value problems are discussed, for instance, in [18].In this paper, the proposed differential equations are integratedq times, where q is the order of the equation, and we make use of the formulae relating the expansion coefficients of integration appearing in these integrated forms of the proposed differential equations to the Jacobi polynomials themselves (see, Doha [19]). An advantage of this approach is that the general equation in the algebraic system then contains a finite number of terms. We, therefore, motivated our interest in integrated forms of these differential equations. The interested reader is referred to Doha and Bhrawy [7, 20].The main aim of this paper is to propose a suitable way to approximate some integrated forms of third- and fifth-order ODEs with constant coefficients using a spectral method, based on the Jacobi polynomials such that it can be implemented efficiently and at the same time has a good convergence property. It is worthy to note here that odd-order problems lack the symmetry of even-order ones, so we propose a Jacobi dual-Petrov-Galerkin (JDPG) method. The method leads to systems with specially structured matrices that can be efficiently inverted. We apply the method for solving the integrated forms of third- and fifth-order ODEs by using compact combinations of the Jacobi polynomials, which satisfy essentially all the underlying homogeneous boundary conditions. To be more precise, for the JDPG we choose the trial functions to satisfy the underlying boundary conditions of the differential equations, and we choose the test functions to satisfy the dual boundary conditions. Extension of the JDPG for polynomial coefficient ODEs is obtained by approximating the weighted inner products in the JDPG by using the Jacobi-Gauss-Lobatto quadrature. Finally, examples are given to illustrate the efficiency and implementation of the method. Comparisons are made to confirm the reliability of the method.The remainder of this paper is organized as follows. In Section2 we give an overview of the Jacobi polynomials and their relevant properties needed hereafter. Section 3 is devoted to the theoretical derivation of the JDPG method for third-order differential equations with constant and variable polynomial coefficients. Section 4 gives the corresponding results for those obtained in Section 3, but for the fifth-order differential equations with constant coefficient and two choices of boundary conditions. In Section 5, we present some numerical results exhibiting the accuracy and efficiency of our numerical algorithms. Some concluding remarks are given in the final section.


## 2. Preliminaries

LetSN(I) be the space of polynomials of degree at most N on the interval I=(-1,1). We set(2.1)WN={u∈SN:u(±1)=u(1)(1)=0},WN*={u∈SN:u(±1)=u(1)(-1)=0},
and let Pn(α,β)(x)(n=0,1,2,…) be the Jacobi polynomials orthogonal with the weight functions wα,β(x)=(1-x)α(1+x)β, where α,β>-1.LetxN,j(α,β), 0≤j≤N, be the zeros of (1-x2)∂xPN(α,β). Denote by ϖN,j(α,β), 0≤j≤N, the weights of the corresponding Gauss-Lobatto quadrature formula, which are arranged in decreasing order. We define the discrete inner product and norm of weighted space Lwα,β2(I) as follows:(2.2)(u,v)wα,β,N=∑k=0Nu(xN,k(α,β))v(xN,k(α,β))ϖN,k(α,β),∥u∥wα,β,N=(u,u)wα,β,N.
Obviously (see, e.g., formula (2.25) of [21]),(2.3)(u,v)wα,β,N=(u,v)wα,β∀u,v∈S2N-1.
Thus, for any u∈SN, the norms ∥u∥wα,β,N and ∥u∥wα,β coincide.For anyu∈C(I), the Jacobi-Gauss-Lobatto interpolation operator INP(α,β)u(x)∈SN, satisfying(2.4)INP(α,β)u(xN,k(α,β))=u(xN,k(α,β)),0≤k≤N.
We also denote by INc=INP(-1/2,-1/2) and INl=INP(0,0) the Chebyshev-Gauss-Lobatto and Legendre-Gauss-Lobatto interpolation operators, respectively.For any real numbersα,β>-1, the set of the Jacobi polynomials forms a complete Lwα,β2(I)-orthogonal system, and(2.5)(Pk(α,β),Pj(α,β))w(α,β)=hkδk,j,
where δk,j is the Kronecker function and(2.6)hk=2α+β+1Γ(k+α+1)Γ(k+β+1)(2k+α+β+1)Γ(k+1)Γ(k+α+β+1).The following special values will be of fundamental importance in what follows (see, [22–24])(2.7)Pn(α,β)(1)=(α+1)nn!,Pn(α,β)(-1)=(-1)n(β+1)nn!,DqPn(α,β)(1)=∏i=0q-1Γ(n+α+1)(n+λ+i)2q(n-q)!Γ(q+α+1),DqPn(α,β)(-1)=(-1)n+qDqPn(β,α)(1),
where (a)k=Γ(a+k)/Γ(a) and λ=1+β+α.If we define theq times repeated differentiation and integration of Pn(α,β)(x) by DqPn(α,β)(x) and In(q,α,β)(x), respectively, then (cf. Doha [19, 22])(2.8)DqPn(α,β)(x)=2-q(n+λ)q∑i=0n-qCn-q,i(q,α,β)Pi(α,β)(x),(2.9)In(q,α,β)(x)=2q(n-q+λ)q∑i=qn+qCn+q,i(-q,α,β)Pi(α,β)(x)+πq-1(x),q≥0,n≥q+1forα=β=-12,q≥0,n≥qforα≠-12orβ≠-12,
where(2.10)Cl,i(q,α,β)=(l+2q+λ)i(i+q+α+1)l-iΓ(i+λ)(l-i)!Γ(2i+λ)×F32(-l+il+2q+i+λi+α+1i+q+α+12i+λ+11),
with πq-1(x) being a polynomial of degree at most (q-1). It is to be noted that In(q,α,β)(x) may be obtained from DqPn(α,β)(x) by replacing q with negative q. In general, the hypergeometric series F32 cannot be summed in explicit form, but it can be summed by Watsons identity [25], if α=β. The following two lemmas will be of fundamental importance in what follows.Lemma 2.1 (see, [19, 26]).
One has(2.11)xl=∑j=0lDlj(γ,δ)Pj(γ,δ)(x),
where
(2.12)Dlj(γ,δ)=2jl!(1+γ+δ+2j)∑i=0l-j(-2)i(1+γ+j)ii!(l-i-j)!(1+γ+δ+j)i+1(2+γ+δ+i+j)j.Lemma 2.2.
If one writes(2.13)Ik(q,α,β)(x)=∑i=qk+qSq(k,i,α,β)Pi(α,β)(x)+πq-1(x),
then
(2.14)Sq(k,i,α,β)=2q(k-q+λ)qCk+q,i(-q,α,β),q=1,2,….Proof.
It is immediately obtained from relation (2.9).


## 3. Third-Order Differential Equation

We are interested in using the JDPG method to solve the third-order differential equation(3.1)u(3)(x)+γ1u(2)(x)+γ2u(1)(x)+γ3u(x)=g(x),inI,
subject to(3.2)u(±1)=u(1)(1)=0,
where γ1, γ2, and γ3 are constants and g is a given source function. In this paper, we consider the fully integrated form of the ODE, given by(3.3)u(x)+γ1∫(1)u(x)(dx)(1)+γ2∫(2)u(x)(dx)(2)+γ3∫(3)u(x)(dx)(3)=f(x)+∑i=02diPi(α,β)(x),u(±1)=u′(1)=0,
where(3.4)∫(q)u(x)(dx)(q)=∫⋯∫︷qtimesu(x)dxdx⋯dx︷qtimes,f(x)=∫(3)g(x)(dx)(3).
In this work we assume that f, the three-fold indefinite integral form of g, can be evaluated analytically. We set(3.5)SN=span⁡{P0(α,β)(x),P1(α,β)(x),…,PN(α,β)(x)},
then the dual-Petrov-Galerkin approximation to (3.3) is to find uN∈WN such that(3.6)(uN,v)wα,β+γ1(∫(1)uN(dx)(1),v)wα,β+γ2(∫(2)uN(dx)(2),v)wα,β+γ3(∫(3)uN(dx)(3),v)wα,β=(f+∑i=02diPi(α,β),v)wα,β∀v∈WN*,
where wα,β(x)=(1-x)α(1+x)β and (u,v)w=∫Iuvwdx is the inner product in the weighted space Lwα,β2(I). The norm in Lwα,β2(I) will be denoted by ∥·∥wα,β.


### 3.1. The Jacobi Dual-Petrov-Galerkin Method

We choose compact combinations of the Jacobi polynomials as basis functions to minimize the bandwidth hoping to improve the condition number of the coefficient matrix corresponding to (3.6). We choose the test basis and trial functions of expansions ϕk(x) and ψk(x) to be of the form(3.7)ϕk(x)=Pk(α,β)(x)+ϵkPk+1(α,β)(x)+ɛkPk+2(α,β)(x)+ζkPk+3(α,β)(x),(3.8)ψk(x)=Pk(α,β)(x)+ρkPk+1(α,β)(x)+ϱkPk+2(α,β)(x)+σkPk+3(α,β)(x),
where ϵk, ɛk, ζk, ρk, ϱk, and σk are the unique constants such that ϕk(x)∈WN and ψk(x)∈WN*. From the boundary conditions, ϕk(±1)=ϕk(1)(1)=0 and (2.7), hence ϵk, ɛk and ζk can be uniquely determined by using mathematica to give (see, [27])(3.9)ϵk=-(k+1)(2k+λ+2)(k-α+2β+1)(k+α+1)(k+β+1)(2k+λ+4),ɛk=-(k+1)2(2k+λ+1)(k-β+2α+3)(k+α+1)2(k+β+1)(2k+λ+5),ζk=(k+1)3(2k+λ+1)2(k+α+1)2(k+β+1)(2k+λ+4)2.
Using (2.7), and ψk(±1)=ψk(1)(-1)=0, one verifies readily that(3.10)ρk=(k+1)(2k+λ+2)(k-β+2α+1)(k+α+1)(k+β+1)(2k+λ+4),ϱk=-(k+1)2(2k+λ+1)(k-α+2β+3)(k+β+1)2(k+α+1)(2k+λ+5),σk=-(k+1)3(2k+λ+1)2(k+β+1)2(k+α+1)(2k+λ+4)2.Now it is clear that (3.6) is equivalent to(3.11)(uN,ψk(x))wα,β+γ1(∫(1)uN(dx)(1),ψk(x))wα,β+γ2(∫(2)uN(dx)(2),ψk(x))wα,β+γ3(∫(3)uN(dx)(3),ψk(x))wα,β=(f(x)+∑i=02diPi(α,β)(x),ψk(x))wα,β,N,k=0,1,…,N,
where (·,·)wα,β,N is the discrete inner product associated with the Jacobi-Gauss-Lobatto quadrature. The constants d0, d1, and d2 would not appear if we take k≥3 in (3.11), therefore we get(3.12)(uN,ψk(x))wα,β+γ1(∫(1)uN(dx)(1),ψk(x))wα,β+γ2(∫(2)uN(dx)(2),ψk(x))wα,β+γ3(∫(3)uN(dx)(3),ψk(x))wα,β=(f(x),ψk(x))wα,β,N,k=3,4,…,N.
If we take ϕk(x) and ψk(x) as defined in (3.7) and (3.8), respectively, and if we denote(3.13)fk=(f,ψk(x))wα,β,N,f=(f3,f4,…,fN)T,uN(x)=∑n=0N-3vnϕn(x),v=(v0,v1,…,vN-3)T,akj=(ϕj-3(x),ψk(x))wα,β,bkj=(∫(1)ϕj-3(x)(dx)(1),ψk(x))wα,β,ckj=(∫(2)ϕj-3(x)(dx)(2),ψk(x))wα,β,dkj=(∫(3)ϕj-3(x)(dx)(3),ψk(x))wα,β,
then(3.14)WN=span⁡{ϕ0(x),ϕ1(x),…,ϕN-3(x)},WN*=span⁡{ψ3(x),ψ4(x),…,ψN(x)},
and the nonzero elements (akj), (bkj), (ckj), and (dkj) for 3≤k,j≤N are given  as  follows:(3.15)akk=ζk-3hk,ak,k+1=ɛk-2hk+ζk-2ρkhk+1,ak,k+2=ϵk-1hk+ɛk-1ρkhk+1+ζk-1ϱkhk+2,ak,k+3=hk+ϵkρkhk+1+ɛkϱkhk+2+ζkσkhk+3,ak,k+4=ρkhk+1+ϵk+1ϱkhk+2+ɛk+1σkhk+3,ak,k+5=ϱkhk+2+ϵk+2σkhk+3,ak,k+6=σkhk+3,bk,j=[R1(j,k,α,β)hk+R1(j,k+1,α,β)ρkhk+1+R1(j,k+2,α,β)ϱkhk+2+R1(j,k+3,α,β)σkhk+3],j=k+l-1,l=0,1,…,8,ck,j=[R2(j,k,α,β)hk+R2(j,k+1,α,β)ρkhk+1+R2(j,k+2,α,β)ϱkhk+2+R2(j,k+3,α,β)σkhk+3],j=k+l-2,l=0,1,…,10,dk,j=R3(j,k,α,β)hk+R3(j,k+1,α,β)ρkhk+1+R3(j,k+2,α,β)ϱkhk+2+R3(j,k+3,α,β)σkhk+3,j=k+l-3,l=0,1,…,12,
where(3.16)Ri(j,k,α,β)=Si(j-3,k,α,β)+ϵj-3Si(j-2,k,α,β)+ɛj-3Si(j-1,k,α,β)+ζj-3Si(j,k,α,β).Hence by setting(3.17)A=(akj),B=(bkj),C=(ckj),D=(dkj),3≤k,j≤N,
then (3.12) is equivalent to the following matrix equation:(3.18)(A+γ1B+γ2C+γ3D)v=f.
All the analytical formulae of the nonzero elements of matrices A, B, C  and  D can be obtained by direct computations using the properties of the Jacobi polynomials (for details see, [27, 28]).In the case ofα,β≠0, α,β∈(-1,∞), we can form explicitly the LU factorization, that is, A+γ1B+γ2C+γ3D=LU. In general, the expense of calculating LU factorization of an N×N dense matrix M is O(N3) operations, and the expense of solving Ax=b, provided that the factorization is known, is O(N2) (see, [27]). However, in the case of banded matrix A of bandwidth r, we need just O(r2N) operations to factorize and O(rN) operations to solve a linear system. In the case of γi≠0, i=0,1,2, the square matrix A+γ1B+γ2C+γ3D has bandwidth of 13. So we need just O(13N) operations to solve the linear system (3.18). If r≪N, this represents a very substantial saving.It is worthy to note that, forα,β≠0, the algebraic system (3.18) resulting from fully integrated reformulation of (3.1) is sparse and is therefore cheaper to solve than those obtained from the differentiated form (see [27, Theorem 3.1]). Moreover, the savings in computational effort increase as the size of the systems grow. Thus, we have demonstrated the advantage of using the integrated forms over the differentiated ones for constant coefficients ODEs.


### 3.2. A Quadrature JDPG Method

The JDPG can be extended for ODEs with polynomial coefficients because of analytical form of a product of an algebraic polynomial, and the Jacobi polynomials are known.Now the formula of the Jacobi coefficients of the moments of one single Jacobi polynomial of any degree (see, Doha [22]) is(3.19)xmPj(α,β)(x)=∑n=02mΘm,n(j)Pj+m-n(α,β)(x)∀m,j≥0,
with P-r(α,β)(x)=0, r≥1, where(3.20)Θm,n(j)=(-1)n2j+m-nm!(2j+2m-2n+λ)Γ(j+m-n+λ)Γ(j+α+1)Γ(j+β+1)Γ(j+m-n+α+1)Γ(j+m-n+β+1)Γ(j+λ)×∑k=max⁡⁡(0,j-n)min⁡⁡(j+m-n,j)(j+m-nk)Γ(j+k+λ)2k(n+k-j)!Γ(3j+2m-2n-k+λ+1)×∑l=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)×F21(j-k-n,j+m+l-n+β+1;3j+2m-2n-k+λ+1;2).
This formula can be used to facilitate greatly the setting up of the algebraic systems to be obtained by applying the spectral methods for solving differential equations with polynomial coefficients of any order.Let us consider the following integrated form of the third-order differential equation:(3.21)u(x)+∫(1)γ1(x)u(x)(dx)(1)+∫(2)γ2(x)u(x)(dx)(2)+∫(3)γ3(x)u(x)(dx)(3)=f(x)+∑i=02diPi(α,β)(x),inI=(-1,1),u(±1)=u(1)(1)=0,
where γ1(x), γ2(x), and γ3(x) are the variable polynomial coefficients of the differential equation. The quadrature dual-Petrov-Galerkin method for (3.21) is to find uN∈WN such that(3.22)(uN,vN)wα,β,N+(∫(1)γ1(x)uN(dx)(1),vN)wα,β,N+(∫(2)γ2(x)uN(dx)(2),vN)wα,β,N+(∫(3)γ3(x)uN(dx)(3),vN)wα,β,N=(f+∑i=02diPi(α,β),vN)wα,β,N∀vN∈WN*,
where (u,v)wα,β,N is the discrete inner product of u and v associated with the Jacobi-Gauss-Lobatto quadrature.Let us consider(3.23)uN=∑k=0N-3ãkϕk,a¯=(ã0,ã1,…,ãN-3)T,f̃k=(f,ψk)wα,β,N,f¯=(f̃3,f̃1,…,f̃N)T,
and using Lemma 2.2 and formula (3.19), we can obtain(3.24)ãij=(ϕj-3,ψi)wα,β,N,b̃ij=(∫(1)γ1(x)ϕj-3(x)(dx)(1),ψi)wα,β,N,c̃ij=(∫(2)γ2(x)ϕj-3(x)(dx)(2),ψi)wα,β,N,d̃ij=(∫(3)γ3(x)ϕj-3(x)(dx)(3),ψi)wα,β,N.
And by setting(3.25)Ã=(ãkj),B̃=(b̃kj),C̃=(c̃kj),D̃=(d̃kj),3≤k,j≤N,
then the linear system of (3.22) becomes(3.26)(Ã+B̃+C̃+D̃)a¯=f¯.


## 3.1. The Jacobi Dual-Petrov-Galerkin Method

We choose compact combinations of the Jacobi polynomials as basis functions to minimize the bandwidth hoping to improve the condition number of the coefficient matrix corresponding to (3.6). We choose the test basis and trial functions of expansions ϕk(x) and ψk(x) to be of the form(3.7)ϕk(x)=Pk(α,β)(x)+ϵkPk+1(α,β)(x)+ɛkPk+2(α,β)(x)+ζkPk+3(α,β)(x),(3.8)ψk(x)=Pk(α,β)(x)+ρkPk+1(α,β)(x)+ϱkPk+2(α,β)(x)+σkPk+3(α,β)(x),
where ϵk, ɛk, ζk, ρk, ϱk, and σk are the unique constants such that ϕk(x)∈WN and ψk(x)∈WN*. From the boundary conditions, ϕk(±1)=ϕk(1)(1)=0 and (2.7), hence ϵk, ɛk and ζk can be uniquely determined by using mathematica to give (see, [27])(3.9)ϵk=-(k+1)(2k+λ+2)(k-α+2β+1)(k+α+1)(k+β+1)(2k+λ+4),ɛk=-(k+1)2(2k+λ+1)(k-β+2α+3)(k+α+1)2(k+β+1)(2k+λ+5),ζk=(k+1)3(2k+λ+1)2(k+α+1)2(k+β+1)(2k+λ+4)2.
Using (2.7), and ψk(±1)=ψk(1)(-1)=0, one verifies readily that(3.10)ρk=(k+1)(2k+λ+2)(k-β+2α+1)(k+α+1)(k+β+1)(2k+λ+4),ϱk=-(k+1)2(2k+λ+1)(k-α+2β+3)(k+β+1)2(k+α+1)(2k+λ+5),σk=-(k+1)3(2k+λ+1)2(k+β+1)2(k+α+1)(2k+λ+4)2.Now it is clear that (3.6) is equivalent to(3.11)(uN,ψk(x))wα,β+γ1(∫(1)uN(dx)(1),ψk(x))wα,β+γ2(∫(2)uN(dx)(2),ψk(x))wα,β+γ3(∫(3)uN(dx)(3),ψk(x))wα,β=(f(x)+∑i=02diPi(α,β)(x),ψk(x))wα,β,N,k=0,1,…,N,
where (·,·)wα,β,N is the discrete inner product associated with the Jacobi-Gauss-Lobatto quadrature. The constants d0, d1, and d2 would not appear if we take k≥3 in (3.11), therefore we get(3.12)(uN,ψk(x))wα,β+γ1(∫(1)uN(dx)(1),ψk(x))wα,β+γ2(∫(2)uN(dx)(2),ψk(x))wα,β+γ3(∫(3)uN(dx)(3),ψk(x))wα,β=(f(x),ψk(x))wα,β,N,k=3,4,…,N.
If we take ϕk(x) and ψk(x) as defined in (3.7) and (3.8), respectively, and if we denote(3.13)fk=(f,ψk(x))wα,β,N,f=(f3,f4,…,fN)T,uN(x)=∑n=0N-3vnϕn(x),v=(v0,v1,…,vN-3)T,akj=(ϕj-3(x),ψk(x))wα,β,bkj=(∫(1)ϕj-3(x)(dx)(1),ψk(x))wα,β,ckj=(∫(2)ϕj-3(x)(dx)(2),ψk(x))wα,β,dkj=(∫(3)ϕj-3(x)(dx)(3),ψk(x))wα,β,
then(3.14)WN=span⁡{ϕ0(x),ϕ1(x),…,ϕN-3(x)},WN*=span⁡{ψ3(x),ψ4(x),…,ψN(x)},
and the nonzero elements (akj), (bkj), (ckj), and (dkj) for 3≤k,j≤N are given  as  follows:(3.15)akk=ζk-3hk,ak,k+1=ɛk-2hk+ζk-2ρkhk+1,ak,k+2=ϵk-1hk+ɛk-1ρkhk+1+ζk-1ϱkhk+2,ak,k+3=hk+ϵkρkhk+1+ɛkϱkhk+2+ζkσkhk+3,ak,k+4=ρkhk+1+ϵk+1ϱkhk+2+ɛk+1σkhk+3,ak,k+5=ϱkhk+2+ϵk+2σkhk+3,ak,k+6=σkhk+3,bk,j=[R1(j,k,α,β)hk+R1(j,k+1,α,β)ρkhk+1+R1(j,k+2,α,β)ϱkhk+2+R1(j,k+3,α,β)σkhk+3],j=k+l-1,l=0,1,…,8,ck,j=[R2(j,k,α,β)hk+R2(j,k+1,α,β)ρkhk+1+R2(j,k+2,α,β)ϱkhk+2+R2(j,k+3,α,β)σkhk+3],j=k+l-2,l=0,1,…,10,dk,j=R3(j,k,α,β)hk+R3(j,k+1,α,β)ρkhk+1+R3(j,k+2,α,β)ϱkhk+2+R3(j,k+3,α,β)σkhk+3,j=k+l-3,l=0,1,…,12,
where(3.16)Ri(j,k,α,β)=Si(j-3,k,α,β)+ϵj-3Si(j-2,k,α,β)+ɛj-3Si(j-1,k,α,β)+ζj-3Si(j,k,α,β).Hence by setting(3.17)A=(akj),B=(bkj),C=(ckj),D=(dkj),3≤k,j≤N,
then (3.12) is equivalent to the following matrix equation:(3.18)(A+γ1B+γ2C+γ3D)v=f.
All the analytical formulae of the nonzero elements of matrices A, B, C  and  D can be obtained by direct computations using the properties of the Jacobi polynomials (for details see, [27, 28]).In the case ofα,β≠0, α,β∈(-1,∞), we can form explicitly the LU factorization, that is, A+γ1B+γ2C+γ3D=LU. In general, the expense of calculating LU factorization of an N×N dense matrix M is O(N3) operations, and the expense of solving Ax=b, provided that the factorization is known, is O(N2) (see, [27]). However, in the case of banded matrix A of bandwidth r, we need just O(r2N) operations to factorize and O(rN) operations to solve a linear system. In the case of γi≠0, i=0,1,2, the square matrix A+γ1B+γ2C+γ3D has bandwidth of 13. So we need just O(13N) operations to solve the linear system (3.18). If r≪N, this represents a very substantial saving.It is worthy to note that, forα,β≠0, the algebraic system (3.18) resulting from fully integrated reformulation of (3.1) is sparse and is therefore cheaper to solve than those obtained from the differentiated form (see [27, Theorem 3.1]). Moreover, the savings in computational effort increase as the size of the systems grow. Thus, we have demonstrated the advantage of using the integrated forms over the differentiated ones for constant coefficients ODEs.


## 3.2. A Quadrature JDPG Method

The JDPG can be extended for ODEs with polynomial coefficients because of analytical form of a product of an algebraic polynomial, and the Jacobi polynomials are known.Now the formula of the Jacobi coefficients of the moments of one single Jacobi polynomial of any degree (see, Doha [22]) is(3.19)xmPj(α,β)(x)=∑n=02mΘm,n(j)Pj+m-n(α,β)(x)∀m,j≥0,
with P-r(α,β)(x)=0, r≥1, where(3.20)Θm,n(j)=(-1)n2j+m-nm!(2j+2m-2n+λ)Γ(j+m-n+λ)Γ(j+α+1)Γ(j+β+1)Γ(j+m-n+α+1)Γ(j+m-n+β+1)Γ(j+λ)×∑k=max⁡⁡(0,j-n)min⁡⁡(j+m-n,j)(j+m-nk)Γ(j+k+λ)2k(n+k-j)!Γ(3j+2m-2n-k+λ+1)×∑l=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)×F21(j-k-n,j+m+l-n+β+1;3j+2m-2n-k+λ+1;2).
This formula can be used to facilitate greatly the setting up of the algebraic systems to be obtained by applying the spectral methods for solving differential equations with polynomial coefficients of any order.Let us consider the following integrated form of the third-order differential equation:(3.21)u(x)+∫(1)γ1(x)u(x)(dx)(1)+∫(2)γ2(x)u(x)(dx)(2)+∫(3)γ3(x)u(x)(dx)(3)=f(x)+∑i=02diPi(α,β)(x),inI=(-1,1),u(±1)=u(1)(1)=0,
where γ1(x), γ2(x), and γ3(x) are the variable polynomial coefficients of the differential equation. The quadrature dual-Petrov-Galerkin method for (3.21) is to find uN∈WN such that(3.22)(uN,vN)wα,β,N+(∫(1)γ1(x)uN(dx)(1),vN)wα,β,N+(∫(2)γ2(x)uN(dx)(2),vN)wα,β,N+(∫(3)γ3(x)uN(dx)(3),vN)wα,β,N=(f+∑i=02diPi(α,β),vN)wα,β,N∀vN∈WN*,
where (u,v)wα,β,N is the discrete inner product of u and v associated with the Jacobi-Gauss-Lobatto quadrature.Let us consider(3.23)uN=∑k=0N-3ãkϕk,a¯=(ã0,ã1,…,ãN-3)T,f̃k=(f,ψk)wα,β,N,f¯=(f̃3,f̃1,…,f̃N)T,
and using Lemma 2.2 and formula (3.19), we can obtain(3.24)ãij=(ϕj-3,ψi)wα,β,N,b̃ij=(∫(1)γ1(x)ϕj-3(x)(dx)(1),ψi)wα,β,N,c̃ij=(∫(2)γ2(x)ϕj-3(x)(dx)(2),ψi)wα,β,N,d̃ij=(∫(3)γ3(x)ϕj-3(x)(dx)(3),ψi)wα,β,N.
And by setting(3.25)Ã=(ãkj),B̃=(b̃kj),C̃=(c̃kj),D̃=(d̃kj),3≤k,j≤N,
then the linear system of (3.22) becomes(3.26)(Ã+B̃+C̃+D̃)a¯=f¯.


## 4. Fifth-Order Differential Equation

In this section, we consider the fifth-order differential equation of the form(4.1)u(5)+γ1u(3)+γ2u(1)+γ3u=g(x),x∈I,
but by considering its integrated form, namely,(4.2)u(x)+γ1∫(2)u(x)(dx)(2)+γ2∫(4)u(x)(dx)(4)+γ3∫(5)u(x)(dx)(5)=f(x)+∑i=04diPi(α,β)(x).


### 4.1. First Choice of Boundary Conditions

Here, we apply the dual-Petrov-Galerkin approximation to (4.2) subject to the boundary conditions(4.3)u(±1)=u(1)(±1)=u(2)(1)=0.We set(4.4)VN={v∈SN:v(±1)=v(1)(±1)=v(2)(1)=0},VN*={v∈SN:v(±1)=v(1)(±1)=v(2)(-1)=0};
then the Jacobi dual-Petrov-Galerkin approximation to (4.2) is to find uN∈VN such that (4.5)(uN,v)wα,β+γ1(∫(2)uN(dx)(2),v)wα,β+γ2(∫(4)uN(dx)(4),v)wα,β+γ3(∫(5)uN(dx)(5),v)wα,β=(f+∑i=04diPi(α,β),v)wα,β,N∀v∈VN*.
We consider the following the Jacobi dual-Petrov-Galerkin procedure for (4.1): find uN∈VN such that(4.6)(uN,v)wα,β+γ1(∫(2)uN(dx)(2),v)wα,β+γ2(∫(4)uN(dx)(4),v)wα,β+γ3(∫(5)uN(dx)(5),v)wα,β=(f+∑i=04diPi(α,β),v)wα,β∀v∈VN*.
Now, we choose the basis and the dual basis functions Φk(x) and Ψk(x) to be of the form(4.7)Φk(x)=Pk(α,β)(x)+ϵ̂kPk+1(α,β)(x)+ɛ̂kPk+2(α,β)(x)+ζ̂kPk+3(α,β)(x)+μ̂kPk+4(α,β)(x)+υ̂kPk+5(α,β)(x),k=0,1,…,N-5,Ψk(x)=Pk(α,β)(x)+ρ̂kPk+1(α,β)(x)+ϱ̂kPk+2(α,β)(x)+σ̂kPk+3(α,β)(x)+ς̂kPk+4(α,β)(x)+τ̂kPk+5(α,β)(x),k=0,1,…,N-5.
It is not difficult to show that the basis functions Φk(x)∈Vk+5 and the dual basis functions Ψk(x)∈Vk+5*.We choose the coefficientsϵ̂k, ɛ̂k, ζ̂k, μ̂k, and υ̂k such that Φk(x) verifies the boundary conditions (4.3). Making use of (2.7), then the boundary conditions (4.3) lead to linear system for these coefficients. The computation of the exact solution of such linear system for the unknown coefficients is extremely tedious by hand, and we have resorted to the symbolic computation software mathematica version 6, hence ϵ̂k, ɛ̂k, ζ̂k, μ̂k, and υ̂k can be uniquely determined to give(4.8)ϵ̂k=-(k+1)(2k+λ+2)(k-2α+3β+1)(k+α+1)(k+β+1)(2k+λ+6),ɛ̂k=-(k+1)2(2k+λ+1)(2k+λ+4)(k+α+1)2(k+β+1)2(2k+λ+6)2×[2k2+4(α+3)k+(-α2+9α+3β2+3(2α-1)β+16)],ζ̂k=(k+1)3(2k+λ+1)2(k+α+1)3(k+β+1)2(2k+λ+7)2×[2k2+4(β+3)k+(-3α2-3α-β2+3(2α+5)β+16)],μ̂k=(k+1)4(k+3α-2β+5)(2k+λ+1)3(k+α+1)3(k+β+1)2(2k+λ+6)2(2k+λ+9),υ̂k=-(k+1)5(2k+λ+1)4(k+α+1)3(k+β+1)2(2k+λ+6)4.
Since the dual basis functions Ψk(x) satisfy the dual boundary conditions, and making use of (2.7) then the unknown coefficients ρ̂k, ϱ̂k, σ̂k, ς̂k, and τ̂k are determined by using Mathematica to give(4.9)ρ̂k=(k+1)(2k+λ+2)(k+3α-2β+1)(k+α+1)(k+β+1)(2k+λ+6),ϱ̂k=-(k+1)2(2k+λ+1)(2k+λ+4)(k+α+1)2(k+β+1)2(2k+λ+6)2×[2k2+4(β+3)k+(-3α2+3α-β2+3(2α+3)β+16)],σ̂k=-(k+1)3(2k+λ+1)2(k+α+1)2(k+β+1)3(2k+λ+7)2×[2k2+4(α+3)k+(-α2+15α-3β2+3(2α-1)β+16)],ς̂k=(k+1)4(k+3β-2α+5)(2k+λ+1)3(k+α+1)2(k+β+1)3(2k+λ+6)2(2k+λ+9),τ̂k=(k+1)5(2k+λ+1)4(k+α+1)2(k+β+1)3(2k+λ+6)4.
It is clear that (4.6) is equivalent to(4.10)(uN,Ψk(x))wα,β+γ1(∫(2)uN(dx)(2),Ψk(x))wα,β+γ2(∫(4)uN(dx)(4),Ψk(x))wα,β+γ3(∫(5)uN(dx)(5),Ψk(x))wα,β=(f(x)+∑i=04diPi(α,β)(x),Ψk(x))wα,β,k=0,1,…,N.
The constants d0, d1, d2, d3, and d4 would not appear if we take k≥5 in (4.10), therefore we get(4.11)(uN,Ψk(x))wα,β+γ1(∫(2)uN(dx)(2),Ψk(x))wα,β+γ2(∫(4)uN(dx)(4),Ψk(x))wα,β+γ3(∫(5)uN(dx)(5),Ψk(x))wα,β=(f(x),Ψk(x))wα,β,k=5,6,…,N.
If we take Φk(x) and Ψk(x) as defined in (4.7) and if we denote(4.12)fk=(f,Ψk(x))wα,β,f=(f5,f6,…,fN)T,uN(x)=∑n=0N-5vnΦn(x),v=(v0,v1,…,vN-5)T,pkj=(Φj-5(x),Ψk(x))wα,β,qkj=(∫(2)Φj-5(x)(dx)(2),Ψk(x))wα,β,skj=(∫(4)Φj-5(x)(dx)(4),Ψk(x))wα,β,tkj=(∫(5)Φj-5(x)(dx)(5),Ψk(x))wα,β,
then(4.13)VN=span⁡{Φ0(x),Φ1(x),…,ΦN-5(x)},VN*=span⁡{Ψ5(x),Ψ6(x),…,ΨN(x)},
and the nonzero elements (pkj), (qkj), (skj), and (tkj) for 5≤k,j≤N are given  as  follows:(4.14)pkk=υ̂k-5hk,pk,k+1=μ̂k-4hk+υ̂k-4ρ̂khk+1,pk,k+2=ζ̂k-3hk+μ̂k-3ρ̂khk+1+υ̂k-3ϱ̂khk+2,pk,k+3=ɛ̂k-2hk+ζ̂k-2ρ̂khk+1+μ̂k-2ϱ̂khk+2+υ̂k-2σ̂khk+3,pk,k+4=ϵ̂k-1hk+ɛ̂k-1ρ̂khk+1+ζ̂k-1ϱ̂khk+2+μ̂k-1σ̂khk+3+υ̂k-1ς̂khk+4,pk,k+5=hk+ϵ̂kρ̂khk+1+ɛ̂kϱ̂khk+2+ζ̂kσ̂khk+3+μ̂kς̂khk+4+υ̂kτ̂khk+5,pk,k+6=ρ̂khk+1+ϵ̂k+1ϱ̂khk+2+ɛ̂k+1σ̂khk+3+ζ̂k+1ς̂khk+4+μ̂k+1τ̂khk+5,pk,k+7=ϱ̂khk+2+ϵ̂k+2σ̂khk+3+ɛ̂k+2ς̂khk+4+ζ̂k+2τ̂khk+5,pk,k+8=σ̂khk+3+ϵ̂k+3ς̂khk+4+ɛ̂k+3τ̂khk+5,pk,k+9=ς̂khk+4+ϵ̂k+4τ̂khk+5,pk,k+10=τ̂khk+5,(4.15)qk,j=R2(j,k,α,β)hk+R2(j,k+1,α,β)ρ̂khk+1+R2(j,k+2,α,β)ϱ̂khk+2+R2(j,k+3,α,β)σ̂khk+3+R2(j,k+4,α,β)ς̂khk+4+R2(j,k+5,α,β)τ̂khk+5,j=k+l-2,l=0,1,…,14,sk,j=R4(j,k,α,β)hk+R4(j,k+1,α,β)ρ̂khk+1+R4(j,k+2,α,β)ϱ̂khk+2+R4(j,k+3,α,β)σ̂khk+3+R4(j,k+4,α,β)ς̂khk+4+R4(j,k+5,α,β)τ̂khk+5,j=k+l-4,l=0,1,…,18,tk,j=R5(j,k,α,β)hk+R5(j,k+1,α,β)ρ̂khk+1+R5(j,k+2,α,β)ϱ̂khk+2+R5(j,k+3,α,β)σ̂khk+3+R5(j,k+4,α,β)ς̂khk+4+R5(j,k+5,α,β)τ̂khk+5,j=k+l-5,l=0,1,…,20,
where(4.16)Ri(j,k,α,β)=Si(j-5,k,α,β)+ϵ̂j-5Si(j-4,k,α,β)+ɛ̂j-5Si(j-3,k,α,β)+ζ̂j-5Si(j-2,k,α,β)+μ̂j-5Si(j-1,k,α,β)+υ̂j-5Si(j,k,α,β).
Then (4.11) is equivalent to the following matrix equation:(4.17)(P+γ1Q+γ2S+γ3T)v=f.


### 4.2. Second Choice of Boundary Conditions

In this subsection, we consider the fifth-order differential equation (4.1) with the following boundary conditions:(4.18)u(±1)=u(1)(±1)=u(3)(-1)=0.
Equation (4.1) subject to the boundary conditions (4.18) has been considered in [29, 30]. Let us denote(4.19)ZN={v∈SN:v(±1)=v(1)(±1)=v(3)(-1)=0},ẐN={v∈SN:v(±1)=v(1)(±1)=v(3)(1)=0};
then the Jacobi dual-Petrov-Galerkin approximation of (4.1) subject to (4.18) consists of finding uN∈ZN such that(4.20)(uN,v)wα,β+γ1(∫(2)uN(dx)(2),v)wα,β+γ2(∫(4)uN(dx)(4),v)wα,β+γ3(∫(5)uN(dx)(5),v)wα,β=(f+∑i=04diPi(α,β),v)wα,β∀v∈ẐN.
We consider the following choice of basis functions:(4.21)φk(x)=Pk(α,β)(x)+ξ1,kPk+1(α,β)(x)+ξ2,kPk+2(α,β)(x)+ξ3,kPk+3(α,β)(x)+ξ4,kPk+4(α,β)(x)+ξ5,kPk+5(α,β)(x),
and dual basis functions:(4.22)φ̂k(x)=Pk(α,β)(x)+δ1,kPk+1(α,β)(x)+δ2,kPk+2(α,β)(x)+δ3,kPk+3(α,β)(x)+δ4,kPk+4(α,β)(x)+δ5,kPk+5(α,β)(x),
where ξ1,k, ξ2,k, ξ3,k, ξ4,k, ξ5,k, δ1,k, δ2,k, δ3,k, δ4,k, and δ5,k  are chosen to be the unique constants such that φk(x)∈ZN and φ̂k(x)∈ẐN, for all k=0,1,…,N-5.Equation (4.20) is equivalent to the following matrix equation:(4.23)(P̂+γ1Q̂+γ2Ŝ+γ3T̂)v=f,
where the elements of the matrices  P̂,  Q̂,  Ŝ, and T̂ can be obtained similarly as in the previous sections, but details are not given here.


## 4.1. First Choice of Boundary Conditions

Here, we apply the dual-Petrov-Galerkin approximation to (4.2) subject to the boundary conditions(4.3)u(±1)=u(1)(±1)=u(2)(1)=0.We set(4.4)VN={v∈SN:v(±1)=v(1)(±1)=v(2)(1)=0},VN*={v∈SN:v(±1)=v(1)(±1)=v(2)(-1)=0};
then the Jacobi dual-Petrov-Galerkin approximation to (4.2) is to find uN∈VN such that (4.5)(uN,v)wα,β+γ1(∫(2)uN(dx)(2),v)wα,β+γ2(∫(4)uN(dx)(4),v)wα,β+γ3(∫(5)uN(dx)(5),v)wα,β=(f+∑i=04diPi(α,β),v)wα,β,N∀v∈VN*.
We consider the following the Jacobi dual-Petrov-Galerkin procedure for (4.1): find uN∈VN such that(4.6)(uN,v)wα,β+γ1(∫(2)uN(dx)(2),v)wα,β+γ2(∫(4)uN(dx)(4),v)wα,β+γ3(∫(5)uN(dx)(5),v)wα,β=(f+∑i=04diPi(α,β),v)wα,β∀v∈VN*.
Now, we choose the basis and the dual basis functions Φk(x) and Ψk(x) to be of the form(4.7)Φk(x)=Pk(α,β)(x)+ϵ̂kPk+1(α,β)(x)+ɛ̂kPk+2(α,β)(x)+ζ̂kPk+3(α,β)(x)+μ̂kPk+4(α,β)(x)+υ̂kPk+5(α,β)(x),k=0,1,…,N-5,Ψk(x)=Pk(α,β)(x)+ρ̂kPk+1(α,β)(x)+ϱ̂kPk+2(α,β)(x)+σ̂kPk+3(α,β)(x)+ς̂kPk+4(α,β)(x)+τ̂kPk+5(α,β)(x),k=0,1,…,N-5.
It is not difficult to show that the basis functions Φk(x)∈Vk+5 and the dual basis functions Ψk(x)∈Vk+5*.We choose the coefficientsϵ̂k, ɛ̂k, ζ̂k, μ̂k, and υ̂k such that Φk(x) verifies the boundary conditions (4.3). Making use of (2.7), then the boundary conditions (4.3) lead to linear system for these coefficients. The computation of the exact solution of such linear system for the unknown coefficients is extremely tedious by hand, and we have resorted to the symbolic computation software mathematica version 6, hence ϵ̂k, ɛ̂k, ζ̂k, μ̂k, and υ̂k can be uniquely determined to give(4.8)ϵ̂k=-(k+1)(2k+λ+2)(k-2α+3β+1)(k+α+1)(k+β+1)(2k+λ+6),ɛ̂k=-(k+1)2(2k+λ+1)(2k+λ+4)(k+α+1)2(k+β+1)2(2k+λ+6)2×[2k2+4(α+3)k+(-α2+9α+3β2+3(2α-1)β+16)],ζ̂k=(k+1)3(2k+λ+1)2(k+α+1)3(k+β+1)2(2k+λ+7)2×[2k2+4(β+3)k+(-3α2-3α-β2+3(2α+5)β+16)],μ̂k=(k+1)4(k+3α-2β+5)(2k+λ+1)3(k+α+1)3(k+β+1)2(2k+λ+6)2(2k+λ+9),υ̂k=-(k+1)5(2k+λ+1)4(k+α+1)3(k+β+1)2(2k+λ+6)4.
Since the dual basis functions Ψk(x) satisfy the dual boundary conditions, and making use of (2.7) then the unknown coefficients ρ̂k, ϱ̂k, σ̂k, ς̂k, and τ̂k are determined by using Mathematica to give(4.9)ρ̂k=(k+1)(2k+λ+2)(k+3α-2β+1)(k+α+1)(k+β+1)(2k+λ+6),ϱ̂k=-(k+1)2(2k+λ+1)(2k+λ+4)(k+α+1)2(k+β+1)2(2k+λ+6)2×[2k2+4(β+3)k+(-3α2+3α-β2+3(2α+3)β+16)],σ̂k=-(k+1)3(2k+λ+1)2(k+α+1)2(k+β+1)3(2k+λ+7)2×[2k2+4(α+3)k+(-α2+15α-3β2+3(2α-1)β+16)],ς̂k=(k+1)4(k+3β-2α+5)(2k+λ+1)3(k+α+1)2(k+β+1)3(2k+λ+6)2(2k+λ+9),τ̂k=(k+1)5(2k+λ+1)4(k+α+1)2(k+β+1)3(2k+λ+6)4.
It is clear that (4.6) is equivalent to(4.10)(uN,Ψk(x))wα,β+γ1(∫(2)uN(dx)(2),Ψk(x))wα,β+γ2(∫(4)uN(dx)(4),Ψk(x))wα,β+γ3(∫(5)uN(dx)(5),Ψk(x))wα,β=(f(x)+∑i=04diPi(α,β)(x),Ψk(x))wα,β,k=0,1,…,N.
The constants d0, d1, d2, d3, and d4 would not appear if we take k≥5 in (4.10), therefore we get(4.11)(uN,Ψk(x))wα,β+γ1(∫(2)uN(dx)(2),Ψk(x))wα,β+γ2(∫(4)uN(dx)(4),Ψk(x))wα,β+γ3(∫(5)uN(dx)(5),Ψk(x))wα,β=(f(x),Ψk(x))wα,β,k=5,6,…,N.
If we take Φk(x) and Ψk(x) as defined in (4.7) and if we denote(4.12)fk=(f,Ψk(x))wα,β,f=(f5,f6,…,fN)T,uN(x)=∑n=0N-5vnΦn(x),v=(v0,v1,…,vN-5)T,pkj=(Φj-5(x),Ψk(x))wα,β,qkj=(∫(2)Φj-5(x)(dx)(2),Ψk(x))wα,β,skj=(∫(4)Φj-5(x)(dx)(4),Ψk(x))wα,β,tkj=(∫(5)Φj-5(x)(dx)(5),Ψk(x))wα,β,
then(4.13)VN=span⁡{Φ0(x),Φ1(x),…,ΦN-5(x)},VN*=span⁡{Ψ5(x),Ψ6(x),…,ΨN(x)},
and the nonzero elements (pkj), (qkj), (skj), and (tkj) for 5≤k,j≤N are given  as  follows:(4.14)pkk=υ̂k-5hk,pk,k+1=μ̂k-4hk+υ̂k-4ρ̂khk+1,pk,k+2=ζ̂k-3hk+μ̂k-3ρ̂khk+1+υ̂k-3ϱ̂khk+2,pk,k+3=ɛ̂k-2hk+ζ̂k-2ρ̂khk+1+μ̂k-2ϱ̂khk+2+υ̂k-2σ̂khk+3,pk,k+4=ϵ̂k-1hk+ɛ̂k-1ρ̂khk+1+ζ̂k-1ϱ̂khk+2+μ̂k-1σ̂khk+3+υ̂k-1ς̂khk+4,pk,k+5=hk+ϵ̂kρ̂khk+1+ɛ̂kϱ̂khk+2+ζ̂kσ̂khk+3+μ̂kς̂khk+4+υ̂kτ̂khk+5,pk,k+6=ρ̂khk+1+ϵ̂k+1ϱ̂khk+2+ɛ̂k+1σ̂khk+3+ζ̂k+1ς̂khk+4+μ̂k+1τ̂khk+5,pk,k+7=ϱ̂khk+2+ϵ̂k+2σ̂khk+3+ɛ̂k+2ς̂khk+4+ζ̂k+2τ̂khk+5,pk,k+8=σ̂khk+3+ϵ̂k+3ς̂khk+4+ɛ̂k+3τ̂khk+5,pk,k+9=ς̂khk+4+ϵ̂k+4τ̂khk+5,pk,k+10=τ̂khk+5,(4.15)qk,j=R2(j,k,α,β)hk+R2(j,k+1,α,β)ρ̂khk+1+R2(j,k+2,α,β)ϱ̂khk+2+R2(j,k+3,α,β)σ̂khk+3+R2(j,k+4,α,β)ς̂khk+4+R2(j,k+5,α,β)τ̂khk+5,j=k+l-2,l=0,1,…,14,sk,j=R4(j,k,α,β)hk+R4(j,k+1,α,β)ρ̂khk+1+R4(j,k+2,α,β)ϱ̂khk+2+R4(j,k+3,α,β)σ̂khk+3+R4(j,k+4,α,β)ς̂khk+4+R4(j,k+5,α,β)τ̂khk+5,j=k+l-4,l=0,1,…,18,tk,j=R5(j,k,α,β)hk+R5(j,k+1,α,β)ρ̂khk+1+R5(j,k+2,α,β)ϱ̂khk+2+R5(j,k+3,α,β)σ̂khk+3+R5(j,k+4,α,β)ς̂khk+4+R5(j,k+5,α,β)τ̂khk+5,j=k+l-5,l=0,1,…,20,
where(4.16)Ri(j,k,α,β)=Si(j-5,k,α,β)+ϵ̂j-5Si(j-4,k,α,β)+ɛ̂j-5Si(j-3,k,α,β)+ζ̂j-5Si(j-2,k,α,β)+μ̂j-5Si(j-1,k,α,β)+υ̂j-5Si(j,k,α,β).
Then (4.11) is equivalent to the following matrix equation:(4.17)(P+γ1Q+γ2S+γ3T)v=f.


## 4.2. Second Choice of Boundary Conditions

In this subsection, we consider the fifth-order differential equation (4.1) with the following boundary conditions:(4.18)u(±1)=u(1)(±1)=u(3)(-1)=0.
Equation (4.1) subject to the boundary conditions (4.18) has been considered in [29, 30]. Let us denote(4.19)ZN={v∈SN:v(±1)=v(1)(±1)=v(3)(-1)=0},ẐN={v∈SN:v(±1)=v(1)(±1)=v(3)(1)=0};
then the Jacobi dual-Petrov-Galerkin approximation of (4.1) subject to (4.18) consists of finding uN∈ZN such that(4.20)(uN,v)wα,β+γ1(∫(2)uN(dx)(2),v)wα,β+γ2(∫(4)uN(dx)(4),v)wα,β+γ3(∫(5)uN(dx)(5),v)wα,β=(f+∑i=04diPi(α,β),v)wα,β∀v∈ẐN.
We consider the following choice of basis functions:(4.21)φk(x)=Pk(α,β)(x)+ξ1,kPk+1(α,β)(x)+ξ2,kPk+2(α,β)(x)+ξ3,kPk+3(α,β)(x)+ξ4,kPk+4(α,β)(x)+ξ5,kPk+5(α,β)(x),
and dual basis functions:(4.22)φ̂k(x)=Pk(α,β)(x)+δ1,kPk+1(α,β)(x)+δ2,kPk+2(α,β)(x)+δ3,kPk+3(α,β)(x)+δ4,kPk+4(α,β)(x)+δ5,kPk+5(α,β)(x),
where ξ1,k, ξ2,k, ξ3,k, ξ4,k, ξ5,k, δ1,k, δ2,k, δ3,k, δ4,k, and δ5,k  are chosen to be the unique constants such that φk(x)∈ZN and φ̂k(x)∈ẐN, for all k=0,1,…,N-5.Equation (4.20) is equivalent to the following matrix equation:(4.23)(P̂+γ1Q̂+γ2Ŝ+γ3T̂)v=f,
where the elements of the matrices  P̂,  Q̂,  Ŝ, and T̂ can be obtained similarly as in the previous sections, but details are not given here.


## 5. Numerical Results

In this section some examples are considered aiming to illustrate how one can apply the proposed algorithms presented in the previous sections. Comparisons between JDPG method and other methods proposed in [29–32] are made.Example 5.1.
Consider the one-dimensional third-order equation(5.1)u(3)(x)+πu(2)(x)+2πu(1)(x)+3πu(x)=f(x),x∈I,u(±1)=±12,u(1)(1)=π42,
with an exact smooth solution
(5.2)u(x)=sin⁡(πx4).Table1 lists the maximum pointwise error of  u-uN, using the JDPG with various choices of α, β, and N.Table 1
Maximum pointwise error using JDPG method forN=8,16,24, for Example 5.1.
NαβJDPGαβJDPG88.278·10-23.645·10-21612121.404·10-11-12123.368·10-12241.249·10-162.220·10-1681.165·10-15.644·10-21612-122.439·10-11-12-126.311·10-12242.220·10-162.220·10-16Example 5.2.
Consider the one-dimensional fifth-order differential problem(5.3)u(5)-γ1u(3)-γ2u(1)-γ3u=f(x),x∈I,u(±1)=0,u(1)(±1)=∓π2,u(2)(1)=-2,
with the exact solution u(x)=x2cos((π/2)x).Table2 lists the maximum pointwise error, using the JDPG method with various choices of α, β, γ1, γ2, γ3, and N. Numerical results of this example show that the JDPG method converges exponentially.Table 2
Maximum pointwise error using JDPG method forN=12,24,36, for Example 5.2.
Nαβγ1γ2γ3JDPGγ1γ2γ3JDPG123.294·10-19.117·10-12412121112.650·10-13NN2N31.052·10-10362.116·10-164.770·10-16123.024·10-18.366·10-124001111.982·10-13NN2N37.947·10-11363.053·10-163.053·10-16122.748·10-17.598·10-124-12-121111.437·10-13NN2N35.828·10-11362.116·10-161.387·10-16Example 5.3.
Consider the following fifth-order boundary value problem (see [29–32]):
(5.4)u(5)(x)-u(x)=-(15+10x)ex,x∈[0,1],u(0)=0,u(1)=0,u(1)(0)=1,u(1)(1)=-e,u(3)(0)=-3.
The analytic solution of this problem is u(x)=x(1-x)ex. Approximate solutions uN(x) (N=10,13,20,26,40) are obtained by using our proposed method. Table 3 exhibits a comparison between the error obtained by using JDPG method and the sixth-degree B-spline [31], sextic spline [32], nonpolynomial sextic spline [29], and the computational method in [30]. The numerical results show that JDPG method is more accurate than the existing methods.Table 3
Maximum pointwise error ofu-uN  for  N=15,18,25,31, for Example 5.3.
N-5αβOur algorithmIn [31]In [32]In [29]In [30]12122.50·10-610002.28·10-60.15702.2593·10-46.29887·10-5-12-122.05·10-612121.35·10-1013001.17·10-101.3767·10-45.91739·10-5-12-121.00·10-1012122.49·10-1620002.22·10-160.07471.3300·10-52.14116·10-6-12-121.87·10-1612122.22·10-1626001.18·10-167.1273·10-63.40705·10-7-12-121.94·10-16Example 5.4.
Consider the following third-order ODE with polynomial coefficients:(5.5)u(3)(x)+x2u(2)(x)+(x3+4x)u(1)(x)+(3x2+x+2)u(x)=g(x),x∈[-1,1],
subject to
(5.6)u(±1)=u(1)(1)=0,
where g is selected such that exact solution is u(x)=(1-x2)(1-x)e2x.Equation (5.5) can be rearranged to take the form(5.7)u(3)(x)+(x2u(x))(2)+(x3u(x))(1)+xu(x)=g(x),
and accordingly its fully integrated form is(5.8)u(x)+∫(1)x2u(x)(dx)(1)+∫(2)x3u(x)(dx)(2)+∫(3)xu(x)(dx)(3)=f(x)+∑i=02diPi(α,β)(x).Using the quadrature dual-Petrov-Galerkin method described in Section3.2, we evaluate the maximum pointwise error of u-uN with various choices of α, β, and N in Table 4. Numerical results show that there is a very good agreement between the approximate solution obtained by the quadrature JDPG method and the exact solution and at the same time ascertain that the JDPG method converges exponentially.Table 4
Maximum pointwise error using quadrature JDPG method forN=12,16,20,24, for Example 5.4.
NαβQuadrature JDPGαβQuadrature JDPG8002.643·10-212123.280·10-2121.496·10-52.060·10-5162.513·10-92.920·10-9203.613·10-145.847·10-14


## 6. Concluding Remarks

In this paper, we described a JDPG method for fully integrated forms of third- and fifth-order ODEs with constant coefficients. Because of the constant coefficients, the matrix elements of the discrete operators are provided explicitly, and this in turn greatly simplifies the steps and the computational effort for obtaining solutions. However, the integrated form of the source function (involving severalfold indefinite integrals) should be known analytically, and the right hand side vector require quadrature approximations. This approach is also considered for ODEs with polynomial coefficients. Numerical results exhibit the high accuracy of the proposed numerical methods of solutions.



---

*Source: 947230-2011-04-10.xml*