Quadratic Finite Volume Element Schemes over Triangular Meshes for a Nonlinear Time-Fractional Rayleigh-Stokes Problem

Yanlong Zhang,Yanhui Zhou and Jiming Wu

1Graduate School of China Academy of Engineering Physics,Beijing,100088,China

2School of Data and Computer Science,Sun Yat-Sen University,Guangzhou,510275,China

3Institute of Applied Physics and Computational Mathematics,Beijing,100088,China

ABSTRACT In this article,we study a 2D nonlinear time-fractional Rayleigh-Stokes problem,which has an anomalous subdiffusion term,on triangular meshes by quadratic finite volume element schemes.Time-fractional derivative,defined by Caputo fractional derivative,is discretized through L2 ?1σ formula,and a two step scheme is used to approximate the time first-order derivative at time tn?α/2,where the nonlinear term is approximated by using a matching linearized difference scheme.A family of quadratic finite volume element schemes with two parameters are proposed for the spatial discretization,where the range of values for two parameters are β1 ∈(0,1/2),β2 ∈(0,2/3).For testing the precision of numerical algorithms,we calculate some numerical examples which have known exact solution or unknown exact solution by several kinds of quadratic finite volume element schemes,and contrast with the results of an existing quadratic finite element scheme by drawing diversified comparison plots and showing the detailed data of L2 error results and convergence orders.Numerical results indicate that,L2 error estimate of one scheme with parametersan d L2 error estimates of other schemes are,where h and Δt denote the spatial and temporal discretization parameters,respectively.

KEYWORDS Quadratic finite volume element schemes;anomalous sub-diffusion term; L2 error estimate;quadratic finite element scheme

1 Introduction

Recently,due to the widespread use of fractional partial differential equations(FPDEs),such as dispersion in a porous medium,statistical mechanics,mathematical biology and so on,numerical solution of FPDEs becomes one of the frontier fields in the research.Fractional partial differential equations can be roughly classified into three categories:Space FPDEs[1–10],time FPDEs[11–29]and space-time FPDEs[30–34].Anomalous sub-diffusion equations,one type of time FPDEs,arise in some physical and biological processes.And the study of FPDEs with anomalous sub-diffusion terms,such as modified anomalous sub-diffusion equations[17,18],fractional Cable equations[11,16]or others,is also meaningful and popular.The problem considered in this article,which belongs to a nonlinear time-fractional Rayleigh-Stokes problem[19–22]applied in some non-Newtonian fluids,is a variant of the Stokes’first problems and Rayleigh-Stokes problems[35–38],and it is important in physics and engineering.

At present,numerical simulation is an important and effective way to solve partial differential equations,and the relevant numerical methods can be finite difference methods[7,8,14–17],finite element methods(FEMs)[1,2,9,11–13,21,29,31–33,39],meshless methods[40,41],finite volume methods[3–6,42–54]and so on.Of course,the research for FPDEs by finite volume element methods(FVEMs)[10,23–28]has no exception for the local conservation and simple implementation.Sayevand et al.[23]presented a spatially semi-discrete piecewise linear FVEM for the time-fractional sub-diffusion problem and obtained some error estimates of the solution in both FEMs and FVEMs.A linear finite volume element scheme for the 2D time-fractional anomalous sub-diffusion equations was studied and analyzed by Karaa et al.[24],where the convergence rate was of orderh2+Δt1+αin theL∞(L2)norm and the results were improved in[25]for both smooth and nonsmooth initial data.Badr et al.[26]proposed a linear FVEM for the timefractional advection diffusion problem in one-dimension,and proved that the fully discrete scheme is unconditionally stable.Furthermore,Yazdani et al.[10]solved a space-fractional advectiondispersion problem in one-dimension by using linear FVEM and proved it is stable when the mesh grid size is small enough.Zhao et al.[27]constructed a mixed finite volume element scheme for the time-fractional reaction-diffusion equation,and showed the unconditional stability analysis for it.Moreover,Zhao et al.[28]proposed a linear FVEM for the nonlinear time-fractional mobile/immobile transport equations on triangular grids,and obtained the optimal priori error estimates inL∞(L2)andL2(H1)norms.To our knowledge,the study of high order finite volume element methods for 2D FPDEs is undiscovered.

