A Legendre tau-Spectral Method for Solving Time-Fractional Heat Equation with Nonlocal Conditions

Abstract

We develop the tau-spectral method to solve the time-fractional heat equation (T-FHE) with nonlocal condition. In order to achieve highly accurate solution of this problem, the operational matrix of fractional integration (described in the Riemann-Liouville sense) for shifted Legendre polynomials is investigated in conjunction with tau-spectral scheme and the Legendre operational polynomials are used as the base function. The main advantage in using the presented scheme is that it converts the T-FHE with nonlocal condition to a system of algebraic equations that simplifies the problem. For demonstrating the validity and applicability of the developed spectral scheme, two numerical examples are presented. The logarithmic graphs of the maximum absolute errors is presented to achieve the exponential convergence of the proposed method. Comparing between our spectral method and other methods ensures that our method is more accurate than those solved similar problem.

1. Introduction

In recent years, many engineering and physical phenomena can be successfully described by models of fractional differential equations (FDEs); see, for instance, [1–7]. Thus many researchers have been interested in studying the properties of fractional calculus and finding stable and robust numerical and analytical schemes for solving FDEs such as spectral tau method [8–10], Crank-Nicolson method [11], compact finite difference approximation [12], Legendre wavelets method [13], Haar wavelet operational matrix method [14], iterative Laplace transform method [15], Lie symmetry analysis method [16], and other methods [17–20].

Recently, spectral methods [21–23] have been applied to solve ordinary FDEs (see [24, 25]) while in [26, 27] the authors introduced the operational matrices of fractional derivatives with the help of the spectral methods to solve FDEs. This is not all; the partial FDEs are also investigated by using the spectral methods. In [28–31], the tau and collocation spectral methods are implemented in combination with the operational matrices of fractional integration for approximating the solution of some classes of space-fractional differential equations.

The T-FHE is a generalization of the classical heat equation obtained by replacing the first order time derivative by a fractional derivative of order 𝜈, 0<𝜈≤1. Ali and Jassim [32] used the homotopy perturbation method to solve the T-FHE, while in [33] the authors introduced a general iteration formula of variational iteration method for a solution of the T-FHE. Moreover, in [34] the differential transform method is applied to solve the T-FHE. In addition, Rostamy and Karimi [35] constructed the Bernstein operational matrix for the fractional derivatives and used it together with spectral method to solve the T-FHE.

In this paper, we consider the T-FHE with the nonlocal condition [36]: 𝜕^𝜈𝑢(𝑥,𝑡)𝜕𝑡^𝜈−𝜕²𝑢(𝑥,𝑡)𝜕𝑥²=𝑞(𝑥,𝑡),0<𝑥≤1,0<𝑡≤1, (1) subject to 𝑢()()(𝑥)𝑢()𝑥,0=𝑢𝑥,1+𝑓,0<𝑥≤1,0,𝑡=𝑔₀(𝑡),𝑢(1,𝑡)=𝑔₁(𝑡),0<𝑡≤1, (2) where 0<𝜈≤1, 𝑢(𝑥,𝑡) is the temperature as a function of space 𝑥 and time 𝑡, and 𝑞(𝑥,𝑡) is known source term. Our main aim is to achieve highly accurate solution of the T-FHE with nonlocal conditions (1) and (2). The tau-spectral method is applied based on the shifted Legendre polynomials as a basis function with the help of the operational matrix of fractional integration of such polynomials. Two numerical examples are introduced and solved using the presented technique to show its accuracy and validity. Also, we introduce comparisons between our numerical results and those obtained using the implicit difference approximation (IDA).

This paper is arranged in the following way: in Section 2 we introduce some definitions and notations of fractional calculus with some properties of Legendre polynomials. In Section 3 we apply our algorithm for the solution of the T-FHE with nonlocal condition. In Section 4 two numerical examples and comparisons between our results and those obtained by the IDA are introduced. Also in Section 5, a conclusion is presented.

2. Preliminaries and Notations

2.1. Fractional Calculus Definitions

