An Extension of the Legendre-Galerkin Method for Solving Sixth-Order Differential Equations with Variable Polynomial Coefficients

We extend the application of Legendre-Galerkin algorithms for sixth-order elliptic problems with constant coefficients to sixth-order elliptic equations with variable polynomial coefficients. The complexities of the algorithm are O(N) operations for a one-dimensional domain with N − 5 unknowns. An efficient and accurate direct solution for algorithms based on the LegendreGalerkin approximations developed for the two-dimensional sixth-order elliptic equations with variable coefficients relies upon a tensor product process. The proposed Legendre-Galerkin method for solving variable coefficients problem is more efficient than pseudospectral method. Numerical examples are considered aiming to demonstrate the validity and applicability of the proposed techniques.


Introduction
Spectral methods are preferable in numerical solutions of ordinary and partial differential equations due to its high-order accuracy whenever it works 1, 2 . Recently, renewed interest in the Galerkin technique has been prompted by the decisive work of Shen 3 , where new Legendre polynomial bases for which the matrices systems are sparse are introduced. We introduce a generalization of Shen's basis to numerically solve the sixth-order differential equations with variable polynomial coefficients.
Sixth-order boundary-value problems arise in astrophysics; the narrow convecting layers bounded by stable layers, which are believed to surround A-type stars, may be modeled by sixth-order boundary-value problems 4, 5 . Further discussion of the sixth-order boundary-value problems is given in 6 . The literature of numerical analysis contains little work on the solution of the sixth-order boundary-value problems 4,5,7,8 . Theorems that 2 Mathematical Problems in Engineering list conditions for the existence and uniqueness of solutions of such problems are thoroughly discussed in 9 , but no numerical methods are contained therein.
From the numerical point of view, Shen 3 , Doha and Bhrawy 10-12 , and Doha et al. 13 have constructed efficient spectral-Galerkin algorithms using compact combinations of orthogonal polynomials for solving elliptic equations of the second and fourth order with constant coefficients in various situations. Recently, the authors in 14, 15 and 16 have developed efficient Jacobi dual-Petrov-Galerkin and Jacobi-Gauss collocation methods for solving some odd-order differential equations. Moreover, the Bernstein polynomials have been applied for the numerical solution of high even-order differential equations see, 17, 18 . For sixth-order differential equations, Twizell and Boutayeb 5 developed finitedifference methods of order two, four, six, and eight for solving such problems. Siddiqi and Twizell 7 used sixth-degree splines, where spline values at the mid knots of the interpolation interval and the corresponding values of the even order derivatives were related through consistency relations. A sixth-degree B-spline functions is used to construct an approximate solution for sixth-order boundary-value problems see 19 . Moreover, Septic spline solutions of sixth-order boundary value problems are introduced in 20 . El-Gamel et al. 8 proposed Sinc-Galerkin method for the solutions of sixth-order boundaryvalue problems. In fact, the decomposition and modified domain decomposition methods to investigate solution of the sixth-order boundary-value problems are introduced in 21 . Recently, Bhrawy 22 developed a spectral Legendre-Galerkin method for solving sixthorder boundary-value problems with constant coefficients. In this work, we introduce an efficient direct solution algorithm to generalize the work in 3, 22 .
The main aim of this paper is to extend the application of Legendre-Galerkin method LGM to solve sixth-order elliptic differential equations with variable coefficients by using the expansion coefficients of the moments of the Legendre polynomials and their highorder derivatives. We present appropriate basis functions for the Legendre-Galerkin method applied to these equations. This leads to discrete systems with sparse matrices that can be efficiently inverted. The complexities of the algorithm is O N operations for a onedimensional domain with N−5 unknowns. The direct solution algorithms developed for the homogeneous problem in two-dimensions with constant and variable coefficients rely upon a tensor product process. Numerical results indicating the high accuracy and effectiveness of these algorithms are presented. This paper is organized as follows. In the next section, we discuss an algorithm for solving the one-dimensional sixth-order elliptic equations with variable polynomial coefficients. In Section 3, we extend our results of Sections 2 to the two-dimensional sixth-order equations with variable polynomial coefficients. In Section 4, we present two numerical examples to exhibit the accuracy and efficiency of the proposed numerical algorithms. Also a conclusion is given in Section 5.

