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

Abstract

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 (𝑁−5) unknowns. An efficient and accurate direct solution for algorithms based on the Legendre-Galerkin 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.

1. 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 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 finite-difference 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 boundary-value 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 sixth-order 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 high-order 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 𝑂(𝑁) operations for a one-dimensional domain with (𝑁−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.

2. One-Dimensional Sixth-Order Equations with Polynomial Coefficients

We first introduce some basic notation which will be used in the sequel. We denote by 𝐿_𝑛(𝑥) the 𝑛th degree Legendre polynomial, and we set 𝑆_𝑁=span𝐿₀(𝑥),𝐿₁(𝑥),…,𝐿_𝑁(𝑥),𝑊_𝑁=𝑣∈𝑆_𝑁∶𝑣^(𝑗)(±1)=0,𝑗=0,1,2, (2.1) where 𝑣^(𝑗)(𝑥) denotes 𝑗th-order differentiation of 𝑣(𝑥) with respect to 𝑥.

We recall that the {𝐿_𝑛(𝑥)} satisfy the orthogonality relation¹₋₁𝐿_𝑚(𝑥)𝐿_𝑛(𝑥)ℎ𝑑𝑥=0,𝑚≠𝑛,_𝑛ℎ,𝑚=𝑛,^𝑛=22𝑛+1,∀𝑚,𝑛≥0. (2.2) We recall also that 𝐿_𝑛(𝑥) is a polynomial of degree 𝑛 and therefore 𝐿_𝑛^(𝑞)(𝑥)∈𝑆_𝑛−𝑞. The following relation (the 𝑞th derivative of 𝐿_𝑘(𝑥)) will be needed for our main results (see Doha [23])𝐷^𝑞𝐿_𝑘(𝑥)=^k−𝑞_{()𝑖=0𝑘+𝑖𝑒𝑣𝑒𝑛}𝐶_𝑞(𝑘,𝑖)𝐿_𝑖(𝑥), (2.3) where𝐶_𝑞()=2𝑘,𝑖^𝑞−1(2𝑖+1)Γ(1/2)𝑞+𝑘−𝑖Γ(1/2)𝑞+𝑘+𝑖+1Γ𝑞Γ)Γ).(1/22−𝑞+𝑘−𝑖(1/23−𝑞+𝑘+𝑖 (2.4)

Some other useful relations are𝐿_𝑛(±1)=(±1)^𝑛,𝐿_𝑛^(𝑞)(±1)=(±1)^𝑛+𝑞𝑛+𝑞!2^𝑞𝑛−𝑞!𝑞!. (2.5)

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:−𝜒₃(𝑥)𝑢⁽⁶⁾(𝑥)+²_𝑖=0(−1)^𝑖𝛼_𝑖𝜒_𝑖(𝑥)𝑢^(2𝑖)(𝑥)=𝑓(𝑥),𝑥∈[−1,1], (2.6) subject to𝑢^(𝑗)(±1)=0,𝑗=0,1,2, (2.7) where 𝛼_𝑖, 𝑖=0,1,2 are constants and 𝜒_𝑖(𝑥), 𝑖=0,1,2,3 are given polynomials. Moreover, 𝑓(𝑥) is a given source function. Without loss of generality, we suppose that 𝜒₃(𝑥)=𝑥^𝜇,𝜒_𝑖(𝑥)=𝑥^𝜈, and 𝜒₀(𝑥)=𝑥^𝜎 where 𝜇, 𝜈, and 𝜎 are positive integers.

2.1. 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 𝜙_𝑘(𝑥) of the form𝜙_𝑘(𝑥)=³_𝑗=0𝜂_𝑗(𝑘)𝐿_𝑘+2𝑗(𝑥);𝜂₀(𝑘)=1,𝑘=0,1,…,𝑁−6. (2.8) We will choose the coefficients {𝜂_𝑗(𝑘)} such that 𝜙_𝑘(𝑥) verifies the boundary conditions (2.7). Making use of (2.5) and (2.8), hence {𝜂_𝑗(𝑘)} can be uniquely determined to obtain𝜂_𝑗(𝑘)=(−1)^𝑗32𝑘+4𝑗+1(𝑘+3/2)₃3−𝑗!𝑗!𝑘+𝑗+1/2₄,𝑗=1,2,3, (2.9) where (𝑎)_𝑏=Γ(𝑎+𝑏)/Γ(𝑎). Now, substitution of (2.9) into (2.8) yields 𝜙_𝑘(𝑥)=𝐿_𝑘(𝑥)+³_𝑗=1(−1)^𝑗32𝑘+4𝑗+1(𝑘+3/2)₃3−𝑗!𝑗!𝑘+𝑗+1/2₄𝐿_𝑘+2𝑗(𝑥),𝑘=0,1,2,…,𝑁−6. (2.10) It is obvious that {𝜙_𝑘(𝑥)} are linearly independent. Therefore by dimension argument we have𝑊_𝑁=span𝜙_𝑘(𝑥)∶𝑘=0,1,2,…,𝑁−6. (2.11)

2.2. 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𝑥^𝑚𝐿_𝑗(𝑥)=^2𝑚_𝑛=0Θ_𝑚,𝑛𝑗𝐿_{𝑗+𝑚−𝑛}(𝑥),∀𝑚,𝑗≥0, (2.12) with 𝐿_−𝑟(𝑥)=0,𝑟≥1, whereΘ_𝑚,𝑛𝑗=(−1)^𝑛2^{𝑗+𝑚−𝑛}𝑚!2𝑗+2𝑚−2𝑛+1Γ𝑗+1Γ×𝑗+𝑚−𝑛+1_{min(𝑗+𝑚−𝑛,𝑗)𝑘=max0,𝑗−𝑛}𝑘𝑗+𝑚−𝑛Γ𝑗+𝑘+12^𝑘×𝑛+𝑘−𝑗!Γ3𝑗+2𝑚−2𝑛−𝑘+2_{𝑗−𝑘ℓ=0}(−1)^ℓΓ2𝑗+𝑚−𝑛−𝑘−ℓ+1Γ𝑗+𝑚+ℓ−𝑛+1Γ()×ℓ!𝑗−𝑘−ℓ!Γ𝑗−ℓ+1𝑘+ℓ+1₂𝐹₁𝑗−𝑘−𝑛,𝑗+𝑚+ℓ−𝑛+1;3𝑗+2𝑚−2𝑛−𝑘+2;2. (2.13) 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.

Lemma 2.1. We have, for arbitrary constants ℓ,𝑗≥0, 𝑥^ℓ𝜙_𝑗^(𝑞)(𝑥)=³_𝑖=0^2ℓ_𝑛=0Θ_ℓ,𝑛𝑗+2𝑖𝐿_{𝑗+2𝑖+ℓ−𝑛}(𝑥)𝜂_𝑖𝑗,𝑥^ℓ𝜙_𝑗^(𝑞)(𝑥)=³_𝑖=0^{𝑗+2𝑖−𝑞}_𝑠=0^2ℓ_𝑛=0Θ_ℓ,𝑛(𝑠)𝐿_{𝑠+ℓ−𝑛}(𝑥)𝐶_𝑞𝑗+2𝑖,𝑠𝜂_𝑖𝑗, (2.14) where 𝜙_𝑗(𝑥), 𝜂_𝑖(𝑗), and Θ_ℓ,𝑛(𝑗) 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 𝑢_𝑁∈𝑊_𝑁 such that −𝑥^𝜇𝑢_𝑁⁽⁶⁾(𝑥),𝑣(𝑥)+²_𝑖=1𝛼_𝑖(−1)^𝑖𝑥^𝜈𝑢_𝑁^(2𝑖)(𝑥),𝑣(𝑥)+𝛼₀𝑥^𝜎𝑢_𝑁(𝑥),𝑣(𝑥)=𝑓(𝑥),𝑣(𝑥)_𝑁,∀𝑣∈𝑊_𝑁, (2.15) where (𝑢,𝑣)=∫¹₋₁𝑢𝑣𝑑𝑥 is the scalar product in 𝐿²(−1,1) and (⋅,⋅)_𝑁 is the inner product associated with the Legendre-Gauss-Lobatto quadrature. It is clear that if we take 𝜙_𝑘(𝑥) as defined in (2.8) and 𝑣(𝑥)=𝜙_𝑘(𝑥), then we find that (2.15) is equivalent to −𝑥^𝜇𝑢_𝑁⁽⁶⁾(𝑥),𝜙_𝑘(𝑥)+²_𝑖=1𝛼_𝑖(−1)^𝑖𝑥^𝜈𝑢_𝑁^(2𝑖)(𝑥),𝜙_𝑘(𝑥)+𝛼₀𝑥^𝜎𝑢_𝑁(𝑥),𝜙_𝑘(𝑥)=𝑓(𝑥),𝜙_𝑘(𝑥)_𝑁,𝑘=0,1,…,𝑁−6. (2.16)
Hence, by setting 𝑓_𝑘=𝑓,𝜙_𝑘(𝑥)_𝑁,𝐟=𝑓₀,𝑓₁,…,𝑓_𝑁−6^𝑇,𝑢_𝑁(𝑥)=^𝑁−6_𝑛=0𝑎_𝑛𝜙_𝑛(𝑥),𝐚=𝑎₀,𝑎₁,…,𝑎_𝑁−6^𝑇,𝑍^ℓ_𝑖=𝑧^𝑖,ℓ_{𝑘𝑗0≤𝑘,𝑗≤𝑁−6};0≤𝑖≤3,ℓisapositiveinteger, (2.17) where 𝑧^𝑟,ℓ_𝑘𝑗=(−1)^𝑟𝑥^ℓ𝜙_𝑗^(2𝑟)(𝑥),𝜙_𝑘(𝑥),𝑟=0,1,2,3, (2.18) then the matrix system associated with (2.16) becomes 𝑍^𝜇₃+²_𝑖=1𝛼_𝑖𝑍^𝜈_𝑖+𝛼₀𝑍^𝜎₀𝐚=𝐟, (2.19) where the nonzero elements of the matrices 𝑍^𝜇₃, 𝑍^𝜈₂, 𝑍^𝜈₁, and 𝑍^𝜎₀ are given explicitly in the following theorem.

Theorem 2.2. If we take 𝜙_𝑘(𝑥) as defined in (2.8), and if we denote 𝑧^3,𝜇_𝑘𝑗=(−𝑥^𝜇𝜙_𝑗⁽⁶⁾(𝑥),𝜙_𝑘(𝑥)), 𝑧^𝑟,𝜈_𝑘𝑗=((−1)^𝑟𝑥^𝜈𝜙_𝑗^(2𝑟)(𝑥),𝜙_𝑘(𝑥)), 𝑟=1,2 and 𝑧^0,𝜎_𝑘𝑗=(𝑥^𝜎𝜙_𝑗(𝑥),𝜙_𝑘(𝑥)) then the nonzero elements (𝑧^3,𝜇_𝑘𝑗), (𝑧^𝑟,𝜈_𝑘𝑗), (𝑧^0,𝜎_𝑘𝑗) for 0≤𝑘, 𝑗≤𝑁−6 are given by 𝑧^3,𝜇_𝑘𝑗=−³_ℓ=0³_𝑖=0^2𝜇_𝑠=0𝐶₆𝑗+2𝑖,𝑘+2ℓ−𝜇+𝑠Θ_𝜇,𝑠𝑘+2ℓ−𝜇+𝑠𝜂_𝑖𝑗𝜂_ℓ(𝑘)ℎ_𝑘+2ℓ,𝑧𝑗=𝑘+2𝑝−𝜇,𝑝=0,1,…,𝜇,_{𝑟,𝑣𝑘𝑗}=(−1)^𝑟3_ℓ=0³_𝑖=0^2𝑣_𝑠=0𝐶_2𝑟𝑗+2𝑖,𝑘+2ℓ−𝑣+𝑠Θ_𝑣,𝑠(𝑘+2ℓ−𝑣+𝑠)𝜂_𝑖𝑗𝜂_ℓ(𝑘)ℎ_𝑘+2ℓ,𝑧𝑗=𝑘+2𝑝−𝑣+2𝑟−6,𝑝=0,1,…,𝑣+6−2𝑟,𝑟=1,2,_{0,𝜎𝑘𝑗}=³_ℓ=0³_𝑖=0Θ_{𝜎,𝑗+2𝑖+𝜎−(𝑘+2ℓ)}𝑗+2𝑖𝜂_𝑖𝑗𝜂_ℓ(𝑘)ℎ_𝑘+2ℓ,𝑗=𝑘+2𝑝−𝜎−6,𝑝=0,1,…,𝜎+6. (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 𝑍^𝜇₃ is a band matrix with an upper bandwidth of 𝜇, lower bandwidth of 𝜇, and an overall bandwidth 2𝜇+1. The sparse matrices 𝑍^𝑣₂, 𝑍^𝑣₁, and 𝑍^𝜎₀ have bandwidths of 2𝑣+5, 2𝑣+9, and 2𝜎+13, respectively.

In general, the expense of calculating an LU factorization of an 𝑁×𝑁 dense matrix 𝐴 is 𝑂(𝑁³) operations, and the expense of solving 𝐴𝑥=𝑏, provided that the factorization is known, is 𝑂(𝑁²). However, in the case that a banded 𝐴 has bandwidth of 𝑟, we need just 𝑂(𝑟²𝑁) operations to factorize and 𝑂(𝑟𝑁) operations to solve a linear system. In the case of 𝛼_𝑖=0, 𝑖=0,1,2, the square matrix 𝑍^𝜇₃ has bandwidth of 2𝜇+1. We need just 𝑂((2𝜇+1)²𝑁) operations to factorize and 𝑂((2𝜇+1)𝑁) operations to solve the linear system (2.19). If 𝜇≪𝑁 this represents a very substantial saving. Notice also that the system (2.19) reduces to a diagonal system for 𝜇=0 and 𝛼_𝑖=0, 𝑖=0,1,2.

2.3. Constant Coefficients

In the special case, (𝜇=𝜈=𝜎=0, i.e., the sixth-order differential equation with constant coefficients), the corresponding matrix system becomes (𝑍⁰₃+∑³_𝑖=1𝛼_𝑖𝑍⁰_𝑖)𝐚=𝐟, where 𝑍⁰_𝑖=(𝑧^𝑖,0_𝑘𝑗)_{0≤𝑘,𝑗≤𝑁−6}; 𝑖=0,1,2,3.

Corollary 2.3. If 𝜇=𝜈=𝜎=0 then the nonzero elements (𝑧^3,0_𝑘𝑗), (𝑧^2,0_𝑘𝑗), (𝑧^1,0_𝑘𝑗), (𝑧^0,0_𝑘𝑗) for 0≤𝑘, 𝑗≤𝑁−6 are given as follows: 𝑧^3,0_𝑘𝑘=2⁶3𝑘+2₂3𝑘+2₃,𝑧^2,0_𝑘𝑗=³_ℓ=0³_𝑖=0𝜂_𝑖𝑗𝐶₄𝑗+2𝑖,𝑘+2ℓ𝜂_ℓ(𝑘)ℎ_𝑘+2ℓ𝑧,𝑗=𝑘+2𝑝−2,𝑝=0,1,2,_1,0𝑘𝑗=−³_ℓ=0³_𝑖=0𝜂_𝑖𝑗𝐶₂𝑗+2𝑖,𝑘+2ℓ𝜂_ℓ(𝑘)ℎ_𝑘+2ℓ𝑧,𝑗=𝑘+2𝑝−4,𝑝=0,1,2,3,4,_{0,0𝑘+2𝑝,𝑘}=𝑧^0,0_{𝑘,𝑘+2𝑝}=^3−𝑝_𝑖=0(−1)^2𝑖+𝑝3𝑖3()𝑖+𝑝2𝑘+4𝑝+4𝑖+1𝑘+3/2₃𝑘+2𝑝+3/2₃2𝑘+𝑝+𝑖+1/2₄𝑘+2𝑝+𝑖+1/2₄,𝑝=0,1,2,3. (2.21)(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 𝛼_𝑖=0, 𝑖=0,1,2, then the nonzero elements of the matrix 𝑍⁰₃ are given by (2.21) and the solution of the linear system is given explicitly by a_𝑘=(𝑓,𝜙_𝑘)_𝑁/𝑧^3,0_𝑘𝑘.

Obviously 𝑍⁰₂, 𝑍⁰₁, and 𝑍⁰₀ are symmetric positive definite matrices. Furthermore, 𝑍⁰₃ is a diagonal matrix, 𝑍⁰₂ can be split into two tridiagonal submatrices, 𝑍⁰₁ can be split into two pentadiagonal submatrices, and 𝐴₀ can be split into two sparse submatrices with bandwidth of 7. Therefore, the system can be efficiently solved. More precisely for 𝑘+𝑗 odd, 𝑧^𝑟,0_𝑘𝑗=0, 𝑟=0,1,2,3. Hence system (𝑍⁰₃+∑³_𝑖=1𝛼_𝑖𝑍⁰_𝑖)𝐚=𝐟 of order 𝑁−5 can be decoupled into two separate systems of order (𝑁/2−2) and (𝑁/2−3), respectively. In this way one needs to solve two systems of order 𝑛 instead of one of order 2𝑛, which leads to substantial savings. Moreover, in the case of 𝛼_𝑖≠0, 𝑖=0,1,2, we can form explicitly the LU factorization, that is, 𝑍⁰₃+∑³_𝑖=1𝛼_𝑖𝑍⁰_𝑖=𝐿𝑈. The special structure of 𝐿 and 𝑈 allows us to obtain the solution in 𝑂(𝑁) operations.

Remark 2.4. If the boundary conditions are nonhomogeneous, one can split the solution 𝑢(𝑥) into the sum of a low-degree polynomial which satisfies the nonhomogeneous boundary conditions plus a sum over the basis functions 𝜙(𝑥) that satisfy the equivalent homogeneous boundary conditions.

3. 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:−𝑋₃(𝑥)𝑌₃𝑦(Δ)³𝑢𝑥,𝑦+²_𝑟=1𝛼_𝑟𝑋_𝑟(𝑥)𝑌_𝑟𝑦(−Δ)^𝑟𝑢𝑥,𝑦+𝛼₀𝑋₀(𝑥)𝑌₀𝑦𝑢𝑥,𝑦=𝑓𝑥,𝑦,inΩ, (3.1) subject to the boundary conditions𝜕^𝑖𝜕𝑥^𝑖𝑢𝜕±1,𝑦=0,𝑖=0,1,2,_𝑖𝜕𝑦^𝑖𝑢(𝑥,±1)=0,𝑖=0,1,2, (3.2) where Ω=(−1,1)×(−1,1), the differential operator Δ is the well-known Laplacian defined by Δ=(𝜕²/𝜕𝑥²)+(𝜕²/𝜕𝑦²), (𝜕⁰/𝜕𝑥⁰)𝑢(±1,𝑦)=𝑢(±1,𝑦), (𝜕⁰/𝜕𝑥⁰)𝑢(𝑥,±1)=𝑢(𝑥,±1) and 𝑓(𝑥,𝑦) is a given source function. Moreover, 𝑋_𝑖(𝑥) and 𝑌_𝑖(𝑦), 𝑖=0,1,2,3 are given polynomials. Without loss of generality, we suppose that 𝑋₃(𝑥)=𝑥^𝜇, 𝑌₃(𝑦)=𝑦^𝜈, 𝑋_𝑖(𝑥)=𝑥^𝜌, 𝑌_𝑖(𝑦)=𝑦^𝜎, 𝑖=1,2, 𝑋₀(𝑥)=𝑥^𝛿, and 𝑌₀(𝑦)=𝑦^𝜀 where 𝜇, 𝜈, 𝜌, 𝜎, 𝛿, and 𝜀 are positive integers.

The Legendre-Galerkin approximation to (3.1)-(3.2) is, to find 𝑢_𝑁∈𝑊²_𝑁 such that𝑥^𝜇𝑦^𝜈(−Δ)³𝑢_𝑁,𝑣+²_𝑟=1𝛼_𝑟𝑥^𝜌𝑦^𝜎(−Δ)^𝑟𝑢_𝑁,𝑣+𝛼₀𝑥^𝛿𝑦^𝜀𝑢_𝑁,𝑣=𝑓,𝑣_𝑁∀𝑣∈𝑊²_𝑁. (3.3) It is clear that if we take 𝜙_𝑘(𝑥) as defined in (2.8), then𝑊²_𝑁=span𝜙_𝑖(𝑥)𝜙_𝑗𝑦,𝑖,𝑗=0,1,…,𝑁−6. (3.4)

We denote𝑢_𝑁=^𝑁−6_𝑘=0^𝑁−6_𝑗=0𝑢_𝑘𝑗𝜙_𝑘(𝑥)𝜙_𝑗𝑦,𝑓_𝑘𝑗=𝑓,𝜙_𝑘(𝑥)𝜙_𝑗(𝑦)_𝑁,𝑈=𝑢_𝑘𝑗,𝐹=𝑓_𝑘𝑗,𝑘,𝑗=0,1,…,𝑁−6. (3.5) Taking 𝑣(𝑥,𝑦)=𝜙_ℓ(𝑥)𝜙_𝑚(𝑦) in (3.3) for ℓ,𝑚=0,1,…,𝑁−6, then one can observe that (3.3) is equivalent to the following equation:^𝑁−6_ℓ,𝑚=0𝑧^3,𝜇_𝑖ℓ𝑢_ℓ𝑚𝑧^0,𝜈_𝑗𝑚+3𝑧^2,𝜇_𝑖ℓ𝑢_ℓ𝑚𝑧^1,𝜈_𝑗𝑚+3𝑧^1,𝜇_𝑖ℓ𝑢_ℓ𝑚𝑧^2,𝜈_𝑗𝑚+𝑧^0,𝜇_𝑖ℓ𝑢_ℓ𝑚𝑧^3,𝜈_𝑗𝑚+𝛼₂𝑧^2,𝜌_𝑖ℓ𝑢_ℓ𝑚𝑧^0,𝜎_𝑗𝑚+2𝑧^1,𝜌_𝑖ℓ𝑢_ℓ𝑚𝑧^1,𝜎_𝑗𝑚+𝑧^0,𝜌_𝑖ℓ𝑢_ℓ𝑚𝑧^2,𝜎_𝑗𝑚+𝛼₁𝑧^1,𝜌_𝑖ℓ𝑢_ℓ𝑚𝑧^0,𝜎_𝑗𝑚+𝑧^0,𝜌_𝑖ℓ𝑢_ℓ𝑚𝑧^1,𝜎_𝑗𝑚+𝛼₀𝑧^0,𝛿_𝑖ℓ𝑢_ℓ𝑚𝑧^0,𝜀_𝑗𝑚=𝑓_𝑖𝑗,𝑖,𝑗=0,1,…,𝑁−6, (3.6) which may be written in the matrix form𝑍^𝜇₃𝑈𝑍^𝜈₀+3𝑍^𝜇₂𝑈𝑍^𝜈₁+3𝑍^𝜇₁𝑈𝑍^𝜈₂+𝑍^𝜇₀𝑈𝑍^𝜈₃+𝛼₂𝑍^𝜌₂𝑈𝑍^𝜎₀+2𝑍^𝜌₁𝑈𝑍^𝜎₁+𝑍^𝜌₀𝑈𝑍^𝜎₂+𝛼₁𝑍^𝜌₁𝑈𝑍^𝜎₀+𝑍^𝜌₀𝑈𝑍^𝜎₁+𝛼₀𝑍^𝛿₀𝑈𝑍^𝜀₀=𝐹, (3.7) where {𝑍^𝜚_𝑖for𝑖=0,1,2,3and𝜚isapositiveinteger} 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 𝑛×𝑛 and 𝑚×𝑚, respectively. Their tensor product𝑃𝑃⊗𝑅=¹¹𝑅⋯𝑃_1𝑛𝑅𝑃⋮⋱⋮_𝑛1𝑅⋯𝑃_𝑛𝑛𝑅, (3.8) is a matrix of size 𝑚𝑛×𝑚𝑛.

We can also rewrite (3.7) in the following form using the Kronecker matrix algebra (See, Graham [24]):𝐿𝐯≡𝑍^𝜇₃⊗𝑍^𝜈₀+3𝑍^𝜇₂⊗𝑍^𝜈₁+3𝑍^𝜇₁⊗𝑍^𝜈₂+𝑍^𝜇₀⊗𝑍^𝜈₃+𝛼₂𝑍^𝜌₂⊗𝑍^𝜎₀+2𝑍^𝜌₁⊗𝑍^𝜎₁+𝑍^𝜌₀⊗𝑍^𝜎₂+𝛼₁𝑍^𝜌₁⊗𝑍^𝜎₀+𝑍^𝜌₀⊗𝑍^𝜎₁+𝛼₀𝑍^𝛿₀⊗𝑍^𝜀₀𝐯=𝐟, (3.9) where 𝐟 and 𝐯 are 𝐹 and 𝑈 written in a column vector, that is,𝐟=𝑓₀₀,𝑓₁₀,…,𝑓_𝑁−6,0;𝑓₀₁,𝑓₁₁,…,𝑓_𝑁−6,1;…;𝑓_0,𝑁−6,…,𝑓_{𝑁−6,𝑁−6}^𝑇,𝐯=𝑢₀₀,𝑢₁₀,…,𝑢_𝑁−6,0;𝑢₀₁,𝑢₁₁,…,𝑢_𝑁−6,1;…;𝑢_0,𝑁−6,…,𝑢_{𝑁−6,𝑁−6}^𝑇, (3.10) and ⊗ denotes the tensor product of matrices, that is, 𝑍^𝜇₃⊗𝑍^𝜈₀=(𝑍^𝜇₃𝑧^0,𝜈_𝑖𝑗)_{𝑖,𝑗=0,1,…,𝑁−6}. In brief, the solution of (3.1) subject to (3.2) can be summarized in Algorithm 1.

Compute the matrices {𝑍𝜇𝑗,𝑍𝜈𝑗for=0,1,2,3}, {𝑍𝜌𝑗,𝑍𝜎𝑗for𝑗=0,1,2} and {𝑍𝛿0,𝑍𝜀0} from Theorem 2.2. Find the tensor product of the matrices in the previous step. Compute 𝐹 and write it in a column vector 𝐟. Obtain a column vector 𝐯 by solving (3.9).

Algorithm 1 

4. 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 practice, since for a general right-hand side function 𝑓 we are not able to compute exactly its representation by Legendre polynomials. In fact, the so-called pseudospectral method is used to treat the right-hand side; that is, we replace 𝑓 by 𝐼_𝑁𝑓 (polynomial interpolation over the set of Gauss-Lobatto points); see for instance; Funaro [25].

Example 4.1. Here we present some numerical results for a one-dimensional sixth-order equation with polynomial coefficients. We only consider the case 𝜒₃(𝑥)=(𝑥⁴+𝑥+1), 𝜒₂(𝑥)=𝑥²+2𝑥+3, 𝜒₁(𝑥)=𝑥²+1 and 𝜒₀(𝑥)=1; that is, we consider the −𝑥⁴+1𝑢⁽⁶⁾(𝑥)+𝛾₂𝑥²+2𝑥+3𝑢⁽⁴⁾(𝑥)−𝛾₁𝑥²+1𝑢^′′(𝑥)+𝛾₀𝑢(𝑥)=𝑓(𝑥), (4.1) subject to 𝑢(±1)=𝑢′(±1)=𝑢^′′(±1)=0,𝑥∈𝐼, (4.2) where 𝑓(𝑥) is chosen such that the exact solution of (4.1) is 𝑢(𝑥)=(1−𝑥²)sin²(2𝜋𝑥).

For LGM, we have 𝑢_𝑁(𝑥)=∑^𝑁−6_𝑘=0𝑎_𝑘𝜙_𝑘(𝑥), the vectors of unknowns 𝐚 is the solution of the following system𝑍⁴₃+𝑍⁰₃+𝛾₂𝑍²₂+2𝑍¹₂+3𝑍⁰₂+𝛾₁𝑍²₁+𝑍⁰₁+𝛾₀𝑍⁰₀𝐚=𝐟, (4.3) where the nonzero elements of the matrices 𝑍⁴₃, 𝑍⁰₃, 𝑍²₂, 𝑍¹₂, 𝑍⁰₂, 𝑍²₁, 𝑍⁰₁ and 𝑍⁰₀ can be evaluated explicitly from Theorem 2.2. Table 1 lists the maximum pointwise error of 𝑢−𝑢_𝑁, using the LGM with various choices of 𝛾₀,𝛾₁,𝛾₂ and 𝑁. Numerical results of this problem show that the Legendre Galerkin method converge exponentially.

Table 1 

Maximum pointwise error of 𝑢−𝑢_𝑁 for 𝑁=20,30,40.

In order to examine the algorithm proposed in Section 3, we will consider a problem for a two-dimensional sixth-order elliptic differential equation with constant coefficients.

Example 4.2. Consider the two-dimensional sixth-order equation −Δ³𝑢+𝛾₂Δ²𝑢−𝛾₁Δ𝑢+𝛾₀𝑢=𝑓𝑥,𝑦, (4.4) subject to the boundary conditions 𝑢±1,𝑦=𝑢(𝑥,±1)=0,𝜕𝑢=𝜕𝑥±1,𝑦𝜕𝑢()𝜕𝜕𝑦𝑥,±1=0,₂𝑢𝜕𝑥²=𝜕±1,𝑦²𝑢𝜕𝑦²(𝑥,±1)=0, (4.5) where 𝑓(𝑥,𝑦) is chosen such that the exact solution of (4.4)-(4.5) is 𝑢𝑥,𝑦=1−𝑥²1−𝑦²sin²(2𝜋𝑥)sin²2𝜋𝑦. (4.6)

In Table 2, we list the maximum pointwise error of 𝑢−𝑢_𝑁 by the LGM with two choices of 𝛾₀, 𝛾₁, 𝛾₂ and various choices of 𝑁. The results indicate that the spectral accuracy is achieved and that the effect of roundoff errors is very limited.

Table 2 

Maximum pointwise error of 𝑢−𝑢_𝑁 for 𝑁=20,30,40.

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

Acknowledgments

The authors are very grateful to the referees for carefully reading the paper and for their comments and suggestions which have improved the paper. This paper was funded by the Deanship of Scientific Research (DSR), King Abdulaziz University, Jeddah. The authors, therefore, acknowledge with thanks DSR technical and financial support.