# The Rational Spectral Method Combined with the Laplace Transform for Solving the Robin Time-Fractional Equation

Authors: Lufeng Yang
Journal: Advances in Mathematical Physics
 (2020)

---

## Abstract

In this paper, the rational spectral method combined with the Laplace transform is proposed for solving Robin time-fractional partial differential equations. First, a time-fractional partial differential equation is transformed into an ordinary differential equation with frequency domain components by the Laplace transform. Then, the spatial derivatives are discretized by the rational spectral method, the linear equation with the parameters is solved, and the approximation Ux,s is obtained. The approximate solution at any given time, which is the numerical inverse Laplace transform, is obtained by the modified Talbot algorithm. Numerical experiments are carried out to demonstrate the high accuracy and efficiency of our method.


---

## Body

## 1. Introduction

Because the fractional-order calculus operator has nonlocality, it is suitable for describing the material involving memory and heredity in real life, and many practical problems can be described by fractional-order differential equations [1, 2]. Several studies have been conducted on the construction of approximate solutions for various fractional partial differential equations. Developing analytical and numerical methods for solving fractional PDEs is a very important task [3]. To solve such problems, we need to introduce special functions to express the exact solutions of fractional differential equations, which can be very difficult. Therefore, attempts have been made to propose numerical methods that approximate the solutions of such equations [4–9].The purpose of this article is to investigate a rational spectral method combined with the Laplace transform, to find approximate solutions for certain classes of Robin time-fractional PDEs with parameters that have derivatives in the sense of Caputo fractional derivatives as follows:(1)0CDtαu=px∂2u∂x2+qx∂u∂x+fx,t,0<x<L,0<t≤T,subject to the following initial and Robin boundary conditions:
(2)ux,0=ϕx,u0,t+a∂xu0,t=v1t,uL,t+b∂xuL,t=v2t,where 0<α≤1, a, b∈R−0, and f, p, q are continuous real-valued functions. This method results in an accurate solution that is continuous in the temporal domain and is computationally efficient. The use of the Laplace transform circumvents the need for time marching in the temporal domain, which is computationally expensive.The rest of this article is organized as follows: In Section2, we introduce some necessary notation and preliminary lemmas. In Section 3, we describe the method of implementing the proposed method. Numerical results are discussed in Section 4, and some conclusions are drawn in Section 5.


## 2. Preliminaries

In this section, we provide some definitions and basic properties of fractional calculus [10] and the Laplace Transform that are required for the subsequent development.


### 2.1. Review of Fractional Calculus

Definition 1.
The Riemann-Liouville derivative of fractional orderα of a function ft is defined as follows:
(3)0RDtαft=1Γn−αdndtn∫0tfτt−τα−n+1dτ,where n−1<α≤n.Definition 2.
The Caputo derivative of fractional orderα is defined as follows:
(4)0CDtαft=1Γn−α∫0tfnτt−τα−n+1dτ,where n−1<α≤n.For these two derivatives, we have the following properties:Lemma 3.
Ifft∈Cn0,∞ and n−1<α≤n, then
(5)0CDtαft=0RDtαft−∑k=0n−1tkk!fk0.Clearly, from Lemma3, it follows that if f0=0, the Riemann-Liouville and Caputo fractional derivatives are equivalent for 0<α<1.Like integer-order differentiation, the Caputo fractional differential has the property of linearity as follows:(6)0CDtαk1⋅f1t+k2⋅f2t=k1⋅0CDtαf1t+k2⋅0CDtαf2t,where k1, k2 are constants.


### 2.2. Laplace Transform

Definition 4.
Suppose thatft is a real- or complex-valued function of the variable t>0 and s is a real or complex parameter. We define the Laplace transform of f as follows:
(7)Fs=Lft=∫0∞e−stftdt.For the fractional derivative, the Laplace transform has the following properties:Lemma 5.
Ifut∈Cn0,∞ and n−1<α≤n, then the Laplace transform of the Caputo derivative is given by the following:
(8)L0CDtαut=sαUs−∑k=0m−1sα−k−1uk0,where Us=Lut.Lemma 6.
[11] If ut∈Cn0,∞ and n−1<α≤n, then the Laplace transform of the Riemann-Liouville derivative given by the following:
(9)L0RDtαut=sαUs−∑k=0m−1sk0RDtα−k−1ut,where Us=Lut.