Riemann-Liouville and Caputo fractional definitions are the two most used from other definitions of fractional derivatives which have been introduced recently.

Definition 1. The integral of order 𝛾≥0 (fractional) according to Riemann-Liouville is given by 𝐼^𝜈𝑓(𝑥)=1Γ(𝜈)^𝑥₀(𝑥−𝑡)^𝜈−1𝑓(𝑡)𝐼𝑑𝑡,𝜈>0,𝑥>0,₀𝑓(𝑥)=𝑓(𝑥), (3) where Γ(𝜈)=^∞₀𝑥^𝜈−1𝑒^−𝑥𝑑𝑥 (4) is gamma function.

The operator 𝐼^𝜈 satisfies the following properties: 𝐼^𝜈𝐼^𝜇𝑓(𝑥)=𝐼^𝜈+𝜇𝑓(𝑥),𝐼^𝜈𝐼^𝜇𝑓(𝑥)=𝐼^𝜇𝐼^𝜈𝑓(𝑥),𝐼^𝜈𝑥^𝛽=Γ𝛽+1Γ𝑥𝛽+1+𝜈^𝛽+𝜈. (5)

Definition 2. The Caputo fractional derivative of order 𝜈 is defined by 𝐷^𝜈𝑓(𝑥)=1Γ(𝑚−𝜈)^𝑥₀(𝑥−𝑡)^{𝑚−𝜈−1}𝑑^𝑚𝑑𝑡^𝑚𝑓(𝑡)𝑑𝑡,𝑚−1<𝜈≤𝑚,𝑥>0, (6) where 𝑚 is the ceiling function of 𝜈.

The operator 𝐷^𝛾 satisfies the following properties: 𝐷^𝜈()𝐼𝐶=0,𝐶isconstant_𝜈𝐷^𝜈𝑓(𝑥)=𝑓(𝑥)−^𝑚−1_𝑖=0𝑓^(𝑖)0⁺𝑥^𝑖,𝐷𝑖!_𝜈𝑥^𝛽=Γ𝛽+1Γ𝑥𝛽+1−𝜈^𝛽−𝜈,𝐷^𝜈𝜆𝑓(𝑥)+𝜇𝑔(𝑥)=𝜆𝐷^𝜈𝑓(𝑥)+𝜇𝐷^𝜈𝑔(𝑥). (7)

2.2. Shifted Legendre Polynomials

Assuming that the Legendre polynomial of degree 𝑗 is denoted by 𝑃_𝑗(𝑧) (defined on the interval (−1,1)), then 𝑃_𝑗(𝑧) may be generated by the recurrence formulae 𝑃_𝑗+1(𝑧)=2𝑗+1𝑗+1𝑧𝑃_𝑗(𝑧)−𝑗𝑃𝑗+1^𝑗−1(𝑧)𝑃,𝑗=1,2,…,₀(𝑧)=1,𝑃₁(𝑧)=𝑧. (8)

Considering 𝑧=2𝑥−1, Legendre polynomials are defined on the interval (0,1) that may be called shifted Legendre polynomials 𝑃^∗_𝑗(𝑥) that were generated using the following recurrence formulae: 𝑃^∗_𝑗+1(𝑥)=2𝑗+1()𝑃𝑗+12𝑥−1^∗_𝑗(𝑥)−𝑗𝑃𝑗+1^∗𝑗−1(𝑥),𝑃𝑗=1,2,…,_∗0(𝑥)=1,𝑃^∗₁(𝑥)=2𝑥−1. (9)

The orthogonality relation is ¹₀𝑃^∗_𝑖(𝑥)𝑃^∗_𝑗(𝑥)1𝑑𝑥=2𝑗+1,for𝑖=𝑗,0,for𝑖≠𝑗. (10)

The explicit analytical form of shifted Legendre polynomial 𝑃^∗_𝑗(𝑥) of degree 𝑗 may be written as 𝑃^∗_𝑗(𝑥)=^𝑗_𝑘=0(−1)^𝑗+𝑘𝑗+𝑘!𝑥^𝑘𝑗−𝑘!(𝑘!)², (11) and this in turn enables one to get 𝑃^∗_𝑗(0)=(−1)^𝑗,𝑃^∗_𝑗(1)=1. (12)

