# Efficient Solutions of Multidimensional Sixth-Order Boundary Value Problems Using Symmetric Generalized Jacobi-Galerkin Method

Authors: E. H. Doha; W. M. Abd-Elhameed
Journal: Abstract and Applied Analysis
 (2012)

---

## Abstract

This paper presents some efficient spectral algorithms for solving linear sixth-order
two-point boundary value problems in one dimension based on the application of the
Galerkin method. The proposed algorithms are extended to solve the two-dimensional
sixth-order differential equations. A family of symmetric generalized Jacobi polynomials
is introduced and used as basic functions. The algorithms lead to linear systems with
specially structured matrices that can be efficiently inverted. The various matrix systems
resulting from the proposed algorithms are carefully investigated, especially their
condition numbers and their complexities. These algorithms are extensions to some of
the algorithms proposed by Doha and Abd-Elhameed (2002) and Doha and Bhrawy (2008) for second- and
fourth-order elliptic equations, respectively. Three numerical results are presented
to demonstrate the efficiency and the applicability of the proposed algorithms.


---

## Body

## 1. Introduction

The classical Jacobi polynomialsPn(α,β)(x) play important roles in mathematical analysis and its applications (see, e.g., [1–4]). In particular, the Legendre, the Chebyshev, and the ultraspherical polynomials have played important roles in spectral methods for partial differential equations (see, e.g., [5, 6]). It is proven that the Jacobi polynomials are precisely the only polynomials arising as eigenfunctions of a singular Sturm-Liouville problem, (see [7, Section 9.2]). This class of polynomials comprises all the polynomial solution to singular Sturm-Liouville problems on [−1,1].Spectral methods have developed rapidly over the past four decades. Their fascinating merit is the high accuracy; they have been applied successfully to numerical simulations of many problems in science and engineering. The spectral methods that are mostly used are based on the Chebyshev and Legendre approximations.Sixth-order boundary value problems arise in astrophysics, that is, the narrow convecting layers bounded by stable layers which are believed to surroundA-type stars may be modelled by sixth-order boundary value problems (see [8]).Chandrasekhar [9] determined that when an infinite horizontal layer of fluid is heated from below and is under the action of rotation, instability sets in. When this instability is as ordinary convection, the ordinary differential equation is sixth-order. Agarwal [10] presented the theorems stating the conditions for the existence and uniqueness of solutions of sixth order boundary value problems, while no numerical methods are contained therein. In [11] Bhrawy discussed the solution of sixth-order boundary value problems using Legendre Galerkin method. Lamnii et al. [12] used Spline collocation method for solving linear sixth-order boundary-value problems. Boutayeb and Twizell [8] developed a family of numerical methods for the solution of special nonlinear sixth-order boundary value problems. Siddiqi and Twizell [13] presented the solution of sixth order boundary value problem using the sextic spline. El-Gamel et al. [14] used Sinc-Galerkin method for the solutions of sixth order boundary value problems.Guo et al. [15] extended the definition of the classical Jacobi polynomials with indexes α,β>-1 to allow α and/or β to be negative integers. They showed also that the generalized Jacobi polynomials, with indexes corresponding to the number of boundary conditions in a given partial differential equation, are the natural basis  functions for the spectral approximation of this equation. Moreover, it is shown that the use of generalized Jacobi polynomials not only simplified the numerical analysis for the spectral approximations of differential equations, but also led to very efficient numerical algorithms.From the numerical point of view, Doha and Abd-Elhameed [16, 17], Doha and Bhrawy [18], and Doha et al. [19, 20]  have constructed efficient spectral Galerkin algorithms using compact combinations of orthogonal polynomials for solving elliptic equations of second-, third-, fourth-, and fifth-order equations in various situations.In this paper we are concerned with the direct solution techniques for sixth-order two-point boundary value problems, using symmetric generalized Jacobi-Galerkin approximations. Our algorithms lead to discrete linear systems with specially structured matrices that can be efficiently inverted.We organize the materials of this paper as follows. In Section2, we give some properties of classical and generalized Jacobi polynomials. In Section 3, we discuss two algorithms for solving the sixth-order elliptic linear differential equations subject to homogeneous and nonhomogeneous boundary conditions using symmetric generalized Jacobi Galerkin method (SGJGM). In Section 4, we explain how the idea of Section 3 can be extended to handle the sixth-order two dimensional differential equations. Three Numerical examples are given in Section 5 to show the efficiency of our algorithms. Some Concluding remarks are given in Section 6.


## 2. Some Properties of Classical and Generalized Jacobi Polynomials

### 2.1. Classical Jacobi Polynomials