There are some research about quadratic finite volume element methods for solving partial differential equations on triangular meshes.Tian et al.[42]presented quadratic element generalized differential methods to solve elliptic equations where two parameters of the quadratic element wereβ1=β2=1/3(referring to the definition Eqs.(4),(5)).Liebau[43]solved one type of elliptic boundary value problems by a quadratic element scheme with parametersβ1=1/4,β2=1/3,and provederror estimate under some assumption conditions.Xu et al.[44]started to study the structure about two parameters of the quadratic element,and improved some existing coercivity results.Chen et al.[45]established a general framework for construction and analysis of the higher-order finite volume methods.Wang et al.[46]established a unified framework to perform theL2error analysis for high order finite volume methods on triangular meshes,and proposed a new quadratic scheme with parametersto achieve the optimalL2convergence order.An unconditionally stable quadratic finite volume scheme with parametersfor elliptic equations was presented by Zou[47],and it had optimal convergence orders underH1norms.For the quadratic finite volume element schemes with parametersZhou et al.[48]obtained an analytic minimum angle condition and an optimalH1error estimates under the improved coercivity result.Moreover,a unified framework for the coercivity analysis of a class of quadratic schemes with parametersβ1∈(0,1/2),β2∈(0,2/3)was established for elliptic boundary value problems[49],which covered all the existing quadratic schemes of Lagrange type,and minimum angle conditions of the existing literatures are improved.All the above papers are mainly confined to the elliptic problems and we have seen some applications to other problems[50–53],but none relevant study for FPDEs.

In this article,the quadratic finite volume element method is proposed to solve one class of FPDEs,that is,a 2D nonlinear time-fractional Rayleigh-Stokes problem with the time-fractional derivative defined by Caputo fractional derivative.In spatial direction,this problem is solved by a class of quadratic finite volume element schemes with two parametersβ1andβ2.Moreover,this problem is discretized at timetn?α/2and the time-fractional derivative is discretized throughL2 ?1σformula in time direction.Numerical experiments indicate the efficiency of the schemes,specifically,theL2error estimate of one scheme isO(h3+Δt2),andL2error estimates of other schemes areWe find that the new finite volume element schemes are comparable with an existing finite element scheme[29].

The outline of this paper is as follows.In Section 2,we describe in details the specific algorithm steps of the quadratic finite volume element schemes over triangular meshes,and finally obtain the fully discrete schemes.In Section 3,some numerical experiments are performed to investigate the performance of the quadratic finite volume element schemes.The numerical results are also compared with those of an existing quadratic finite element scheme.A brief conclusion ends this article in last section.

2 Quadratic Finite Volume Element Schemes

2.1 Preliminary

In this article,we construct a family of quadratic finite volume element schemes to solve the following 2D nonlinear time-fractional Rayleigh-Stokes problem:

denotes the Caputo fractional derivative of orderα∈(0,1),which is an anomalous sub-diffusion term.f(u)is a nonlinear term subjected to the following conditions,|f(u)|≤C1|u|and |f′(u)|≤C2whereC1andC2are positive constants,andg(x,t)is the source term.

For the numerical solution of Eqs.(1)–(3),the space domainΩis first triangulated to get the so-called primary meshsee the solid line segments in Fig.1.The trial function spaceUhis then defined with respect togiven by

whereKdenotes a generic triangular element andis the set of all polynomials of degree less than or equal to 2.It is easy to see thatUhis the same one as that of the standard quadratic finite element method,but the test function space here is different.

Figure 1:The primary mesh and dual mesh

In order to define the test function space,we need to construct a dual mesh associated withsee the dashed line segments in Fig.1.For any triangular elementK=ΔP1P2P3withP1=(x1,y1),P2=(x2,y2)andP3=(x3,y3),we denote the barycenter ofKasQ,and the midpoints of the three sides asM1,M2andM3,respectively,see Fig.2.Gk,k+1andGk+1,kare the two points onPkPk+1,satisfying

wherekis a periodic index with Period 3.Using the above notations,Kcan be further partitioned into six subcells,i.e.,three quadrilaterals and three pentagons,see Fig.2.Letbe the set of all vertices and edge midpoints on the primary meshThen,the dual cell associated withis defined as the union of the subcells sharingP,and the dual mesh is defined asNow,the test function space is chosen as