Any square integrable function 𝑢(𝑥) defined on the interval (0,1) may be expressed in terms of shifted Legendre polynomials as 𝑢(𝑥)=^∞_𝑗=0𝑎_𝑗𝑃^∗_𝑗(𝑥), (13) from which the coefficients 𝑎_𝑗 are given by 𝑎_𝑗=2𝑗+1¹₀𝑢(𝑥)𝑃^∗_𝑗(𝑥)𝑑𝑥,𝑗=0,1,…. (14)

If we approximate 𝑢(𝑥) by the first (𝑁+1)-terms, then we can write 𝑢_𝑁(𝑥)=^𝑁_𝑗=0𝑎_𝑗𝑃^∗_𝑗(𝑥), (15) which alternatively may be written in the matrix form 𝑢_𝑁(𝑥)≃𝐀^𝑇Ψ_𝑁(𝑥), (16) with 𝐀^𝑇≡𝑎₀,𝑎₁,…,𝑎_𝑁,Ψ_𝑁(𝑥)≡𝑃^∗₀(𝑥),𝑃^∗₁(𝑥),…,𝑃^∗_𝑁(𝑥)^𝑇. (17) Similarly, let 𝑢(𝑥,𝑡) be an infinitely differentiable function defined on 0<𝑥≤1 and 0<𝑡≤1. Then it may be expressed as 𝑢_𝑀,𝑁(𝑥,𝑡)≃^𝑀𝑁_𝑖=0𝑗=0𝑢_𝑖𝑗𝑃^∗_𝑖(𝑡)𝑃^∗_𝑗(𝑥)=Ψ^𝑇_𝑀(𝑡)𝐔Ψ_𝑁(𝑥), (18) with 𝑢𝐔=⁰⁰𝑢₀₁⋯𝑢_0𝑁𝑢₁₀𝑢₁₁⋯𝑢_1𝑁𝑢⋮⋮⋯⋮_𝑀0𝑢_𝑀1⋯𝑢_𝑀𝑁,𝑢_𝑖𝑗=(2𝑖+1)2𝑗+1¹₀𝑢(𝑥,𝑡)𝑃^∗_𝑖(𝑡)𝑃^∗_𝑗(𝑥)𝑑𝑥𝑑𝑡,𝑖=0,1,…,𝑀,𝑗=0,1,…,𝑁. (19)(20)

Theorem 3. The first derivative of the shifted Legendre vector Ψ_𝑁(𝑥) may be expressed as 𝑑Ψ_𝑁(𝑥)𝑑𝑥=𝐃Ψ_𝑁(𝑥), (21) where 𝐃 is the (𝑁+1)×(𝑁+1) operational matrix of derivative given by 𝐃=𝑑_𝑖𝑗=22𝑗+1,𝑓𝑜𝑟𝑗=𝑖−𝑘,𝑘=1,3,…,𝑁,𝑖𝑓𝑁𝑖𝑠𝑜𝑑𝑑,𝑘=1,3,…,𝑁−1,𝑖𝑓𝑁𝑖𝑠𝑒𝑣𝑒𝑛,0,𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒. (22)

Repeated use of (21) enables one to write 𝑑^𝑞Ψ_𝑁(𝑥)𝑑𝑥^𝑞=𝐃^𝑞Ψ_𝑁(𝑥), (23) where 𝑞 is a natural number and 𝐃^𝑞 means matrix power.