## 2.1. Review of Fractional Calculus

Definition 1.
The Riemann-Liouville derivative of fractional orderα of a function ft is defined as follows:
(3)0RDtαft=1Γn−αdndtn∫0tfτt−τα−n+1dτ,where n−1<α≤n.Definition 2.
The Caputo derivative of fractional orderα is defined as follows:
(4)0CDtαft=1Γn−α∫0tfnτt−τα−n+1dτ,where n−1<α≤n.For these two derivatives, we have the following properties:Lemma 3.
Ifft∈Cn0,∞ and n−1<α≤n, then
(5)0CDtαft=0RDtαft−∑k=0n−1tkk!fk0.Clearly, from Lemma3, it follows that if f0=0, the Riemann-Liouville and Caputo fractional derivatives are equivalent for 0<α<1.Like integer-order differentiation, the Caputo fractional differential has the property of linearity as follows:(6)0CDtαk1⋅f1t+k2⋅f2t=k1⋅0CDtαf1t+k2⋅0CDtαf2t,where k1, k2 are constants.


## 2.2. Laplace Transform

Definition 4.
Suppose thatft is a real- or complex-valued function of the variable t>0 and s is a real or complex parameter. We define the Laplace transform of f as follows:
(7)Fs=Lft=∫0∞e−stftdt.For the fractional derivative, the Laplace transform has the following properties:Lemma 5.
Ifut∈Cn0,∞ and n−1<α≤n, then the Laplace transform of the Caputo derivative is given by the following:
(8)L0CDtαut=sαUs−∑k=0m−1sα−k−1uk0,where Us=Lut.Lemma 6.
[11] If ut∈Cn0,∞ and n−1<α≤n, then the Laplace transform of the Riemann-Liouville derivative given by the following:
(9)L0RDtαut=sαUs−∑k=0m−1sk0RDtα−k−1ut,where Us=Lut.


## 3. The Coupling Scheme

Applying of the Laplace transform to models (1) and (2), we can obtain the following:
(10)sαUx,s−ϕx=pxUxxx,s+qxUxx,s+Fx,s,U0,s+aUx0,s=V1s,UL,s+bUxL,s=V2s,where Fx,s=Lfx,t,V1s=Lv1t,V2s=Lv2t.


### 3.1. Rational Spectral Method

To solve Robin boundary value problems of ordinary differential equation (10), the rational spectral method is chosen to obtain a highly accurate numerical solution in the transformed domain.The rational functionpNx in a barycentric form which interpolates the function ux at N+1 distinct points xkk=0N can be expressed as follows [12]:
(11)ux≈pNx=∑k=0Nωk/x−xkuxk∑k=0Nωk/x−xk,where ωkk=0N are nonzero numbers called barycentric weights. In particular, for Chebyshev-Gauss-Lobatto points xk=−coskπ/N, the barycentric weights are chosen as follows [13]:
(12)ω0=12,ωk=−1k,k=1,2,⋯,N−1,ωN=−1N2.The derivatives ofpN can serve to determine the nth order differentiation matrix Djknj,k=0N associated with pN represented by (10) at the point xj:
(13)pNnxj=∑k=0Ndndxnωk/x−xkuxk∑l=0Nωl/x−xlx=xj=∑k=0NDjknuxk,j=0,1,⋯,N.The advantage of the barycentric representation is the simplicity of the formula for the entries of first- and second-order differentiation matrices [14]
(14)Djk1=ωkωjxj−xk,j≠k,−∑i≠kDji1,j=k,Djk2=2Djk1Djj1−1xj−xk,j≠k,−∑i≠kDji2,j=k.By introducing the transformx=L/2y+1, x∈0,L, and y∈−1,1 and defining U^y=Ux=UL/2y+1, then U′x=2/LU∧′y and U″x=4/L2U∧″y, and we can rewrite (10) as follows:
(15)sαU^is−ϕ^i=p^i4L2∑j=0NDij2U^js+q^i2L∑j=0NDij1U^js+F^is,U^−1+2aL∑j=0ND0j1U^js=V1s,U^1+2bL∑j=0NDNj1U^js=V2s.Equation (15) is a linear system of equations dependent on the transformed parameter s. When the parameter s is defined as the contour of integration, the numerical inversion of the Laplace transform can be used to solve system (15) to obtain the approximate solution of each point on the contour of integration in the complex domain.


