2007 JCP Finite difference spectral approximations for the timefractional

Finite difference/spectral approximations for the time-fractional diffusion equation qYumin Lin,Chuanju Xu*School of Mathematical Sciences,Xiamen University,361005Xiamen,ChinaReceived 9October 2006;received in revised form 18January 2007;accepted 2February 2007Available online 14February 2007AbstractIn this paper,we consider the numerical resolution of a time-fractional diffusion equation,which is obtained from the standard diffusion equation by replacing the first-order time derivative with a fractional derivative (of order a ,with 06a 61).The main purpose of this work is to construct and analyze stable and high order scheme to efficiently solve the time-fractional diffusion equation.The proposed method is based on a finite difference scheme in time and Legendre spectral methods in space.Stability and convergence of the method are rigourously established.We prove that the full dis-cretization is unconditionally stable,and the numerical solution converges to the exact one with order O ðD t 2Àa þN Àm Þ,where D t ;N and m are the time step size,polynomial degree,and regularity of the exact solution respectively.Numerical experiments are carried out to support the theoretical claims.Ó2007Elsevier Inc.All rights reserved.AMS Classification:65M12;65M06;65M70;35S10Keywords:Fractional diffusion equation;Spectral approximation;Stability;Convergence1.IntroductionThe use of fractional partial differential equations (FPDEs)in mathematical models has become increas-ingly popular in recent years.Different models using FPDEs have been proposed [4,14,16],and there has been significant interest in developing numerical schemes for their solution.Roughly speaking,FPDEs can be classified into two principal kinds:space-fractional differential equation and time-fractional one.One of the simplest examples of the former is fractional order diffusion equations,which are generalizations of classical diffusion equations,treating super-diffusive flow processes.Much of0021-9991/$-see front matter Ó2007Elsevier Inc.All rights reserved.doi:10.1016/j.jcp.2007.02.001qThis work was partially supported by National NSF of China under Grant 10531080,973High Performance Scientific Computation Research Program 2005CB321703,and Program of 985Innovation Engineering on Information in Xiamen University.*Corresponding author.Tel.:+865922580713;fax:+865922580608.E-mail address:cjxu@ (C.Xu).Journal of Computational Physics 225(2007)1533–1552/locate/jcpthe work published to date has been concerned with this kind of FPDEs(see e.g.[1,4–7,13,15,20]for a non-exhaustive list of references).In this paper,we consider the time-fractional diffusion equation(TFDE),obtained from the standard dif-fusion equation by replacing thefirst-order time derivative with a fractional derivative of order a,with 06a61.Time-fractional diffusion or wave equations are derived by considering continuous time random walk problems,which are in general non-Markovian processes.The physical interpretation of the fractional derivative is that it represents a degree of memory in the diffusing material[9].These models have been inves-tigated in analytical and numerical frames by a number of authors[8,11,18,19,21].Some of these authors have tried to construct analytical solutions to problems of time-fractional differential equations.For example, Schneider and Wyss[18]and Wyss[21]considered the time-fractional diffusion-wave equations.The corre-sponding Green functions and their properties are obtained in terms of Fox functions.Gorenflo et al.[8,9] used the similarity method and the method of Laplace transform to obtain the scale-invariant solution of TFDE in terms of the wright function.A similar construction of the solution to the time-fractional advec-tion–dispersion equations or TFDE in whole-space and half-space has been given in[10,11]by using the Fou-rier–Laplace transforms.However,published papers on the numerical solution of the TFDE equations are very sparse.Liu et al.[12] use afirst-orderfinite difference scheme in both time and space directions for this equation,where some sta-bility conditions are derived.Herein,we examine a practicalfinite difference/Legendre spectral method to solve the initial-boundary value time-fractional diffusion problem on afinite domain.An approach based on the backward differentiation combined with spatial collocation method is used to obtain estimates of (2-a)-order convergence in time and exponential convergence in space.It is also shown that the time-stepping scheme is unconditionally stable for all a2½0;1 .A series of numerical examples is presented and compared with the exact solutions to support the theoretical claims.Let us emphasize that the computation of the numerical solution of time-fractional differential equations has been generally limited to simple cases(low spatial dimension or small time integration)due to the‘‘glo-bal dependence’’ually,by the definition of the fractional derivative,the solution at a time t k depends on the solutions at all previous time levels t<t k.The fact that all previous solutions have to be saved to compute the solution at the current time level would make the storage very expensive if low-order methods are employed for spatial discretization.Contrarily,use of the spectral method can relax this stor-age limit because,as compared to a low-order method,the spectral method needs fewer grid points to pro-duce highly accurate solution.This is one of main advantages of the spectral method for FPDES.Another difficult point in solving TFDE lies in the discretization of the time-fractional derivative.The fact that the time-fractional derivative uses Caputo integral makes the standard schemes and corresponding numerical analysis not applicable.One of our main goals here is to propose an approach and provide an error analysis for this problem.The outline of the paper is as follows:first,we provide in Section2an analytical solution of the TFDE in a bounded domain.Second,afinite difference scheme for temporal discretization of this problem is pro-posed in Section3,where the stability and convergence analysis is given.A detailed error analysis is car-ried out for the semi-discrete problem,showing that the temporal accuracy is of(2-a)-order.In Section4, we construct a Legendre spectral collocation method for the spatial discretization of the TFDE.Error esti-mates are provided for the full discrete problem.Finally numerical experiments are presented in Section5 which support the theoretical error estimates.Some concluding remarks are given in thefinal section.2.Analytical solution of the TFDE in a bounded domainIn this section,wefirst describe the problem of fractional differential equation studied in this paper,and present some analytical solutions which will be found helpful in the comprehension of the nature of such a problem.Let L>0;T>0;K¼ð0;LÞ,consider the time-fractional diffusion equation of the formo a uðx;tÞo t a Ào2uðx;tÞo x2¼fðx;tÞ;x2K;0<t6Tð2:1Þ1534Y.Lin,C.Xu/Journal of Computational Physics225(2007)1533–1552subject to the following initial and boundary conditions:uðx;0Þ¼gðxÞ;x2K;ð2:2Þuð0;tÞ¼uðL;tÞ¼0;06t6T;ð2:3Þwhere a is the order of the time-fractional derivative.Here,we consider the case06a61.o a uðx;tÞo t a in(2.1)isdefined as the Caputo fractional derivatives of order a[16],given byo a uðx;tÞo t a ¼1Cð1ÀaÞZ to uðx;sÞo sd sðtÀsÞa;0<a<1:ð2:4ÞWhen a¼1,Eq.(2.1)is the classical diffusion equation:o uðx;tÞo t Ào2uðx;tÞo x2¼fðx;tÞ;x2K;0<t6T;ð2:5Þwhile the case a¼0corresponds to the classical Helmholtz elliptic equation.In fact the time derivative of inte-ger order in(2.5)can be obtained by taking the limit a!1in(2.4).In the case0<a<1,the definition of the Caputo fractional derivatives uses the information of the stan-dard derivatives at all previous time levels(non-Markovian process).If f 0,then by applying thefinite sine and Laplace transforms to(2.1),the analytical solution for the problem(2.1)–(2.3)can be obtained[1]asuðx;tÞ¼2LX1n¼0E aðÀa2n2t aÞsinðanxÞZ LgðrÞsinðanrÞd r;ð2:6Þwhere a¼p andE aðzÞ¼X1m¼0z mCða mþ1Þis the Mittag–Leffler function.We list here some special Mittag–Leffler-type functions having explicit expressions as follows:E1ðÀzÞ¼eÀz;E2ðzÞ¼coshðffiffizpÞ;E12ðzÞ¼X1m¼0z mC mþ1ÀÁ¼e z2erfcðÀzÞ;where erfc(z)is the error function complement defined byerfcðzÞ¼1ffiffiffip pZ1zeÀt2d t:Although expressed in a compact form,we note that using(2.6)to compute the exact solution is gen-erally difficult(especially for large time)due to slow convergence of the series E aðzÞfor a2ð0;1Þand big z.The exact solutions for a number of a are plotted in Fig.1,showing the smoothness of the solu-tions.The regularity of these solutions can be analyzed by using some properties of the Mittag–Leffler function[2].In the analysis of the numerical method that follows,we will assume that problem(2.1)–(2.3)has a unique and sufficiently smooth solution.3.Discretization in time:afinite difference schemeIn order to simplify the notations and without lose of generality,we consider the case f 0in the scheme construction and its numerical analysis.First,we introduce afinite difference approximation to discretize the time-fractional derivative.Lett k:¼k D t,k¼0;1;...;K,where D t:¼TK is the time step.To motivate the construction of the scheme,weuse the following formulation:for all06k6KÀ1,Y.Lin,C.Xu/Journal of Computational Physics225(2007)1533–15521535o a u ðx ;t k þ1Þo t a ¼1C ð1Àa ÞX k j ¼0Z t j þ1t j o u ðx ;s Þo s d sðt k þ1Às Þa¼1C ð1Àa ÞX k j ¼0u ðx ;t j þ1ÞÀu ðx ;t j ÞD t Z t j þ1t jd s ðt k þ1Às Þa þr k þ1D t ;ð3:1Þwhere r k þ1D t is the truncation error.It can be verified that the truncation error takes the following form:r k þ1D t6c u 1C ð1Àa ÞX k j ¼0Z t j þ1t j t j þ1þt j À2s ðt k þ1Às Þa d s þO ðD t 2Þ"#;ð3:2Þwhere c u is a constant depending only on u .For the first term in RHS of (3.2),we have1C ð1Àa ÞX k j ¼0Z t j þ1t j t j þ1þt j À2sðt k þ1Às Þa d s¼À1C ð1Àa ÞX k j ¼011Àa ð2j þ1ÞD t 2Àa½ðk Àj Þ1Àa Àðk þ1Àj Þ1Àa þ1C ð1Àa ÞX k j ¼021Àa D t 2Àa ½ðj þ1Þðk Àj Þ1Àa Àj ðk þ1Àj Þ1Àaþ1C ð1Àa ÞX k j ¼02ð1Àa Þð2Àa ÞD t 2Àa ½ðk Àj Þ2Àa Àðk þ1Àj Þ2Àa ¼D t 2Àa C ð2Àa Þðk þ1Þ1Àa þ2ðk 1Àaþðk À1Þ1Àa þðk À2Þ1Àa þÁÁÁþ11Àa Þh i À2D t 2ÀaC ð3Àa Þðk þ1Þ2Àa ¼D t 2Àa C ð2Àa Þðk þ1Þ1Àa þ2ðk 1Àa þðk À1Þ1Àa þðk À2Þ1Àa þÁÁÁþ11Àa ÞÀ22Àa ðk þ1Þ2Àa!:LetS ðk Þ¼ðk þ1Þ1Àaþ2ðk 1Àa þðk À1Þ1Àaþðk À2Þ1ÀaþÁÁÁþ11Àa ÞÀ22Àaðk þ1Þ2Àa:Incidentally,we find that j S ðk Þj is bounded for all a 2½0;1 and all k P 1,as proven in the following lemma.1536Y.Lin,C.Xu /Journal of Computational Physics 225(2007)1533–1552Lemma3.1.For all a2½0;1 and all K P1,it holdsj SðKÞj6c;where c is a constant independent of a;K.Proof.First,for a¼0,a direct calculation shows SðKÞ¼0for all K P1.Now we prove the lemma for a2ð0;1 .It can be verified thatSðKÞ¼ðKþ1Þ1Àaþ2ðK1ÀaþðKÀ1Þ1ÀaþðKÀ2Þ1ÀaþÁÁÁþ11ÀaÞÀ22ÀaðKþ1Þ2Àa¼X Kk¼0a k;wherea k¼ðkþ1Þ1Àaþk1ÀaÀ22Àaððkþ1Þ2ÀaÀk2ÀaÞ:This observation leads us to prove that the series P1k¼0a k converges.It is well known that the seriesP1k¼11k bconverges for all b>1.By consequence,it suffices to prove j a k j61k1þa for big enough k.In fact we have,for k P2,j a k j¼k1Àa1þ1 k1Àaþ1À2k2Àa1þ1k2ÀaÀ1!¼k1Àa1þ1þð1ÀaÞ1kþð1ÀaÞðÀaÞ2!1k2þð1ÀaÞðÀaÞðÀaÀ1Þ3!1k3þÁÁÁÀ2k2ÀaÀ1þ1þð2ÀaÞ1kþð2ÀaÞð1ÀaÞ2!1kþð2ÀaÞð1ÀaÞðÀaÞ3!1kþð2ÀaÞð1ÀaÞðÀaÞðÀaÀ1Þ4!1kþÁÁÁ¼k1Àa12!À23!ð1ÀaÞðÀaÞ1kþ13!À24!ð1ÀaÞðÀaÞðÀaÀ1Þ1kþÁÁÁ6k1Àa 13!ð1ÀaÞa1k21þ2ðaþ1Þ41kþ3ðaþ1Þðaþ2Þ201k2þÁÁÁ613!ð1ÀaÞa1k1þa1þ1kþ1k2þÁÁÁ623!ð1ÀaÞa1k1þa61k1þa:The proof is completed.hMore precise bound of SðkÞcan be obtained by numerical computations.Indeed our numerical tests show thatÀ16SðkÞ60806a61;k¼1;2;...,as observed in Fig.2,where evolution of SðkÞas a function of k for several typical values of a is plotted in the leftfigure,while limit of SðkÞas k tends to infinity as a function of a is shown in the rightfigure.Y.Lin,C.Xu/Journal of Computational Physics225(2007)1533–15521537Using Lemma3.1,and taking into account the fact that1Cð2ÀaÞ62for all a2½0;1 ,we have1Cð1ÀaÞX kj¼0Z t jþ1t jt jþ1þt jÀ2sðt kþ1ÀsÞad s62D t2Àa:As a result,it holdsr kþ1 D t K cuD t2Àa:ð3:3ÞOn the other side,a straightforward calculation of thefirst term in RHS of(3.1)gives1Cð1ÀaÞX kj¼0uðx;t jþ1ÞÀuðx;t jÞD tZ t jþ1t jd sðt kþ1ÀsÞa¼1Cð1ÀaÞX kj¼0uðx;t jþ1ÞÀuðx;t jÞD tZ t kþ1Àjt kÀjd tt a¼1Cð1ÀaÞX kj¼0uðx;t kþ1ÀjÞÀuðx;t kÀjÞD tZ t jþ1t jd tt¼1Cð2ÀaÞX kj¼0uðx;t kþ1ÀjÞÀuðx;t kÀjÞD t aðjþ1Þ1ÀaÀj1Àah i:For the sake of simplification,let us introduce the notations b j:¼ðjþ1Þ1ÀaÀj1Àa;j¼0;1;...;k,and define the discrete fractional differential operator L atbyL a t uðx;t kþ1Þ:¼1Cð2ÀaÞX kj¼0b juðx;t kþ1ÀjÞÀuðx;t kÀjÞD t a:Then(3.1)readso a uðx;t kþ1Þo t a ¼L atuðx;t kþ1Þþr kþ1D t:ð3:4ÞUsing L at uðx;t kþ1Þas an approximation of o a uðx;t kþ1Þaleads to the followingfinite difference scheme to(2.1):L a t u kþ1ðxÞ¼o2u kþ1ðxÞo x2;k¼0;1;...;KÀ1;ð3:5Þwhere u kþ1ðxÞis an approximation to uðx;t kþ1Þ.By virtue of(3.3),this scheme is formally of(2-a)-order accu-racy.A rigorous analysis of the convergence rate will be provided later for both semi-discrete and full-discrete cases,where we will prove that the temporal accuracy of the scheme(3.5)is globally of order2-a.Scheme(3.5)can be rewritten into,with simplification by omitting the dependence of u kþ1ðxÞon x:b0u kþ1ÀCð2ÀaÞD t a o2u kþ1o x2¼b0u kÀX kj¼1b j½u kþ1ÀjÀu kÀj ¼b0u kÀX kÀ1j¼0b jþ1u kÀjþX kj¼1b j u kÀj:ð3:6ÞIt is worthwhile to noting that the second term in the RHS of(3.6)automatically vanishes when k61,while the last two terms of(3.6)vanish when k¼0.It is direct to check thatb j>0;j¼0;1;...;k;1¼b0>b1>ÁÁÁ>b k;b k!0as k!1;X k j¼0ðb jÀb jþ1Þþb kþ1¼ð1Àb1ÞþX kÀ1j¼1ðb jÀb jþ1Þþb k¼1:ð3:7ÞLet us introduce the parameter a0:a0:¼Cð2ÀaÞD t a;and note that b0¼1,then by reformulating the right-hand side of(3.6),we obtain an equivalent form to scheme(3.5):1538Y.Lin,C.Xu/Journal of Computational Physics225(2007)1533–1552u kþ1Àa0o2u kþ1o x2¼ð1Àb1Þu kþX kÀ1j¼1ðb jÀb jþ1Þu kÀjþb k u0;k P1:ð3:8ÞHere again,when k¼1scheme(3.8)becomesu2Àa0o2u2o x2¼ð1Àb1Þu1þb1u0:For the special case k¼0,that is thefirst time step,the scheme simply readsu1Àa0o2u1o x2¼u0:ð3:9ÞEqs.(3.8)and(3.9),together with the boundary conditionsu kþ1ð0Þ¼u kþ1ðLÞ¼0;k P0;ð3:10Þand the initial conditionu0ðxÞ¼gðxÞ;x2Kð3:11Þform a complete set of the semi-discrete problem.It will be useful to define the error term r kþ1byr kþ1:¼a0o a uðx;t kþ1Þo t aÀL atuðx;t kþ1Þ!:ð3:12ÞThen we have from(3.3)and(3.4)j r kþ1j¼Cð2ÀaÞD t a j r kþ1D tj6c u D t2:ð3:13ÞTo introduce the variational formulation of the problem(3.8),we define some functional spaces endowed with standard norms and inner products that will be used hereafter.H1ðKÞ:¼v2L2ðKÞ;d vd x2L2ðKÞ&';H10ðKÞ:¼v2H1ðKÞ;v jo K¼0ÈÉ;H mðKÞ:¼v2L2ðKÞ;d k vk2L2ðKÞfor all positive integer k6m&';where L2ðKÞis the space of measurable functions whose square is Lebesgue integrable in K.The inner products of L2ðKÞand H1ðKÞare defined,respectively,byðu;vÞ¼ZK uv d x;ðu;vÞ1¼ðu;vÞþd ud x;d vd xand the corresponding norms byk v k0¼ðv;vÞ1=2;k v k1¼ðv;vÞ1=21:The norm kÁkmof the space H mðKÞis defined byk v km ¼X mk¼0d k vd x k2!1=2:In this paper,instead of using the above standard H1-norm,we prefer to define kÁk1byk v k1¼k v k2þa0d ud x2!1=2:ð3:14ÞY.Lin,C.Xu/Journal of Computational Physics225(2007)1533–15521539It is well known that the standard H1-norm and the norm defined by(3.14)are equivalent,the latter will be used in what follows.The variational(weak)formulation of the Eq.(3.8)subject to the boundary condition(3.10)reads:find u kþ12H1ðKÞ,such thatðu kþ1;vÞþa0o u kþ1o x;o vo x¼ð1Àb1Þðu k;vÞþX kÀ1j¼1ðb jÀb jþ1Þðu kÀj;vÞþb kðu0;vÞ;8v2H1ðKÞ:ð3:15ÞFor this weak semi-discretized problem,we have the following stability result.Theorem3.1.The semi-discretized problem(3.15)is unconditionally stable in the sense that for all D t>0,it holdsk u kþ1k16k u0k;k¼0;1;...;KÀ1:Proof.We will prove the result by induction.First when k¼0,we haveðu1;vÞþa0o u1o x;o vo x¼ðu0;vÞ8v2H1ðKÞ:Taking v¼u1and using the inequality k v k06k v k1and Schwarz inequality,we obtain immediatelyk u1k16k u0k:Suppose now we have provenk u j k16k u0k;j¼1;2;...;k;ð3:16Þwe want to prove k u kþ1k16k u0k.Taking v¼u kþ1in(3.15)givesðu kþ1;u kþ1Þþa0o u kþ1;o u kþ1¼ð1Àb1Þðu k;u kþ1ÞþX kÀ1j¼1ðb jÀb jþ1Þðu kÀj;u kþ1Þþb kðu0;u kþ1Þ:Hence,by using(3.16),we havek u kþ1k216ð1Àb1Þk u k kk u kþ1kþX kÀ1j¼1ðb jÀb jþ1Þk u kÀj kk u kþ1kþb k k u0kk u kþ1k6½ð1Àb1ÞþX kÀ1j¼1ðb jÀb jþ1Þþb k k u0kk u kþ1k1:Finally,the last equality of(3.7)yieldsk u kþ1k16k u0k:ÃNow we carry out an error analysis for the solution of the semi-discretized problem.We denote from now by c a generic constant which may not be the same at different occurrences.Theorem3.2.Let u be the exact solution of(2.1)–(2.3),f u k g K k¼0be the time-discrete solution of(3.15)with the initial condition u0ðxÞ¼uðx;0Þ,then we have the following error estimates:(1)when06a<1,k uðt kÞÀu k k16cu;aT a D t2Àa;k¼1;2;...;K;ð3:17Þwhere c u;a:¼c u=ð1ÀaÞ,with c u:constant defined in(3.3);(2)when a!1,k uðt kÞÀu k k16cuT D t;k¼1;2;...;K:ð3:18Þ1540Y.Lin,C.Xu/Journal of Computational Physics225(2007)1533–1552Proof.(1)First consider the case06a<1.We start by proving the following estimate:k uðt jÞÀu j k16cubÀ1jÀ1D t2;j¼1;2;...;K:ð3:19ÞAs in the proof of Theorem3.1,we will use the mathematical induction.Let e k¼u x;t kðÞÀu kðxÞ.For j¼1,we have,by combining(2.1),(3.9)and(3.12),the error equation:ð e1;vÞþa0o e1o x;o vo x¼ð e0;vÞþðr1;vÞ¼ðr1;vÞ8v2H1ðKÞ:Taking v¼ e1yieldsk e1k216k r1kk e1k:This,together with(3.13),givesk uðt1ÞÀu1k16cubÀ1D t2:Therefore,(3.19)is proven for the case j¼1.Suppose now(3.19)holds for all j¼1;2;...;k,we need then to prove that it holds also for j¼kþ1.By combining(2.1),(3.12)and(3.15),we derive,8v2H1ðKÞð e kþ1;vÞþa0o e kþ1o x;o vo x¼ð1Àb1Þð e k;vÞþX kÀ1j¼1ðb jÀb jþ1Þð e kÀj;vÞþb kð e0;vÞþðr kþ1;vÞ:ð3:20ÞLet v¼ e kþ1in(3.20),thenk e kþ1k216ð1Àb1Þk e k kk e kþ1kþX kÀ1j¼1ðb jÀb jþ1Þk e kÀj kk e kþ1kþb k k e0kk e kþ1kþk r kþ1kk e kþ1k:Dividing by k e kþ1k1at both sides,using the induction assumption and the fact that bÀ1jbÀ1jþ1<1for all non-negativeinteger j,we obtaink e kþ1k16ð1Àb1ÞbÀ1kÀ1þX kÀ1j¼1ðb jÀb jþ1ÞbÀ1kÀjÀ1"#c u D t2þc u D t26ð1Àb1ÞþX kÀ1j¼1ðb jÀb jþ1Þþb k"#c u bÀ1kD t2:Using(3.7)in the above inequality givesk e kþ1k16cubÀ1kD t2:The estimate(3.19)is proved.Now,by the definition of b k,a direct computation shows thatkÀa bÀ1kÀ1¼kÀak1ÀaÀðkÀ1Þ!11Àa;as k!1;and the function/ðxÞ:¼xÀax1ÀaÀðxÀ1Þ1Àais increasing on x for all x>1since/0ðxÞ¼1xðxÀ1Þa1À1À1xaÀax!P08x>1;06a61:This means kÀa bÀ1kÀ1increasingly tends to11Àaas1<k!1.Note that kÀa bÀ1kÀ1¼1when k¼1,hence we havekÀa bÀ1kÀ1611Àa;k¼1;2;...;K:ð3:21ÞConsequently we obtain,for all k such that k D t6T,Y.Lin,C.Xu/Journal of Computational Physics225(2007)1533–15521541k uðt kÞÀu k k16cubÀ1kÀ1D t2¼c u kÀa bÀ1kÀ1k a D t26c u11Àaðk D tÞa D t2Àa6c u;a T a D t2Àa:(2)Now we consider the case a!1.Note that in this case,the estimate(3.17)has no meaning since c u;a tends to infinity as a!1.Therefore,we need to look for an estimate of other form.Taking into account the fact j D t6T for all j¼1;2;...;K,we are led to establish:k uðt jÞÀu j k16cuj D t2;j¼1;2;...;K:ð3:22ÞOnce again,we derive this estimate by induction.(3.22)is obvious for j¼1.Suppose now(3.22)holds for j¼1;2;...;k,we want to prove that it remains true for j¼kþ1.A similar procedure as in case(1)leads tok e kþ1k16ð1Àb1Þk e k kþX kÀ1j¼1ðb jÀb jþ1Þk e kÀj kþb k k e0kþk r kþ1k6ð1Àb1Þ½c u k D t2 þX kÀ1j¼1ðb jÀb jþ1Þ½c uðkÀjÞD t2 þc u D t2¼ð1Àb1Þkkþ1þX kÀ1j¼1ðb jÀb jþ1ÞkÀjkþ1þ1kþ1"#c uðkþ1ÞD t2¼ð1Àb1ÞþX kÀ1j¼1ðb jÀb jþ1ÞÀð1Àb1Þ1kþ1ÀX kÀ1j¼1ðb jÀb jþ1Þjþ1kþ1þ1kþ1"#c uðkþ1ÞD t2:ð3:23ÞIt is observed thatð1Àb1Þ1kþ1þX kÀ1j¼1ðb jÀb jþ1Þjþ1kþ1þb k P1kþ1ð1Àb1ÞþX kÀ1j¼1ðb jÀb jþ1Þþb k"#¼1kþ1;which is equivalent toÀð1Àb1Þ1kþ1ÀX kÀ1j¼1ðb jÀb jþ1Þjþ1kþ1þ1kþ16bk:Bringing this inequality into(3.23)givesk e kþ1k16ð1Àb1ÞþX kÀ1j¼1ðb jÀb jþ1Þþb k"#c uðkþ1ÞD t2¼c uðkþ1ÞD t2:Hence(3.22)is proven,i.e.(3.18)holds.The proof is completed.h4.Full discretization4.1.A Galerkin spectral method in spaceTo simplify the notations,we let K¼ðÀ1;1Þhereafter.The Galerkin spectral discretization proceeds by approximating the solution by the polynomials of high degree.To this end,we define P NðKÞthe space ofall polynomials of degree6N with respect to x.Then the discrete space,denoted by P0N ðKÞ,is defined as fol-lows:P0N ðKÞ:¼H1ðKÞ\P NðKÞ.Now consider the spectral discretization to the problem(3.15)as follows:find u kþ1N 2P0NðKÞ,such that forall v N2P0NðKÞðu kþ1N ;v NÞþa0oo xu kþ1N;oo xv N¼ð1Àb1Þðu kN;v NÞþX kÀ1j¼1ðb jÀb jþ1Þðu kÀj N;v NÞþb kðu0N;v NÞ:ð4:1ÞFor f u jN g kj¼0given,the well-posedness of the problem(4.1)is guaranteed by the well-known Lax–MilgramLemma.We are interested in this section in deriving an error estimate for the full-discrete solution f u kN g K k¼0.1542Y.Lin,C.Xu/Journal of Computational Physics225(2007)1533–1552。