Theorem 4. The Riemann-Liouville fractional integral of order 𝜈 of the shifted Legendre polynomial vector Ψ_𝑀(𝑡) is given by 𝐼^𝜈Ψ_𝑀(𝑡)=𝐏_𝜈Ψ_𝑀(𝑡), (24) where 𝐏_𝜈 is the (𝑀+1)×(𝑀+1) operational matrix of fractional integral of order 𝜈 and is defined by 𝐏_𝜈=^0𝑘=0𝜉(0,0,𝑘)⁰_𝑘=0𝜉(0,1,𝑘)⋯⁰_𝑘=0𝜉(0,𝑀,𝑘_1𝑘=0)𝜉(1,0,𝑘)¹_𝑘=0𝜉(1,1,𝑘)⋯¹_𝑘=0𝜉(1,𝑀,𝑘)⋮⋮⋯⋮_𝑖𝑘=0𝜉(𝑖,0,𝑘)^𝑖_𝑘=0𝜉(𝑖,1,𝑘)⋯^𝑖_𝑘=0𝜉(𝑖,𝑀,𝑘)⋮⋮⋯⋮_𝑀𝑘=0𝜉(𝑀,0,𝑘)^𝑀_𝑘=0𝜉(𝑀,1,𝑘)⋯^𝑀_𝑘=0𝜉(),𝑀,𝑀,𝑘 (25) where 𝜉=𝑖,𝑗,𝑘2𝑗+1^𝑗_𝑙=0(−1)^{𝑖+𝑗+𝑘+𝑙}()!!×()()!×()𝑖+𝑘𝑙+𝑗𝑖−𝑘!𝑘!Γ𝑘+𝜈+1𝑗−𝑙𝑙!²(𝑘+𝑙+𝜈+1)⁻¹. (26)

(See [37] for proof.)

3. Legendre tau-Spectral Method

In this section, the Legendre operational matrix of fractional integrals is applied with the help of Legendre tau-spectral method to solve the T-FHE with the nonlocal condition.

Consider the T-FHE with the nonlocal condition 𝜕^𝜈𝑢(𝑥,𝑡)𝜕𝑡^𝜈−𝜕²𝑢(𝑥,𝑡)𝜕𝑥²(),𝑢()()(𝑥)𝑢()=𝑞𝑥,𝑡0<𝑥≤1,0<𝑡≤1,0<𝜈≤1,𝑥,0=𝑢𝑥,1+𝑓,0<𝑥≤1,0,𝑡=𝑔₀(𝑡),𝑢(1,𝑡)=𝑔₁(𝑡),0<𝑡≤1. (27) We integrate (27) of order 𝜈 and making use of (7), we have 𝑢(𝑥,𝑡)−𝑢(𝑥,1)+𝑓(𝑥)−𝐼^𝜈_𝑡𝜕²𝑢(𝑥,𝑡)𝜕𝑥²=𝐼^𝜈_𝑡𝑞()𝑢()𝑥,𝑡,0<𝜈≤1,0,𝑡=𝑔₀(𝑡),𝑢(1,𝑡)=𝑔₁(𝑡),0<𝑡≤1. (28)

In order to use tau-spectral method based on the shifted Legendre operational matrix for fractional integrals to solve the fully integrated problem (28), we approximate (𝑥,𝑡),𝑓(𝑥), and 𝑞(𝑥,𝑡) by the shifted Legendre polynomials as 𝑢_𝑀,𝑁(𝑥,𝑡)=Ψ^𝑇_𝑀(𝑡)𝐔Ψ_𝑁(𝑥),𝑓_𝑁(𝑥)=Ψ^𝑇_𝑀(𝑡)𝐅Ψ_𝑁(𝑥),𝑞_𝑀,𝑁(𝑥,𝑡)=Ψ^𝑇_𝑀(𝑡)𝐐Ψ_𝑁(𝑥), (29) where 𝐔 is the unknown coefficients (𝑀+1)×(𝑁+1) matrix and 𝐅 and 𝐐 are known matrices that can be written as 𝑓𝐅=⁰𝑓₁⋯𝑓_𝑁,𝑞00⋯0⋮⋮⋯⋮00⋯0𝐐=⁰⁰𝑞₀₁⋯𝑞_0𝑁𝑞₁₀𝑞₁₁⋯𝑞_1𝑁𝑞⋮⋮⋯⋮_𝑀0𝑞_𝑀1⋯𝑞_𝑀𝑁, (30) where 𝑓_𝑗 and 𝑞_𝑖𝑗 are given as in (14) and (20), respectively.