One-Dimensional Sixth-Order Equations with Polynomial Coefficients
We first introduce some basic notation which will be used in the sequel. We denote by L n x the nth degree Legendre polynomial, and we set where v j x denotes jth-order differentiation of v x with respect to x. We recall also that L n x is a polynomial of degree n and therefore L q n x ∈ S n−q . The following relation the qth derivative of L k x will be needed for our main results see Doha 23 Some other useful relations are In this section, we are interested in using the Legendre-Galerkin method to solve the variable polynomial coefficients sixth-order differential equation in the form: where α i , i 0, 1, 2 are constants and χ i x , i 0, 1, 2, 3 are given polynomials. Moreover, f x is a given source function. Without loss of generality, we suppose that χ 3 x x μ , χ i x x ν , and χ 0 x x σ where μ, ν , and σ are positive integers.

Basis of Functions
The problem of approximating solutions of ordinary or partial differential equations by Galerkin approximation involves the projection onto the span of some appropriate set of basis functions, typically arising as the eigenfunctions of a singular Sturm-Liouville problem.
The members of the basis may satisfy automatically the boundary conditions imposed on the problem. As suggested in 3, 10-12 , one should choose compact combinations of orthogonal polynomials as basis functions to minimize the bandwidth and condition number of the resulting system. As a general rule, for one-dimensional sixth-order differential equations with six boundary conditions, one can choose the basis functions of expansion φ k x of the form We will choose the coefficients {η j k } such that φ k x verifies the boundary conditions 2.7 .
Making use of 2.5 and 2.8 , hence {η j k } can be uniquely determined to obtain where a b Γ a b /Γ a . Now, substitution of 2.9 into 2.8 yields It is obvious that {φ k x } are linearly independent. Therefore by dimension argument we have

Treatment of Variable Polynomial Coefficients
A more general situation which often arises in the numerical solution of differential equations with polynomial coefficients by using the Legendre Galerkin method is the evaluation of the expansion coefficients of the moments of high-order derivatives of infinitely differentiable functions. The formula of Legendre coefficients of the moments of one single Legendre polynomials of any degree is For more details about the above formula, the reader is referred to Doha 23 . This formula can be used to facilitate greatly the setting up of the algebraic systems to be obtained by applying the LGM for solving differential equations with polynomial coefficients of any order. The following lemma is very important and needed in what follows.

2.14
where φ j x , η i j , and Θ ,n j are as defined in 2.8 , 2.9 , and 2.13 , respectively.
Proof. Immediately obtained from relations 2.3 , 2.8 , and 2.12 , the Legendre-Galerkin approximation to 2.6 -2.7 is, to find u N ∈ W N such that where u, v 1 −1 uvdx is the scalar product in L 2 −1, 1 and ·, · N is the inner product associated with the Legendre-Gauss-Lobatto quadrature. It is clear that if we take φ k x as defined in 2. 8 and v x φ k x , then we find that 2.15 is equivalent to Hence, by setting n 0 a n φ n x , a a 0 , a 1 , . . . , a N−6 T , where the nonzero elements of the matrices Z μ 3 , Z ν 2 , Z ν 1 , and Z σ 0 are given explicitly in the following theorem.

2.20
Proof. The proof of this theorem is rather lengthy, but it is not difficult once Lemma 2.1 is applied.
From Theorem 2.2, we see that Z μ 3 is a band matrix with an upper bandwidth of μ, lower bandwidth of μ, and an overall bandwidth 2μ 1. The sparse matrices Z v 2 , Z v 1 , and Z σ 0 have bandwidths of 2v 5, 2v 9, and 2σ 13, respectively. In general, the expense of calculating an LU factorization of an N × N dense matrix A is O N 3 operations, and the expense of solving Ax b, provided that the factorization is known, is O N 2 . However, in the case that a banded A has bandwidth of r, we need just O r 2 N 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 Z μ 3 has bandwidth of 2μ 1. We need just O 2μ 1 2 N operations to factorize and O 2μ 1 N operations to solve the linear system 2.19 . If μ N this represents a very substantial saving. Notice also that the system 2.19 reduces to a diagonal system for μ 0 and α i 0, i 0, 1, 2.