### 3.2. Numerical Inversion of the Laplace Transform

Definition 7.
The inversion of the Laplace transform ofFs is defined as follows:
(16)ft=L−1Fs=12πi∫ΓestFsds.The numerical inversion of the Laplace transform is performed by approximating the Bromwich integral. Using a conformal mapping, the contour of integration may be deformed to obtain an integral that may be approximated by standard quadrature techniques.Here, we consider the modified Talbot contour in the form [15]
(17)Γ:sθ=σ+λθcotαθ+γiθ,−π≤θ≤π,where σ, λ, γ, and α are constants to be specified by the user. Under these conformal maps, the Bromwich integral becomes as follows:
(18)ux,t=12πi∫ΓestUx,sds=12πi∫−ππesθtUx,sθs′θdθ,where s′θ=λcotαθ−αθcsc2αθ+γi, which can be approximated by the 2M-panel midpoint rule as follows:
(19)ut≈uMx,t=12Mi∑n=−MM−1esθnts′θnUx,sθn,where θn=2n+1π/2M.If contour (18) is symmetric with respect to the real axis, then half of the transform evaluations can be saved. That is, one needs to consider only quadrature nodes in the upper (or lower) half-plane. We can obtain the following:
(20)uMx,t=Im1M∑n=0M−1esntsn′Ux,sn.The following modified Talbot contour [16]
(21)Γ:sθ=Mt−0.2407+0.2378θcot0.7409θ+0.1349iθis used in this paper; it has a convergence rate of Oe−2.716M.


## 3.1. Rational Spectral Method

To solve Robin boundary value problems of ordinary differential equation (10), the rational spectral method is chosen to obtain a highly accurate numerical solution in the transformed domain.The rational functionpNx in a barycentric form which interpolates the function ux at N+1 distinct points xkk=0N can be expressed as follows [12]:
(11)ux≈pNx=∑k=0Nωk/x−xkuxk∑k=0Nωk/x−xk,where ωkk=0N are nonzero numbers called barycentric weights. In particular, for Chebyshev-Gauss-Lobatto points xk=−coskπ/N, the barycentric weights are chosen as follows [13]:
(12)ω0=12,ωk=−1k,k=1,2,⋯,N−1,ωN=−1N2.The derivatives ofpN can serve to determine the nth order differentiation matrix Djknj,k=0N associated with pN represented by (10) at the point xj:
(13)pNnxj=∑k=0Ndndxnωk/x−xkuxk∑l=0Nωl/x−xlx=xj=∑k=0NDjknuxk,j=0,1,⋯,N.The advantage of the barycentric representation is the simplicity of the formula for the entries of first- and second-order differentiation matrices [14]
(14)Djk1=ωkωjxj−xk,j≠k,−∑i≠kDji1,j=k,Djk2=2Djk1Djj1−1xj−xk,j≠k,−∑i≠kDji2,j=k.By introducing the transformx=L/2y+1, x∈0,L, and y∈−1,1 and defining U^y=Ux=UL/2y+1, then U′x=2/LU∧′y and U″x=4/L2U∧″y, and we can rewrite (10) as follows:
(15)sαU^is−ϕ^i=p^i4L2∑j=0NDij2U^js+q^i2L∑j=0NDij1U^js+F^is,U^−1+2aL∑j=0ND0j1U^js=V1s,U^1+2bL∑j=0NDNj1U^js=V2s.Equation (15) is a linear system of equations dependent on the transformed parameter s. When the parameter s is defined as the contour of integration, the numerical inversion of the Laplace transform can be used to solve system (15) to obtain the approximate solution of each point on the contour of integration in the complex domain.