Using (29), it is easy to write 𝑢_𝑀,𝑁(𝑥,1)=Ψ^𝑇_𝑀(1)𝐔Ψ_𝑁(𝑥),=^𝑀_𝑘=0𝑢^𝑀_𝑘0𝑘=0𝑢_𝑘1⋯^𝑀_𝑘=0𝑢_𝑘𝑁Ψ_𝑁(𝑥),=Ψ^𝑇_𝑀(𝑡)𝐕Ψ_𝑁(𝑥), (31) where 𝐕 is a (𝑀+1)×(𝑁+1) matrix that can be written as 𝐕=^𝑀𝑘=0𝑢^𝑀_𝑘0𝑘=0𝑢_𝑘1⋯^𝑀_𝑘=0𝑢_𝑘𝑁.00⋯0⋮⋮⋯⋮00⋯0 (32) Making use of (23), (24), and (29) enables one to write 𝐼^𝜈_𝑡𝜕²𝑢(𝑥,𝑡)𝜕𝑥²≃𝐼^𝜈_𝑡Ψ^𝑇_𝑀(𝑡)𝐔𝐷²Ψ_𝑁(𝑥),=Ψ^𝑇_𝑀(𝑡)𝐏^𝑇_𝜈𝐔𝐃²Ψ_𝑁(𝑥). (33) In addition, if we use (24) and (29), we obtain 𝐼^𝜈_𝑡𝑞_𝑀,𝑁(𝑥,𝑡)=Ψ^𝑇_𝑀(𝑡)𝐏^𝑇_𝜈𝐐Ψ_𝑁(𝑥). (34) Equations (31) and (34) enable one to write the residual 𝑅_𝑀,𝑁(𝑥,𝑡) for (28) in the form 𝑅_𝑀,𝑁(𝑥,𝑡)=Ψ^𝑇_𝑀(𝑡)𝐔−𝐕−𝐅−𝐏^𝑇_𝜈𝐔𝐃²−𝐏^𝑇_𝜈𝐐Ψ_𝑁(𝑥). (35) As in a typical tau method (see [22, 38, 39]) we generate (𝑀+1)(𝑁−1) linear algebraic equations in the unknown expansion coefficients, 𝑢_𝑖𝑗, 𝑖=0,1,…, 𝑀;𝑗=0,1,…,𝑁−2, namely; ¹₀𝑅_𝑀,𝑁(𝑥,𝑡)𝑃^∗_𝑖(𝑡)𝑃^∗_𝑗(𝑥)𝑑𝑥𝑑𝑡=0,𝑖=0,1,…,𝑀,𝑗=0,1,…,𝑁−2, (36) and the rest of linear algebraic equations are obtained from the boundary conditions, as Ψ^𝑇_𝑀(𝑡)𝐔Ψ_𝑁(0)=𝑔₀𝑡_𝑖,Ψ^𝑇_𝑀(𝑡)𝐔Ψ_𝑁(1)=𝑔₁𝑡_𝑖,𝑖=0,1,…,𝑀, (37) where 𝑡_𝑖,𝑖=0,1,…,𝑀 are the roots of 𝑃^∗_𝑀+1(𝑡). The number of the unknown coefficients 𝑢_𝑖𝑗 is equal to (𝑀+1)(𝑁+1) and can be obtained from (36) and (37). Consequently 𝑢_𝑀,𝑁(𝑥,𝑡) given in (29) can be calculated.

4. Numerical Experiments

In order to highlight the accuracy of the presented scheme, we implement it to solve two numerical examples, and also comparisons between their exact solutions with the approximate solutions achieved using the presented scheme and with those achieved using other methods are made.

Example 1. We consider the following problem [36]: 𝜕^0.5𝑢(𝑥,𝑡)𝜕𝑡^0.5−𝜕²𝑢(𝑥,𝑡)𝜕𝑥²=2𝑡^1.5sin(2𝜋𝑥)Γ(1.5)+4𝜋²𝑡²(),𝑢()()()𝑢()()sin2𝜋𝑥𝑥,0=𝑢𝑥,1−sin2𝜋𝑥,0<𝑥≤1,0,𝑡=0,𝑢1,𝑡=0,0<𝑡≤1, (38) with exact solution 𝑢(𝑥,𝑡)=𝑡²sin(2𝜋𝑥).