Constant Coefficients
In the special case, μ ν σ 0, i.e., the sixth-order differential equation with constant coefficients , the corresponding matrix system becomes Z 0

2.22
Note that the results of Corollary 2.3 can be obtained immediately as a special case from Theorem 2.2. For more details see 22 . It is worthy to note here that if α i 0, i 0, 1, 2, then the nonzero elements of the matrix Z 0 3 are given by 2.21 and the solution of the linear system is given explicitly by a k f, φ k N /z 3,0 kk . Obviously Z 0 2 , Z 0 1 , and Z 0 0 are symmetric positive definite matrices. Furthermore, Z 0 3 is a diagonal matrix, Z 0 2 can be split into two tridiagonal submatrices, Z 0 1 can be split into two pentadiagonal submatrices, and A 0 can be split into two sparse submatrices with bandwidth of 7. Therefore, the system can be efficiently solved. More precisely for k j odd, z r,0 kj 0, r 0, 1, 2, 3. Hence system Z 0 i a f of order N − 5 can be decoupled into two separate systems of order N/2 − 2 and N/2 − 3 , respectively. In this way one needs to solve two systems of order n instead of one of order 2n, which leads to substantial savings. Moreover, in the case of α i / 0, i 0, 1, 2, we can form explicitly the LU factorization, that is,

Two-Dimensional Sixth-Order Equations with Polynomial Coefficients
In this section, we extend the results of Section 2 to deal with the two-dimensional sixth-order differential equations with variable polynomial coefficients: subject to the boundary conditions where Ω −1, 1 × −1, 1 , the differential operator Δ is the well-known Laplacian defined by and f x, y is a given source function. Moreover, X i x and Y i y , i 0, 1, 2, 3 are given polynomials. Without loss of generality, we suppose that x δ , and Y 0 y y ε where μ, ν, ρ, σ, δ, and ε are positive integers. The Legendre-Galerkin approximation to 3.1 -3.2 is, to find u N ∈ W 2 N such that

3.3
It is clear that if we take φ k x as defined in 2.8 , then We denote 3.5 Taking v x, y φ x φ m y in 3.3 for , m 0, 1, . . . , N − 6, then one can observe that 3.3 is equivalent to the following equation: which may be written in the matrix form where {Z i for i 0, 1, 2, 3 and is apositive integer} are the matrices defined in Theorem 2.2.
The direct solution algorithm here developed for the sixth-order elliptic differential equation in two dimensions relies upon a tensor product process, which is defined as follows. Let P and R be two matrices of size n × n and m × m, respectively. Their tensor product is a matrix of size mn × mn. We can also rewrite 3.7 in the following form using the Kronecker matrix algebra See, Graham 24 : where f and v are F and U written in a column vector, that is,

Numerical Results
We report on two numerical examples by using the algorithms presented in the previous sections. It is worthy to mention that the pure spectral-Galerkin technique is rarely used in Find the tensor product of the matrices in the previous step. Compute F and write it in a column vector f. Obtain a column vector v by solving 3.9 . 1.67 · 10 −9 1.67 · 10 −9 40 6.10 · 10 −16 6.10 · 10 −16 where f x, y is chosen such that the exact solution of 4.4 -4.5 is u x, y 1 − x 2 1 − y 2 sin 2 2πx sin 2 2πy . 4.6 In Table 2, we list the maximum pointwise error of u−u N by the LGM with two choices of γ 0 , γ 1 , γ 2 and various choices of N. The results indicate that the spectral accuracy is achieved and that the effect of roundoff errors is very limited.

Concluding Remarks
We have presented stable and efficient spectral Galerkin method using Legendre polynomials as basis functions for sixth-order elliptic differential equations in one and two dimensions. We concentrated on applying our algorithms to solve variable polynomials coefficients differential equations by using the expansion coefficients of the moments of the Legendre polynomials and their high-order derivatives. Numerical results are presented which exhibit the high accuracy of the proposed algorithms.