Figure 2:Partition of the triangular element K

2.2 Semi-Discrete Schemes

In this part,we propose the following spatial discrete formulation of Eq.(1),

By the divergence theorem and the definition ofVh,we have,

where ?is the gradient operator andis the unit normal vector outward toThe left hand-side of Eq.(6)can be rewritten as

For the computation of the terms in Eq.(7),we introduce the following affine mapping that transformsKontoin(λ1,λ2)plane,

where

where

Figure 3:The reference element and its associated subcells

andλ3=1 ?λ1?λ2.Now we rewrite Eq.(7)to get the following semi-discrete schemes,

2.3 Fully Discrete Schemes

Next we introduce the fully discrete schemes at timetn?α/2.Let 0=t0

Lemma 2.1.([14],Lemma 2)Supposez(t)∈C3[0,T].Then,we have

Lemma 2.2.([11],Lemma 2)Assume thatf(t)∈C2[0,T].Then,the following second-order formula for the approximation of the nonlinear term at timetn?α/2holds,

Lemma 2.3.([15],Lemma 2)Supposez(t)∈C3[0,T].We have the followingL2 ?1σformula at timetn?α/2,

with

Based on Lemmas 2.1–2.3,we propose the following fully discrete schemes by the Eq.(9)at timetn?α/2:forn=1,

Let the basis functions of the trial function spaceUhbe denoted as?k(x,y),(k=1,2,...,m)wheremis the number of unknownsthen the numerical solutionMoreover,we knowwhereare transformed from?k(x,y),(k=1,2,...,m)by the affine mapping Eq.(8).By simplifying and synthesizing the Eqs.(13),(14),we obtain the following matrix form of the fully discrete schemes:forn=1,

forn≥2,

where

For the above finite volume element schemes,we emphasize that the nonlinear term is approximated by using the linearized difference scheme in Lemma 2.2,and none nonlinear iteration is involved.The whole algorithm is summarized below.

Algorithm 1:The quadratic finite volume element schemes(15)–(20)Step 1:Set n=0 and compute u0h=u0(x) by using the initial condition in(3);Step 2:Set n=1 and compute f0 by the expression f(u) and(17).Then,solve the linear system(15)to obtain u1h;Step 3:Do n=2 to N Compute fn?1 and fn?2 by using the expression f(u) and(17);Obtain unh by solving the linear system(16);End do

3 Numerical Examples

In this section,we use Eqs.(15)–(20)to solve four examples on uniform triangular mesh(Mesh I)and random triangular mesh(Mesh II),respectively,see Figs.4,5.The first level of Mesh II is constructed from Mesh I by the following random distortion of the interior vertices,

whereω∈(0,0.5)is the disturbance coefficient,his the spatial mesh size,ηxandηyare two random numbers located in[?1,1].The subsequent level is refined by the standard bisection procedure from its previous level.In our numerical examples,we choose the disturbance coefficientω=0.2 for Mesh II(see Fig.5).Recall that our quadratic finite volume element schemes have two parametersβ1andβ2.Here we just investigate the following specific schemes:

?Third scheme(QFVE-3):β1=1/4,β2=1/3;

?Fourth scheme(QFVE-4):β1=β2=1/3.;

Figure 4:Uniform triangular mesh(Mesh I)

Figure 5:Random triangular mesh(Mesh II)

We remark that the counterparts of the above schemes for elliptic problems have been studied in[42,43,46,47],respectively.TheL2errorsEuand convergence ordersare defined as

respectively,whereh1andh2are the spatial mesh sizes of two successive meshes,andΔt1andΔt2are the mesh sizes of two successive time levels.Moreover,the results of the quadratic finite element scheme(QFE)[29]are also employed for comparison.

3.1 Example 1

Solve Eqs.(1)–(3)with the nonlinear termf(u)=u2and the source term

Figure 6:Comparison of the numerical solutions and the exact solution for Example 1 on Mesh I.(a)A full profile;(b)A local enlarged profile

Table 1:Error results and spatial convergence orders with Δt=1/2000 on Mesh I in Example 1