Karatay et al. [36] introduced this problem and applied the IDA method to approximate its solution at various choices of time and space nodes 𝑀 and 𝑁.

We apply our numerical scheme for this problem. In order to show that our scheme is more accurate than the IDA method, in Table 1, we compare the maximum absolute errors (MAEs) achieved using our scheme with those obtained using the IDA [36] method at different values of 𝑁, (𝑁=𝑀). Moreover, Figure 1 plots the absolute error function at 𝑀=𝑁=18, while Figure 2 plots the absolute error function for 𝑡=0.5 at 𝑀=𝑁=18.

Table 1 

Comparison of our scheme with the IDA [36] at various choices of 𝑁, (𝑁=𝑀) for Example 1.

Figure 1 

Absolute error function at 𝑁=𝑀=18 for Example 1.

Figure 2 

Absolute error function at 𝑁=𝑀=18 with 𝑡=0.5 for Example 1.

Example 2. Consider the following problem: 𝜕^𝜈𝑢(𝑥,𝑡)𝜕𝑡^𝜈−𝜕²𝑢(𝑥,𝑡)𝜕𝑥²=2𝑡^2−𝜈ln(1+𝑥(1−𝑥))Γ()+𝑡3−𝜈₂2𝑥²−2𝑥−1𝑥²−𝑥+1²𝑢()()((𝑢()(),0<𝜈≤1,𝑥,0=𝑢𝑥,1−ln1+𝑥1−𝑥)),0<𝑥≤1,0,𝑡=0,𝑢1,𝑡=0,0<𝑡≤1, (39) with exact solution 𝑢(𝑥,𝑡)=𝑡²ln(1+𝑥(1−𝑥)).

Karatay et al. [36] introduced this problem and solved it for two choices of 𝜈, 𝜈=0.45,0.95 at different values of 𝑁 and 𝑀. Table 2 lists the MAEs for 𝜈=0.45,0.95 using our scheme at 𝑁=𝑀=4,6,8,10,12,14,16,18 and a comparison with those obtained in [36] at 𝑀=16,𝑁=2,4,8,16,32,64,128,256 and 𝑁=16, 𝑀=2,4,8,16,32,64,128,256. Figures 3 and 4 plot the absolute error functions at 𝑀=𝑁=20 with 𝜈₁=0.45 and 𝜈=0.95, respectively. Finally, in order to demonstrate the convergence of the proposed method, in Figure 5, we plot the logarithmic graphs of the maximum absolute errors (log₁₀Error) at two choices of 𝜈,𝜈=0.45,0.95 and various choices of 𝑁, (𝑁=𝑀), by using the presented algorithm.

Table 2 

Comparison of our scheme with the IDA [36] at various choices of 𝑁 and 𝑀 for Example 2.

Figure 3 

Absolute error function at 𝑁=𝑀=20 with 𝛾=0.45 for Example 2.

Figure 4 

Absolute error function at 𝑁=𝑀=20 with 𝛾=0.95 for Example 2.

Figure 5 

Absolute error function at 𝑁=𝑀=20 with 𝛾=0.95 for Example 2.

From Tables 1 and 2 and Figures 1 and 2 introduced above, it is shown that the proposed scheme is more accurate than the IDA method introduced by Karatay et al. [36].

5. Conclusion

An effective and accurate numerical scheme was developed to approximate the solution of the T-FHE with the nonlocal condition. The developed approach is based on the Legendre tau-spectral method combined with the operational matrix of fractional integration (described in the Riemann-Liouville sense) for orthogonal polynomials. A good approximation of the exact solution was achieved by using a limited number of the basis function.

The logarithmic graphs of the maximum absolute errors were presented to achieve the exponential convergence of the proposed method. Comparisons between our approximate solutions of test problems with their exact solutions and the approximate solutions achieved by the IDA method were introduced to confirm the validity and accuracy of our scheme.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

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.