The classical Jacobi polynomials, associated with the real parameters(α>-1,β>-1) (see [4, 21]), are a sequence of polynomials, Pn(α,β)(x),x∈(-1,1)(n=0,1,2,…), each, respectively, of degree n. For our present purposes, it is more convenient to introduce the normalized orthogonal polynomials Rn(α,β)(x)=(Pn(α,β)(x))/(Pn(α,β)(1)). This means that Rn(α,β)(x)=(n!Γ(α+1))/(Γ(n+α+1))Pn(α,β)(x). In such case Rn(α-(1/2),α-(1/2))(x) is identical to the ultraspherical polynomials Cn(α)(x), and the polynomials Rn(α,β)(x) may be generated using the recurrence relation (2.1)2(n+λ)(n+α+1)(2n+λ-1)Rn+1(α,β)(x)=(2n+λ-1)3xRn(α,β)(x)+(α2-β2)(2n+λ)Rn(α,β)(x)-2n(n+β)(2n+λ+1)Rn-1(α,β)(x),n=1,2,…,
starting from R0(α,β)(x)=1 and R1(α,β)(x)=(1/2(α+1))[α-β+(λ+1)x], or obtained from Rodrigues' formula(2.2)Rn(α,β)(x)=(-12)nΓ(α+1)Γ(n+α+1)(1-x)-α(1+x)-βDn[(1-x)α+n(1+x)β+n],
where (2.3)λ=α+β+1,(a)k=Γ(a+k)Γ(a),D=ddx,
and satisfy the orthogonality relation
(2.4)∫-11(1-x)α(1+x)βRm(α,β)(x)Rn(α,β)(x)dx={0,m≠n,hnα,β,m=n,
where(2.5)hnα,β=2λn!Γ(n+β+1)[Γ(α+1)]2(2n+λ)Γ(n+λ)Γ(n+α+1).
These polynomials are eigenfunctions of the following singular Sturm-Liouville equation: (2.6)(1-x2)ϕ′′(x)+[β-α-(λ+1)x]ϕ′(x)+n(n+λ)ϕ(x)=0.
The following relations will be of important use later:  (2.7)(1-x2)Rk-1(α+1,α+1)(x)=2(α+1)2k+2α+1[Rk-1(α,α)(x)-Rk+1(α,α)(x)],(2.8)DRk(α,α)(x)=k(k+2α+1)2(α+1)Rk-1(α+1,α+1)(x),k=1,2,….
The following two theorems are needed hereafter.Theorem 2.1.
The qth derivative of the symmetric normalized Jacobi polynomialRn(α,α)(x) is given explicitly by
(2.9)DqRk(α,α)(x)=2qk!(q-1)!Γ(k+2α+1)×∑m=0(k+m-q)evenk-q(m+α+1/2)Γ(m+2α+1)((k-m+q-2)/2)!Γ((k+m+q+2α+1)/2)m!((k-q-m)/2)!Γ((k+m-q+2α+3)/2)×Rm(α,α)(x),k≥q.(For the proof of Theorem2.1, see [22].)Theorem 2.2.
If one defines theq times repeated integration of the symmetric normalized Jacobi polynomials Rk(α,α)(x) by
(2.10)Ik(q,α)(x)=∬…∫︷qtimesRk(α,α)(x)dxdx⋯dx︷qtimes,
then
(2.11)Ik(q,α)(x)=2-qk!Γ(k+2α+1)∑j=0q(-1)j(qj)Γ(k-j+α+1/2)Γ(k+q-2j+2α+1)(k+q-2j)!Γ(k+q-j+α+3/2)×(k+q-2j+α+12)Rk+q-2j(α,α)(x),q≥0,k≥q+1forα=0,q≥0,k≥qforα≠0.(For the proof of Theorem2.2, see [23].)Also, the following two lemmas are needed in the sequel.Lemma 2.3.
For allk≥0, one has
(2.12)-D2[(1-x2)Rk(1,1)(x)]=(k+1)(k+2)Rk(1,1)(x).(For the proof of Lemma2.3, see [16].)Lemma 2.4.
For allk≥0, one has
(2.13)D4[(1-x2)2Rk(2,2)(x)]=(k+1)4Rk(2,2)(x).Proof.
Settingα=β=2 in relation (2.7), we get
(2.14)(1-x2)Rk(2,2)(x)=42k+5[Rk(1,1)(x)-Rk+2(1,1)(x)].
Making use of this relation and with the aid of Lemma 2.3, we obtain
(2.15)D4[(1-x2)2Rk(2,2)(x)]=42k+5D2[(k+3)(k+4)Rk+2(1,1)(x)-(k+1)(k+2)Rk(1,1)(x)],
which in turn gives with the aid of relation (2.8)
(2.16)D4[(1-x2)2Rk(2,2)(x)]=12k+5D[(k+2)4Rk+1(2,2)(x)-(k)4Rk-1(2,2)(x)].
Finally, from relation (2.9) (for q=1 and α=2), we get
(2.17)D4[(1-x2)2Rk(2,2)(x)]=(k+1)4Rk(2,2)(x).
This completes the proof of Lemma 2.4.


### 2.2. Generalized Jacobi Polynomials

Following Guo et al. [15], we define a family of generalized Jacobi polynomials/functions with indexes α, β∈ℝ.Letwα,β(x)=(1-x)α(1+x)β. We denote by Lwα,β2(-1,1) the weighted L2  space with inner product: (2.18)(u,v)wα,β(x):=∫Iu(x)v(x)wα,β(x)dx,
and the associated norm ∥u∥wα,β=(u,u)wα,β1/2. We are interested in defining Jacobi polynomials with indexes α and/or β≤-1, referred hereafter as generalized Jacobi polynomials (GJPs), in such a way that they satisfy some selected properties that are essentially relevant to spectral approximations. In this work, we will restrict our attention to the cases when α and β are negative integers.Letℓ,m∈ℤ (the set of all integers), and define (2.19)Jk(l,m)(x)={(1-x)-l(1+x)-mRk-k0(-l,-m)(x),k0=-(l+m),l,m≤-1,(1-x)-lRk-k0(-l,m)(x),k0=-l,l≤-1,m>-1,(1+x)-mRk-k0(l,-m)(x),k0=-m,l>-1,m≤-1,Rk-k0(l,m)(x),k0=0,l,m>-1.
An important property of the GJPs is that for  ℓ,m∈ℤ+,(2.20)DiJk(-l,-m)(1)=0,i=0,1,…,l-1,DjJk(-l,-m)(-1)=0,j=0,1,…,m-1.
Using relation (2.7), and after performing some manipulation, Jk(-3,-3)(x) can be written in terms of Legendre polynomials as:(2.21)Jk(-3,-3)(x)=48(2k-9)(2k-7)(2k-5)+[Lk-6(x)-3(2k-7)2k-3Lk-4(x)+3(2k-9)2k-1Lk-2(x)-(2k-7)(2k-9)(2k-3)(2k-1)Lk(x)].


## 2.1. Classical Jacobi Polynomials

The classical Jacobi polynomials, associated with the real parameters(α>-1,β>-1) (see [4, 21]), are a sequence of polynomials, Pn(α,β)(x),x∈(-1,1)(n=0,1,2,…), each, respectively, of degree n. For our present purposes, it is more convenient to introduce the normalized orthogonal polynomials Rn(α,β)(x)=(Pn(α,β)(x))/(Pn(α,β)(1)). This means that Rn(α,β)(x)=(n!Γ(α+1))/(Γ(n+α+1))Pn(α,β)(x). In such case Rn(α-(1/2),α-(1/2))(x) is identical to the ultraspherical polynomials Cn(α)(x), and the polynomials Rn(α,β)(x) may be generated using the recurrence relation (2.1)2(n+λ)(n+α+1)(2n+λ-1)Rn+1(α,β)(x)=(2n+λ-1)3xRn(α,β)(x)+(α2-β2)(2n+λ)Rn(α,β)(x)-2n(n+β)(2n+λ+1)Rn-1(α,β)(x),n=1,2,…,
starting from R0(α,β)(x)=1 and R1(α,β)(x)=(1/2(α+1))[α-β+(λ+1)x], or obtained from Rodrigues' formula(2.2)Rn(α,β)(x)=(-12)nΓ(α+1)Γ(n+α+1)(1-x)-α(1+x)-βDn[(1-x)α+n(1+x)β+n],
where (2.3)λ=α+β+1,(a)k=Γ(a+k)Γ(a),D=ddx,
and satisfy the orthogonality relation
(2.4)∫-11(1-x)α(1+x)βRm(α,β)(x)Rn(α,β)(x)dx={0,m≠n,hnα,β,m=n,
where(2.5)hnα,β=2λn!Γ(n+β+1)[Γ(α+1)]2(2n+λ)Γ(n+λ)Γ(n+α+1).
These polynomials are eigenfunctions of the following singular Sturm-Liouville equation: (2.6)(1-x2)ϕ′′(x)+[β-α-(λ+1)x]ϕ′(x)+n(n+λ)ϕ(x)=0.
The following relations will be of important use later:  (2.7)(1-x2)Rk-1(α+1,α+1)(x)=2(α+1)2k+2α+1[Rk-1(α,α)(x)-Rk+1(α,α)(x)],(2.8)DRk(α,α)(x)=k(k+2α+1)2(α+1)Rk-1(α+1,α+1)(x),k=1,2,….
The following two theorems are needed hereafter.Theorem 2.1.
The qth derivative of the symmetric normalized Jacobi polynomialRn(α,α)(x) is given explicitly by
(2.9)DqRk(α,α)(x)=2qk!(q-1)!Γ(k+2α+1)×∑m=0(k+m-q)evenk-q(m+α+1/2)Γ(m+2α+1)((k-m+q-2)/2)!Γ((k+m+q+2α+1)/2)m!((k-q-m)/2)!Γ((k+m-q+2α+3)/2)×Rm(α,α)(x),k≥q.(For the proof of Theorem2.1, see [22].)Theorem 2.2.
If one defines theq times repeated integration of the symmetric normalized Jacobi polynomials Rk(α,α)(x) by
(2.10)Ik(q,α)(x)=∬…∫︷qtimesRk(α,α)(x)dxdx⋯dx︷qtimes,
then
(2.11)Ik(q,α)(x)=2-qk!Γ(k+2α+1)∑j=0q(-1)j(qj)Γ(k-j+α+1/2)Γ(k+q-2j+2α+1)(k+q-2j)!Γ(k+q-j+α+3/2)×(k+q-2j+α+12)Rk+q-2j(α,α)(x),q≥0,k≥q+1forα=0,q≥0,k≥qforα≠0.(For the proof of Theorem2.2, see [23].)Also, the following two lemmas are needed in the sequel.Lemma 2.3.
For allk≥0, one has
(2.12)-D2[(1-x2)Rk(1,1)(x)]=(k+1)(k+2)Rk(1,1)(x).(For the proof of Lemma2.3, see [16].)Lemma 2.4.
For allk≥0, one has
(2.13)D4[(1-x2)2Rk(2,2)(x)]=(k+1)4Rk(2,2)(x).Proof.
Settingα=β=2 in relation (2.7), we get
(2.14)(1-x2)Rk(2,2)(x)=42k+5[Rk(1,1)(x)-Rk+2(1,1)(x)].
Making use of this relation and with the aid of Lemma 2.3, we obtain
(2.15)D4[(1-x2)2Rk(2,2)(x)]=42k+5D2[(k+3)(k+4)Rk+2(1,1)(x)-(k+1)(k+2)Rk(1,1)(x)],
which in turn gives with the aid of relation (2.8)
(2.16)D4[(1-x2)2Rk(2,2)(x)]=12k+5D[(k+2)4Rk+1(2,2)(x)-(k)4Rk-1(2,2)(x)].
Finally, from relation (2.9) (for q=1 and α=2), we get
(2.17)D4[(1-x2)2Rk(2,2)(x)]=(k+1)4Rk(2,2)(x).
This completes the proof of Lemma 2.4.


## 2.2. Generalized Jacobi Polynomials

Following Guo et al. [15], we define a family of generalized Jacobi polynomials/functions with indexes α, β∈ℝ.Letwα,β(x)=(1-x)α(1+x)β. We denote by Lwα,β2(-1,1) the weighted L2  space with inner product: (2.18)(u,v)wα,β(x):=∫Iu(x)v(x)wα,β(x)dx,
and the associated norm ∥u∥wα,β=(u,u)wα,β1/2. We are interested in defining Jacobi polynomials with indexes α and/or β≤-1, referred hereafter as generalized Jacobi polynomials (GJPs), in such a way that they satisfy some selected properties that are essentially relevant to spectral approximations. In this work, we will restrict our attention to the cases when α and β are negative integers.Letℓ,m∈ℤ (the set of all integers), and define (2.19)Jk(l,m)(x)={(1-x)-l(1+x)-mRk-k0(-l,-m)(x),k0=-(l+m),l,m≤-1,(1-x)-lRk-k0(-l,m)(x),k0=-l,l≤-1,m>-1,(1+x)-mRk-k0(l,-m)(x),k0=-m,l>-1,m≤-1,Rk-k0(l,m)(x),k0=0,l,m>-1.
An important property of the GJPs is that for  ℓ,m∈ℤ+,(2.20)DiJk(-l,-m)(1)=0,i=0,1,…,l-1,DjJk(-l,-m)(-1)=0,j=0,1,…,m-1.
Using relation (2.7), and after performing some manipulation, Jk(-3,-3)(x) can be written in terms of Legendre polynomials as:(2.21)Jk(-3,-3)(x)=48(2k-9)(2k-7)(2k-5)+[Lk-6(x)-3(2k-7)2k-3Lk-4(x)+3(2k-9)2k-1Lk-2(x)-(2k-7)(2k-9)(2k-3)(2k-1)Lk(x)].


## 3. Spectral-Galerkin Algorithms for One-Dimensional Sixth-Order Equations

In this section, we are interested in using SGJGM to solve the sixth-order two-point boundary value problems in one dimension subject to homogeneous and nonhomogeneous boundary conditions.


### 3.1. Homogeneous Boundary Conditions

Let us consider the sixth-order differential equation(3.1)-u(6)(x)+∑q=16δqη6-qu(6-q)(x)=f(x),x∈(-1,1),
subject to the homogeneous boundary conditions(3.2)u(j)(±1)=0,j=0,1,2,
where u(j)(x) denotes the jth derivative of u(x) with respect to x and {η6-q,q=1,…,6} are positive constants, and(3.3)δq={(-1)1+q/2,qeven,1,qodd.
Let us denote Hwr(I)(r=0,1,2,…), as the weighted Sobolev spaces, whose inner products and norms are denoted by (·,·)r,w and∥·∥r,w, respectively. To account for homogeneous boundary conditions, we define (3.4)H0,w3(I)={v∈Hw3(I):v(j)(±1)=0,0≤j≤2},
where v(j)(x)=(djv)/(dxj). The superscript w will be omitted in case of w=1.LetPN be the space of all polynomials of degree less than or equal to N. Setting VN=PN∩H03(I), then (3.5)VN:=span⁡{J6(-3,-3)(x),J7(-3,-3)(x),…,JN(-3,-3)(x)}.
The symmetric generalized Jacobi-Galerkin procedure for solving (3.1)-(3.2) is to find uN∈VN such that(3.6)(-D6uN(x),v(x))+∑q=16δqη6-q(D6-quN(x),v(x))=(f(x),v(x)),∀v∈VN,
where (u,v)=∫-11uvdx is the scalar inner product in the space L2(-1,1).


### 3.2. The Choice of Basis Functions

We choose the basis functions of expansion to be(3.7)ϕk(x)=Jk+6(-3,-3)(x)=(1-x2)3Rk(3,3)(x),k=0,1,…,N-6,
which fulfills the boundary conditions (3.2).It is obvious that{ϕk(x)} are linearly independent. Therefore, we have(3.8)VN=span⁡{ϕk(x):k=0,1,2,…,N-6}.
Now, the following two lemmas are needed hereafter.Lemma 3.1.
For allk≥0, one has
(3.9)-D6[Jk+6(-3,-3)(x)]=(k+1)6Rk(3,3)(x).Proof.
Settingα=β=3 in relation to (2.7), we get
(3.10)(1-x2)Rk(3,3)(x)=62k+7[Rk(2,2)(x)-Rk+2(2,2)(x)].
Making use of this relation and with the aid of Lemma 2.4, we obtain
(3.11)-D6[Jk+6(-3,-3)(x)]=62k+7D2[(k+3)4Rk+2(2,2)(x)-(k+1)4Rk(2,2)(x)].
The last relation with the aid of the two relations (2.8) and (2.9) yields
(3.12)-D6[Jk+6(-3,-3)(x)]=(k+1)6Rk(3,3)(x).
This completes the proof of Lemma 3.1.Lemma 3.2.
For allk≥0, one has
(3.13)D6-q[Jk+6(-3,-3)(x)]=∑j=0qdj,k,qRk-2j+q(3,3),1≤q≤6,
where
(3.14)dj,k,q=(-1)j+1(qj)Γ(k-j+(7/2))(k+q-2j+6)!(k+q-2j+(7/2))2q(k+q-2j)!Γ(k+q-j+(9/2)).Proof.
Integrating formula (3.9) q times, q∈{1,2,…,6}, and with the aid of relation (2.11) (in case of α=3), we obtain the desired formula.Based on the results of the two Lemmas3.1 and 3.2, we are able to state and prove the following two theorems.Theorem 3.3.
One has, for arbitrary constantsak,
(3.15)-D6[∑k=0N-6akJk+6(-3,-3)(x)]=∑k=0N-6bkRk(3,3)(x),
where
(3.16)bk=(k+1)6ak.Theorem 3.4.
One has, for arbitrary constantsak, and  1≤q≤6,
(3.17)D6-q[∑k=0N-6akJk+6(-3,-3)(x)]=∑k=0N+q-6rk,qRk(3,3)(x),where(3.18)rk,q=∑i=0qdq-i,k+q-2i,qak+q-2i.Now, the application of Galerkin method to (3.1) gives  (3.19)(-D6uN(x)+∑q=16δqη6-qD6-quN(x),ϕk(x))=(f(x),ϕk(x)),where(3.20)uN(x)=∑k=0N-6akϕk(x),ϕk(x)=Jk+6(-3,-3)(x),k=0,1,…,N-6.
The variational formulation (3.19) is equivalent to(3.21)(-D6uN(x),Rk(3,3)(x))w3,3(x)+∑q=16δqη6-q(D6-quN(x),Rk(3,3)(x))w3,3(x)=(f(x),Rk(3,3)(x))w3,3(x),w3,3(x)=(1-x2)3.
Substitution of formulae (3.15) and (3.17) into (3.21) yields (3.22)(∑j=0N-6bjRj(3,3)(x)+∑q=16∑j=0N+q-6δqη6-qrj,qRj(3,3)(x),Rk(3,3)(x))w3,3(x)=(f(x),Rk(3,3)(x))w3,3(x),
where bk and rk,q are as given by (3.16) and (3.18), respectively.Now, if we apply the orthogonality relation ofRk(3,3)(x) on (3.22), we obtain the following linear system of equations:  (3.23)(bk+∑q=16δqη6-qrk,q)hk3,3=fk;k=0,1,…,N-6,
where (3.24)fk=(f(x),Rk(3,3)(x))w3,3(x),(3.25)hk3,3=4608(k+1)6(2k+7).
The linear system (3.23) may be put in the form(3.26)(bk+∑q=16δqη6-qrk,q)=fk*,k=0,1,…,N-6,
where(3.27)fk*=fkhk3,3.
This system of equations may be put in the matrix form(3.28)(B+∑q=16η6-qG6-q)a=f*,
where (3.29)a=(a0,a1,…,aN-6)T,f*=(f0*,f1*,…,fN-6*)T,
and the nonzero elements of the matrices B and G6-q,1≤q≤6 are given explicitly in the following theorem.Theorem 3.5.
IfuN(x)=∑k=0N-6akJk+6(-3,-3)(x) is the symmetric generalized Jacobi-Galerkin approximation to (3.1)-(3.2), then the expansion coefficients {ak:k=0,1,…,N-6} satisfy the matrix system (3.28), where the nonzero elements of the matrices  B=(bkj) and G6-q=(gkj6-q)=(δqrkj6-q), 0≤k,j≤N-6, 1≤q≤6 are given as:
(3.30)bkk=(k+1)6,gk,k+q-2i6-q=δq(-1)q+i+1(2k+7)(k+6)!q!Γ(-i+k+(7/2))2q+1i!k!(q-i)!Γ(k-i+q+(9/2)),0≤i≤q.It is worthy to note here that the case corresponding toη6-q=0, 1≤q≤6 leads to a linear system with diagonal matrix. The result for such case is summarized in the following important corollary.Corollary 3.6.
IfuN(x)=∑k=0N-6akJk+6(-3,-3)(x) and  η6-q=0, 1≤q≤6, is the symmetric generalized Jacobi-Galerkin approximation to problem (3.1)-(3.2), then the expansion coefficients {ak:k=0,1,…,N-6} are given explicitly by
(3.31)ak=2k+74608∫-11(1-x2)3f(x)Rk(3,3)(x),k=0,1,…,N-6.Remark 3.7.
Each of the matricesG6-q, q=1,2,…,6 in (3.28) is a band matrix whose total number of nonzero diagonals upper or lower the main diagonal does not exceed q. Thus the coefficient matrix D=B+∑q=16η6-qG6-q is seven-band matrix at most. This special structure of D simplifies greatly the solution of the linear system (3.28). The system in such case can be factorized by LU-decomposition and the number of operations necessary to construct this factorization is of order 78(N-5), and the number of operations needed to solve the two triangular systems is of order 25(N-5).Note
The total number of operations mentioned in the previous discussion includes the number of all subtractions, additions, divisions, and multiplications (see [24]).


### 3.3. Condition Number

For the direct collocation method, the condition number behaves likeO(N12) (N: maximal degree of polynomials). In this paper we obtain an improved condition number with O(N6). The advantages with respect to propagation of rounding errors are demonstrated.For GJGM, the resulting system from the equation-u(6)=f(x) is Ba=f*, where the matrix B is a diagonal matrix whose elements are bkk=(k+1)6. Thus we note that the condition number of the matrix B behaves like O(k6) for large values of k. Moreover, if we add ∑q=16η6-qG6-q to the matrix B, then we find that the eigenvalues of matrix D=B+∑q=16η6-qG6-q(η6-q=1,(1≤q≤6)) are all real positive and the effect of these additions does not significantly change the values of the condition number for the system. This means that matrix B, which resulted from the highest derivatives of the differential equations under investigation, play the most important role in the propagation of the roundoff errors. In Table 1 we list the values of the condition numbers for the two matrices B and D.Table 1
Condition numbers for the two matricesB and D.
NCond(B)Cond(D)Cond(B)/N6Cond(D)/N61680087627.894.773·10-44.547·10-4329061928607398.439·10-48.016·10-4647.497·1077.117·1071.091·10-31.035·10-31285.424·1095.148·1091.233·10-31.171·10-32563.685·10113.497·10111.309·10-31.243·10-3


### 3.4. Nonhomogeneous Boundary Conditions

Let us consider the sixth-order differential equation,(3.32)-u(6)(x)+∑q=16δqη6-qu(6-q)(x)=f(x),x∈(-1,1),
subject to the nonhomogeneous boundary conditions(3.33)u(j)(±1)=α±j,j=0,1,2.
In such case we can proceed as.Set(3.34)V(x)=u(x)+∑i=05cixi,
where ci,i=0,1,…,5 are coefficients to be determined such that V(x) satisfies the homogeneous boundary conditions, namely:  (3.35)V(j)(±1)=0,j=0,1,2.
Therefore the set of coefficients {ci,i=0,…,5} are determined by solving the following system of six equations:  (3.36)∑i=j5(±1)i-j(i-j+1)jci=-α±j,j=0,1,2.
The system of (3.36) is equivalent to the following matrix equation:  (3.37)Mc=-α,
where c=(c0,c1,…,c5)T is the vector of unknowns, α=(α+0,α-0,…,α+5,α-5)T, and M=(mkj)0≤k,j≤5 is a nonsingular matrix of order six and is given explicitly by (3.38)(1111111-11-11-101234501-23-4500261220002-612-20).
The transformation (3.34) turns the nonhomogeneous boundary conditions (3.36) into the homogeneous boundary conditions (3.35). Hence it suffices to solve the following modified one-dimensional sixth-order differential equation:  (3.39)-V(6)(x)+∑q=16δqη6-qV(6-q)(x)=f*(x),x∈(-1,1),
subject to the homogeneous boundary conditions (3.35), and
(3.40)f*(x)=f(x)+∑q=16δqη6-q∑j=6-q5cj(j+q-5)6-qxj+q-6.
Now, with the aid of the relation (see, [25])
(3.41)xl=∑j=0lDl,jRj(3,3)(x),Dl,j=2jl!(2j+7)(j+1)3n!∑i=0l-j(-2)i(j+4)ii!(l-i-j)!(j+7)i+1(i+j+8)j,
we can write f*(x) as (3.42)f*(x)=f(x)+∑q=16δqη6-q∑j=6-q5cj(j+q-5)6-q∑k=0j+q-6Dj+q-6,kRk(3,3)(x).
If we apply Galerkin method to the modified equation (3.39), we get the following system of equations:(3.43)(B+∑q=16η6-qG6-q)a=f*,
where f*=(f0*,f1*,…,fN-6*)T; fk*=(f*(x),Rk(3,3)(x))w/hk3,3, and the nonzero elements of the matrices B and G2n-q,1≤q≤6 are given explicitly as in Theorem 3.5.


## 3.1. Homogeneous Boundary Conditions

Let us consider the sixth-order differential equation(3.1)-u(6)(x)+∑q=16δqη6-qu(6-q)(x)=f(x),x∈(-1,1),
subject to the homogeneous boundary conditions(3.2)u(j)(±1)=0,j=0,1,2,
where u(j)(x) denotes the jth derivative of u(x) with respect to x and {η6-q,q=1,…,6} are positive constants, and(3.3)δq={(-1)1+q/2,qeven,1,qodd.
Let us denote Hwr(I)(r=0,1,2,…), as the weighted Sobolev spaces, whose inner products and norms are denoted by (·,·)r,w and∥·∥r,w, respectively. To account for homogeneous boundary conditions, we define (3.4)H0,w3(I)={v∈Hw3(I):v(j)(±1)=0,0≤j≤2},
where v(j)(x)=(djv)/(dxj). The superscript w will be omitted in case of w=1.LetPN be the space of all polynomials of degree less than or equal to N. Setting VN=PN∩H03(I), then (3.5)VN:=span⁡{J6(-3,-3)(x),J7(-3,-3)(x),…,JN(-3,-3)(x)}.
The symmetric generalized Jacobi-Galerkin procedure for solving (3.1)-(3.2) is to find uN∈VN such that(3.6)(-D6uN(x),v(x))+∑q=16δqη6-q(D6-quN(x),v(x))=(f(x),v(x)),∀v∈VN,
where (u,v)=∫-11uvdx is the scalar inner product in the space L2(-1,1).


## 3.2. The Choice of Basis Functions

We choose the basis functions of expansion to be(3.7)ϕk(x)=Jk+6(-3,-3)(x)=(1-x2)3Rk(3,3)(x),k=0,1,…,N-6,
which fulfills the boundary conditions (3.2).It is obvious that{ϕk(x)} are linearly independent. Therefore, we have(3.8)VN=span⁡{ϕk(x):k=0,1,2,…,N-6}.
Now, the following two lemmas are needed hereafter.Lemma 3.1.
For allk≥0, one has
(3.9)-D6[Jk+6(-3,-3)(x)]=(k+1)6Rk(3,3)(x).Proof.
Settingα=β=3 in relation to (2.7), we get
(3.10)(1-x2)Rk(3,3)(x)=62k+7[Rk(2,2)(x)-Rk+2(2,2)(x)].
Making use of this relation and with the aid of Lemma 2.4, we obtain
(3.11)-D6[Jk+6(-3,-3)(x)]=62k+7D2[(k+3)4Rk+2(2,2)(x)-(k+1)4Rk(2,2)(x)].
The last relation with the aid of the two relations (2.8) and (2.9) yields
(3.12)-D6[Jk+6(-3,-3)(x)]=(k+1)6Rk(3,3)(x).
This completes the proof of Lemma 3.1.Lemma 3.2.
For allk≥0, one has
(3.13)D6-q[Jk+6(-3,-3)(x)]=∑j=0qdj,k,qRk-2j+q(3,3),1≤q≤6,
where
(3.14)dj,k,q=(-1)j+1(qj)Γ(k-j+(7/2))(k+q-2j+6)!(k+q-2j+(7/2))2q(k+q-2j)!Γ(k+q-j+(9/2)).Proof.
Integrating formula (3.9) q times, q∈{1,2,…,6}, and with the aid of relation (2.11) (in case of α=3), we obtain the desired formula.Based on the results of the two Lemmas3.1 and 3.2, we are able to state and prove the following two theorems.Theorem 3.3.
One has, for arbitrary constantsak,
(3.15)-D6[∑k=0N-6akJk+6(-3,-3)(x)]=∑k=0N-6bkRk(3,3)(x),
where
(3.16)bk=(k+1)6ak.Theorem 3.4.
One has, for arbitrary constantsak, and  1≤q≤6,
(3.17)D6-q[∑k=0N-6akJk+6(-3,-3)(x)]=∑k=0N+q-6rk,qRk(3,3)(x),where(3.18)rk,q=∑i=0qdq-i,k+q-2i,qak+q-2i.Now, the application of Galerkin method to (3.1) gives  (3.19)(-D6uN(x)+∑q=16δqη6-qD6-quN(x),ϕk(x))=(f(x),ϕk(x)),where(3.20)uN(x)=∑k=0N-6akϕk(x),ϕk(x)=Jk+6(-3,-3)(x),k=0,1,…,N-6.
The variational formulation (3.19) is equivalent to(3.21)(-D6uN(x),Rk(3,3)(x))w3,3(x)+∑q=16δqη6-q(D6-quN(x),Rk(3,3)(x))w3,3(x)=(f(x),Rk(3,3)(x))w3,3(x),w3,3(x)=(1-x2)3.
Substitution of formulae (3.15) and (3.17) into (3.21) yields (3.22)(∑j=0N-6bjRj(3,3)(x)+∑q=16∑j=0N+q-6δqη6-qrj,qRj(3,3)(x),Rk(3,3)(x))w3,3(x)=(f(x),Rk(3,3)(x))w3,3(x),
where bk and rk,q are as given by (3.16) and (3.18), respectively.Now, if we apply the orthogonality relation ofRk(3,3)(x) on (3.22), we obtain the following linear system of equations:  (3.23)(bk+∑q=16δqη6-qrk,q)hk3,3=fk;k=0,1,…,N-6,
where (3.24)fk=(f(x),Rk(3,3)(x))w3,3(x),(3.25)hk3,3=4608(k+1)6(2k+7).
The linear system (3.23) may be put in the form(3.26)(bk+∑q=16δqη6-qrk,q)=fk*,k=0,1,…,N-6,
where(3.27)fk*=fkhk3,3.
This system of equations may be put in the matrix form(3.28)(B+∑q=16η6-qG6-q)a=f*,
where (3.29)a=(a0,a1,…,aN-6)T,f*=(f0*,f1*,…,fN-6*)T,
and the nonzero elements of the matrices B and G6-q,1≤q≤6 are given explicitly in the following theorem.Theorem 3.5.
IfuN(x)=∑k=0N-6akJk+6(-3,-3)(x) is the symmetric generalized Jacobi-Galerkin approximation to (3.1)-(3.2), then the expansion coefficients {ak:k=0,1,…,N-6} satisfy the matrix system (3.28), where the nonzero elements of the matrices  B=(bkj) and G6-q=(gkj6-q)=(δqrkj6-q), 0≤k,j≤N-6, 1≤q≤6 are given as:
(3.30)bkk=(k+1)6,gk,k+q-2i6-q=δq(-1)q+i+1(2k+7)(k+6)!q!Γ(-i+k+(7/2))2q+1i!k!(q-i)!Γ(k-i+q+(9/2)),0≤i≤q.It is worthy to note here that the case corresponding toη6-q=0, 1≤q≤6 leads to a linear system with diagonal matrix. The result for such case is summarized in the following important corollary.Corollary 3.6.
IfuN(x)=∑k=0N-6akJk+6(-3,-3)(x) and  η6-q=0, 1≤q≤6, is the symmetric generalized Jacobi-Galerkin approximation to problem (3.1)-(3.2), then the expansion coefficients {ak:k=0,1,…,N-6} are given explicitly by
(3.31)ak=2k+74608∫-11(1-x2)3f(x)Rk(3,3)(x),k=0,1,…,N-6.Remark 3.7.
Each of the matricesG6-q, q=1,2,…,6 in (3.28) is a band matrix whose total number of nonzero diagonals upper or lower the main diagonal does not exceed q. Thus the coefficient matrix D=B+∑q=16η6-qG6-q is seven-band matrix at most. This special structure of D simplifies greatly the solution of the linear system (3.28). The system in such case can be factorized by LU-decomposition and the number of operations necessary to construct this factorization is of order 78(N-5), and the number of operations needed to solve the two triangular systems is of order 25(N-5).Note
The total number of operations mentioned in the previous discussion includes the number of all subtractions, additions, divisions, and multiplications (see [24]).


## 3.3. Condition Number

For the direct collocation method, the condition number behaves likeO(N12) (N: maximal degree of polynomials). In this paper we obtain an improved condition number with O(N6). The advantages with respect to propagation of rounding errors are demonstrated.For GJGM, the resulting system from the equation-u(6)=f(x) is Ba=f*, where the matrix B is a diagonal matrix whose elements are bkk=(k+1)6. Thus we note that the condition number of the matrix B behaves like O(k6) for large values of k. Moreover, if we add ∑q=16η6-qG6-q to the matrix B, then we find that the eigenvalues of matrix D=B+∑q=16η6-qG6-q(η6-q=1,(1≤q≤6)) are all real positive and the effect of these additions does not significantly change the values of the condition number for the system. This means that matrix B, which resulted from the highest derivatives of the differential equations under investigation, play the most important role in the propagation of the roundoff errors. In Table 1 we list the values of the condition numbers for the two matrices B and D.Table 1
Condition numbers for the two matricesB and D.
NCond(B)Cond(D)Cond(B)/N6Cond(D)/N61680087627.894.773·10-44.547·10-4329061928607398.439·10-48.016·10-4647.497·1077.117·1071.091·10-31.035·10-31285.424·1095.148·1091.233·10-31.171·10-32563.685·10113.497·10111.309·10-31.243·10-3


## 3.4. Nonhomogeneous Boundary Conditions

Let us consider the sixth-order differential equation,(3.32)-u(6)(x)+∑q=16δqη6-qu(6-q)(x)=f(x),x∈(-1,1),
subject to the nonhomogeneous boundary conditions(3.33)u(j)(±1)=α±j,j=0,1,2.
In such case we can proceed as.Set(3.34)V(x)=u(x)+∑i=05cixi,
where ci,i=0,1,…,5 are coefficients to be determined such that V(x) satisfies the homogeneous boundary conditions, namely:  (3.35)V(j)(±1)=0,j=0,1,2.
Therefore the set of coefficients {ci,i=0,…,5} are determined by solving the following system of six equations:  (3.36)∑i=j5(±1)i-j(i-j+1)jci=-α±j,j=0,1,2.
The system of (3.36) is equivalent to the following matrix equation:  (3.37)Mc=-α,
where c=(c0,c1,…,c5)T is the vector of unknowns, α=(α+0,α-0,…,α+5,α-5)T, and M=(mkj)0≤k,j≤5 is a nonsingular matrix of order six and is given explicitly by (3.38)(1111111-11-11-101234501-23-4500261220002-612-20).
The transformation (3.34) turns the nonhomogeneous boundary conditions (3.36) into the homogeneous boundary conditions (3.35). Hence it suffices to solve the following modified one-dimensional sixth-order differential equation:  (3.39)-V(6)(x)+∑q=16δqη6-qV(6-q)(x)=f*(x),x∈(-1,1),
subject to the homogeneous boundary conditions (3.35), and
(3.40)f*(x)=f(x)+∑q=16δqη6-q∑j=6-q5cj(j+q-5)6-qxj+q-6.
Now, with the aid of the relation (see, [25])
(3.41)xl=∑j=0lDl,jRj(3,3)(x),Dl,j=2jl!(2j+7)(j+1)3n!∑i=0l-j(-2)i(j+4)ii!(l-i-j)!(j+7)i+1(i+j+8)j,
we can write f*(x) as (3.42)f*(x)=f(x)+∑q=16δqη6-q∑j=6-q5cj(j+q-5)6-q∑k=0j+q-6Dj+q-6,kRk(3,3)(x).
If we apply Galerkin method to the modified equation (3.39), we get the following system of equations:(3.43)(B+∑q=16η6-qG6-q)a=f*,
where f*=(f0*,f1*,…,fN-6*)T; fk*=(f*(x),Rk(3,3)(x))w/hk3,3, and the nonzero elements of the matrices B and G2n-q,1≤q≤6 are given explicitly as in Theorem 3.5.


## 4. Two-Dimensional Sixth-Order Equations

In this section, we consider the basis functionsϕk(x) as defined in (3.7) to solve numerically the two-dimensional even-order equations(4.1)-Δ3u(x,y)+∑r=02γr(-Δ)ru(x,y)=f(x,y),inΩ,
subject to the homogeneous boundary conditions  (4.2)∂i∂xiu(±1,y)=∂i∂yiu(x,±1)=0,i=0,1,2,
where Ω=(-1,1)×(-1,1), the differential operator Δ is the well-known Laplacian defined by Δ≡∂2/∂x2+∂2/∂y2, and γr,1≤r≤n-1 are constant, and f(x,y) is a given function.The symmetric generalized Jacobi-Galerkin approximation to (4.1)-(4.2) is to find uN∈VN2 such that  (4.3)(-Δ3uN,v)+∑r=02γr((-Δ)ruN,v)=(f,v),∀v∈VN2.
It is clear that if we take ϕk(x) as defined in (3.7), then(4.4)VN2=span⁡{ϕi(x)ϕj(y),i,j=0,1,…,N-6}.
Let us denote(4.5)uN=∑k=0N-6∑j=0N-6ukjϕk(x)ϕj(y),fkj*=1hk3,3hj3,3(f(x,y),Rk(3,3)(x)Rj(3,3)(y))w(x,y),(4.6)U=(ukj),F*=(fkj*),k,j=0,1,…,N-6,w(x,y)=w(3,3)(x)w(3,3)(y)=(1-x2)3(1-y2)3.
Taking v(x,y)=ϕℓ(x)ϕm(y) in (4.3) for ℓ,m=0,1,…,N-6, then we find that (4.3) is equivalent to the following equation:  (4.7)∑l,m=0N-6{{bilulmgjm0+3gil4ulmgjm2+3gil2ulmgjm4+glm0ulmbjm}+γ2{gil4ulmgjm0+2gil2ulmgjm2+gil0ulmgjm4}+γ1{gil2ulmgjm0+gil0ulmgjm2}+γ0gil0ulmgjm0}=fij*,i,j=0,1…,N-6,
which may be written in the matrix form: (4.8)BUG0T+3G4UG2T+3G2UG4T+G0UBT+γ2{G4UG0T+2G2UG2T+G0UG4T}+γ1{G2UG0T+G0UG2T}+γ0G0UG0T=F*,
where U and F* are as defined in (4.6) and the nonzero elements of the matrices B and G6-q,1≤q≤6, are those given as in Theorem 3.5.We can also rewrite (4.7) in the following form using the Kronecker matrix algebra (see, [26]):  (4.9)Lv≡[B⊗G0+3G4⊗G2+3G2⊗G4+G0⊗B+γ2{G4⊗G0+2G2⊗G2+G0⊗G4}+γ1{G2⊗G0+G0⊗G2}+γ0G0⊗G0]v=f*,
where f and v are F* and U written in a column vector, that is,(4.10)f*=(f00*,f10*,…,fN-6,0*;f01*,f11*,…,fN-6,1*;…;f0,N-6*,…,fN-6,N-6*)T,v=(u00,u10,…,uN-6,0;u01,u11,…,uN-6,1;…;u0,N-6,…,uN-6,N-6)T,
and ⊗ denotes the tensor product of matrices, that is, A⊗B=(Abij)i,j=0,1,…,N-6.A good review for the properties of the Kronecker product can be found in Graham [26] and Horn and Johnson [27].In summary, the solution of (4.1)-(4.2) consists of the following six steps.  (i)
Compute the matricesF*,B,G0,G2, and G4.(ii)
Compute the tensor products which appear in (4.9).(iii)
WriteF* in a column vector f*.(iv)
Obtain a column vectorv by solving (4.9).(v)
Rewrite a column vectorv in the form U.(vi)
FinduN.Remark 4.1.
SinceB is a diagonal matrix and each of the matrices G6-q,1≤q≤6, is seven-band at most, so the matrix L in system (4.9) is 6(N-4)-band at most, thus this system can be factorized by LU-decomposition and the number of operations necessary to construct this factorization is of order 6(N-5)(N-4)[12(N-4)+1], and the number of operations needed to solve the two triangular systems is of order (N-5)[24(N-4)-3].


## 5. Numerical Results

We consider here three different examples.Example 5.1.
Consider the following one-dimensional sixth-order equation:(5.1)-u(6)(x)+∑q=16δqη6-qu(6-q)(x)=f(x),x∈(-1,1),u(j)(±1)=0,j=0,1,2,
where f(x) is chosen such that the exact solution of (5.1) is u(x)=(1-x2)3cosx. The approximate spectral solution of (5.1) is given by
(5.2)uN(x)=∑k=0N-6akJk+6(-3,-3)(x),
and the vector of unknowns a=(a0,a1,…,aN-6)T is the solution of the system
(5.3)(B+∑q=16η6-qG6-q)a=f*,
where the nonzero elements of the matrices B and G6-q, 1≤q≤6, are given explicitly as in Theorem 3.5 and fk*=(f(x),Rk(3,3)(x))w/hk3,3.
Table2 lists the maximum pointwise error E for u-uN to (5.1), using SGJGM for various values N and the set of coefficients {ηi,0≤i≤5}.Table 2
Maximum pointwise error ofu-uN for N=8,12,16,20.
NηiEηiEηiE806.661·10-419.391·10-48i+12.391·10-3122.791·10-82.272·10-1812i+11.674·10-7161.342·10-131.269·10-1316i+11.241·10-12202.512·10-163.703·10-1620i+15.490·10-16Example 5.2.
Consider the following BVP (see [28]):
(5.4)y(6)(x)+y(x)=6(2xcos⁡x+5sin⁡x),x∈[-1,1],y(-1)=y(1)=0,y(1)(-1)=y(1)(1)=2sin⁡(1),y(2)(-1)=-y(2)(1)=-4cos⁡(1)-2sin⁡(1).
The exact solution of the above problem is
(5.5)y(x)=(x2-1)sin⁡x.
Table 3 lists the maximum pointwise error E=u-uN using SGJGM for various valuesof N. This table shows that the best accuracy obtained by our method is [2.256·10-16 for N=16] which is much better than the best accuracy obtained in Akram and Siddiqi [28] [3.81·10-16].Table 3
Maximum pointwise error ofu-uN for N=8,10,12,14,16.
NE88.301·10-6102.247·10-8124.499·10-11147.707·10-14162.256·10-16Example 5.3.
Consider the two-dimensional sixth-order equation:(5.6)-Δ3u+γ2Δ2u-γ1Δu+γ0u=f(x,y),∂iu∂xi(±1,y)=∂iu∂yi(x,±1)=0,i=0,1,2,
where  f(x,y) is chosen such that the exact solution of (5.6) is
(5.7)u(x,y)=(1-x2)(1-y2)sin⁡2⁡(2πx)sin⁡2⁡(2πy).
In Table 4, we list the maximum pointwise errors of  u-uN, using SGJGM with various choices of N.Table 4
Maximum pointwise error ofu-uN for N=20,30,40.
Nγ0γ1γ2SGJGMγ0γ1γ2SGJGM203111.18·10-30001.18·10-3302.39·10-91.18·10-9402.11·10-142.08·10-14


## 6. Concluding Remarks

We have presented some efficient direct solvers for sixth-order equations in one- and two-dimensions using the symmetric generalized Jacobi-Galerkin method. The algorithms are very efficient. In particular, we have found that, for some particular differential equations, the resulting systems of linear equations are diagonal. This, of course greatly simplify the numerical computations for these special cases. The use of symmetric generalized Jacobi polynomials leads to simplified analysis and very efficient numerical algorithms. Numerical results are presented which exhibit the high accuracy of the proposed algorithms.



---

*Source: 749370-2012-05-22.xml*