Table 2:Error results and spatial convergence orders with Δt=1/2000 on Mesh II in Example 1

Figure 7: L2 errors of the numerical solution on Mesh I in Example 1.(a) α=0.1;(b) α=0.2;(c)α=0.5;(d) α=0.9

3.2 Example 2

where

Figure 8: L2 errors of the numerical solution on Mesh II in Example 1.(a) α=0.1;(b) α=0.2;(c) α=0.5;(d) α=0.9

The exact solution to this example is:

Analogously,we calculate theL2errors and spatial convergence orders for several kinds of quadratic finite volume element schemes and quadratic finite element scheme,see Tabs.5,6.The convergence behavior is similar to that forExample 1.Moreover,whenα=0.1,0.2,0.5,0.9 we draw the log–log plots for these five schemes on Mesh II,and find that the converge orders don’t change withα,see Fig.9.

3.3 Example 3

We take the space-time domain×[0,T]=[0,1]2×[0,1],the nonlinear termf(u)=u3?uand the source term

Table 3:Error results and temporal convergence orders with h=1/160 on Mesh I in Example 1

Table 4:Error results and temporal convergence orders with h=1/160 on Mesh II in Example 1

Table 5:Error results and spatial convergence orders with Δt=1/2000 on Mesh I in Example 2

Table 6:Error results and spatial convergence orders with Δt=1/2000 on Mesh II in Example 2

Figure 9: L2 errors of the numerical solution on Mesh II in Example 2.(a)QFE;(b)QFVE-1;(c)QFVE-2;(d)QFVE-3;(e)QFVE-4

Figure 10:Contour plots of |u?uh|on Mesh I in Example 3.(a)QFE;(b)QFVE-1;(c)QFVE-2;(d)QFVE-3;(e)QFVE-4

Table 7:Numerical results with Δt=1/2000 on Mesh I in Example 3

Table 8:Numerical results with Δt=1/2000 on Mesh II in Example 3

Table 9:Numerical results with Δt=1/100 on Mesh I in Example 4

where Mittag–Leffler functionE1,2?α(t)is defined by

In numerical calculation of this example,we useto approximate the Mittag–Leffler functionE1,2?α(t).The exact solution to the model isu(x,t)=etsin(2πx)sin(2πy).In Fig.10,we draw contour plots for the absolute value of error between exact solutionuand numerical solutionuh,i.e.,|u?uh|,withα=0.5,h=1/10,Δt=1/2000 at timeT=1 on Mesh I.It is obvious that accuracy of the first three schemes(QFE,QFVE-1,QFVE-2)is better than that of QFVE-3 and QFVE-4,see Tabs.7,8 for some detailed data.

3.4 Example 4

In the last example,we choose the nonlinear termf(u)=u3?uand the source termg(x,t)=0 with initial conditionu0(x)=x(1 ?x3)y(1 ?y3),where×[0,T]=[0,1]2×[0,1].Because of unknown exact solution,we take the numerical solution withh=1/160,Δt=1/100 as the ‘exact’solution when computing the errors.The results are given in Tab.9 where one can see that these schemes still work in this situation.

4 Conclusions

In this article,we study a nonlinear time-fractional Rayleigh-Stokes problem by using the quadratic finite volume element method combined with a specific time discretization.In temporal direction,we use a two step scheme to approximate the equation at timetn?α/2,whereL2?1σformula is used to approximate the time-fractional derivative.The fully discrete schemes of quadratic finite volume element are suggested and we find that only one of these schemes achieves the optimal convergence order inL2norm in space direction.We calculate some numerical examples by several kinds of quadratic finite volume element schemes and quadratic finite element scheme,spaceL2error orders of the QFE and QFVE-1 schemes reach 3.Meanwhile,numerical results of other three quadratic finite volume element schemes(QFVE-2,QFVE-3,QFVE-4)are nearly 2 and lower than the optimal order of QFVE-1.The future work includes the stability analysis and error estimates by following the related results on elliptic problems[46,48,54].

Acknowledgement:The authors would like to thank the editor and the anonymous reviewers for their valuable suggestions.

Funding Statement:This work was partially supported by the National Natural Science Foundation of China(No.11871009).

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.