## 3.2. Numerical Inversion of the Laplace Transform

Definition 7.
The inversion of the Laplace transform ofFs is defined as follows:
(16)ft=L−1Fs=12πi∫ΓestFsds.The numerical inversion of the Laplace transform is performed by approximating the Bromwich integral. Using a conformal mapping, the contour of integration may be deformed to obtain an integral that may be approximated by standard quadrature techniques.Here, we consider the modified Talbot contour in the form [15]
(17)Γ:sθ=σ+λθcotαθ+γiθ,−π≤θ≤π,where σ, λ, γ, and α are constants to be specified by the user. Under these conformal maps, the Bromwich integral becomes as follows:
(18)ux,t=12πi∫ΓestUx,sds=12πi∫−ππesθtUx,sθs′θdθ,where s′θ=λcotαθ−αθcsc2αθ+γi, which can be approximated by the 2M-panel midpoint rule as follows:
(19)ut≈uMx,t=12Mi∑n=−MM−1esθnts′θnUx,sθn,where θn=2n+1π/2M.If contour (18) is symmetric with respect to the real axis, then half of the transform evaluations can be saved. That is, one needs to consider only quadrature nodes in the upper (or lower) half-plane. We can obtain the following:
(20)uMx,t=Im1M∑n=0M−1esntsn′Ux,sn.The following modified Talbot contour [16]
(21)Γ:sθ=Mt−0.2407+0.2378θcot0.7409θ+0.1349iθis used in this paper; it has a convergence rate of Oe−2.716M.


## 4. Numerical Experiments

To demonstrate the accuracy and efficiency of the rational spectral method combined with the Laplace transform, we use this novel method to solve certain classes of Robin time-fractional PDEs with respect to the Caputo fractional derivative in this section.In our computations, all experiments are performed using MATLAB (version R2014a) on a personal computer with a 2.5 Hz central processing unit (Intel Core i5-2450M), 4.00 GB of memory, and Windows 7 operating system.Example 1.
Consider the following linear fractional advection-diffusion PDE on a finite domainΩ=0,1⊗0,1, where α∈0,1:
(22)0CDtαu=∂2u∂x2−∂u∂x+Γα+6120ext5,subject to the following initial and Robin boundary conditions:
(23)ux,0=0,u0,t−∂xu0,t=0,u1,t−∂xu1,t=0.
Here, the exact solution of this problem isux,t=extα+5. The rational spectral method combined with the Laplace transform and the RKA (reproducing kernel algorithm) in [17] is used to solve this problem for various values of x,t∈Ω when α∈0.25,0.5. The errors between the numerical solutions and the exact solutions are listed in Table 1. And a comparison of the obtained exact and numerical solutions by the rational spectral method combined with the Laplace transform is tabulated in Table 2.Table 1
Maximum pointwise errors for the new method and the RKA [17].
x
t
α=0.25
α=0.5
This paper
RKA
This paper
RKA
0.5
0.5
1.0456e-07
1.8464e-05
1.7119e-07
9.2059e-06
1.0
3.9790e-06
8.5399e-06
7.7470e-06
1.7199e-05
1.0
0.5
1.7239e-07
1.9596e-05
2.8224e-07
4.1319e-05
1.0
6.5603e-06
2.1198e-05
1.2773e-06
9.6374e-06Table 2
Numerical and exact solutions of Equation (22) and absolute errors when M=12 and N=16.
x
t
α=0.25
α=0.5
Numerical
Exact
Error
Numerical
Exact
Error
0.5
0.2
3.5282e-04
3.5282e-04
8.52e-10
2.3594e-04
2.3595e-04
1.11e-09
0.4
1.3426e-02
1.3426e-02
3.24e-08
1.0678e-02
1.0678e-02
5.02e-08
0.6
1.1283e-01
1.1283e-01
2.72e-07
9.9306e-02
9.9307e-02
4.67e-07
0.8
5.1094e-01
5.1094e-01
1.23e-06
4.8321e-01
4.8322e-01
2.27e-06Compared with the existing RKA withp=q=25 as a two-dimensional partition of the domain Ω, the accuracy of the new method with M=12 and N=16 has been significantly improved for the case with different α. Figure 1 shows the variation in the maximum absolute error between the numerical solutions and the exact solution against M for various values of N=8, 12, and 16 and α=0.5. It can be seen that the maximum error decays exponentially with the increase in M and remains stable when M reaches approximately 15. The maximum absolute errors for a fixed value of N=16 and α=0.5 and for various values of M=8, 12, and 16 are shown in Figure 2.Example 2.
Consider the following linear fractional Navier-Stokes equation on a finite domainΩ=0,1⊗0,1, when α∈0,1:
(24)0CDtαu=x2+1∂2u∂x2+2x∂u∂x+gx,t,subject to the following initial and Robin boundary conditions:
(25)ux,0=cosπx,u0,t+∂xu0,t=tα+3+1,u1,t−∂xu1,t=−tα+3+1,where
(26)gx,t=cosπxΓ4+α6t3+πtα+3+12xsinπx+πcosπxx2+1.
Here, the exact solution of this problem isux,t=cosπxtα+3+1.Figure 1
Log-linear plots of the maximum absolute error againstM with different N.Figure 2
The maximum absolute error with differentM and N=16.We compare the performance of the new method against that of the box-type difference scheme in [18] to compute the linear fractional Navier-Stokes equation. Table 3 lists the results of the rational spectral method combined with the Laplace transform using N=16 and M=12 and of the box-type difference scheme with τ=1/32 and h=1/3000 for the cases with α∈0.2,0.5,0.8. A comparison of the obtained exact and numerical solutions by the rational spectral method combined with the Laplace transform is tabulated in Table 4 for the cases with α∈0.2,0.8, where xj stand for the jth ChebyshevGauss-Lobatto point in [0,1].Table 3
Maximum pointwise errors for the new method and for the box-type difference scheme [18].
α=0.2
α=0.5
α=0.8
This paper
Box-type
This paper
Box-type
This paper
Box-type
1.8238e-08
1.3948e-04
4.3960e-08
7.2848e-04
1.0382e-07
1.2356e-03Table 4
Numerical and exact solutions of Equation (24) and absolute errors when M=12 and N=16.
t
x
α=0.2
α=0.8
Numerical
Exact
Error
Numerical
Exact
Error
0.5
x3
1.1009
1.1009
6.7798e-10
1.0641
1.0641
2.8309e-09
x6
0.8494
0.8494
5.2259e-10
0.8211
0.8211
2.1837e-09
x9
0.0000
0.0000
1.5615e-13
0.0000
0.0000
3.7472e-13
x12
-0.8494
-0.8494
5.2280e-10
-0.8211
-0.8211
2.1833e-09
x15
-1.1009
-1.1009
6.7764e-10
-1.0641
-1.0641
2.8297e-09Figures3 and 4 show the variation of the maximum absolute error between the numerical solutions and the exact solution against M for various values of N=8, 12, and 16 and α=0.2. From the above calculation results, we can reach a conclusion similar to that of Example 1.Figure 3
Log-linear plots of the maximum absolute error againstM with different N.Figure 4
The maximum absolute error with differentM and N=16.


## 5. Conclusion

In this paper, the rational spectral method combined with the Laplace transform has been proposed to solve various certain classes of Robin time-fractional PDEs with parameters that have derivatives in the sense of Caputo fractional derivatives. Numerical experiments show that our method is very accurate and reliable in time-fractional problems; moreover, the solution is valid at any point in time on the prescribed temporal domain. Compared with a traditional time stepping method (FDM), the new method avoids the restrictions of the time step by the Courant condition. It is particularly suitable for the numerical simulation of long-time evolution problems.Here, we restricted our approach to solve one-dimensional time-fractional problems. Actually, the theoretical and numerical frameworks presented in this paper are essential for extension to more complicated problems. In the future work, we expect to expand our method to solve multidimensional problems and nonlinear problems.



---

*Source: 9865682-2020-01-09.xml*