On an explicit finite difference method for fracti

arXiv:cs/0311011v1[cs.NA]10Nov2003OnanexplicitfinitedifferencemethodforfractionaldiffusionequationsS.B.Yuste∗L.AcedoDepartamentodeF´ısica,UniversidaddeExtremadura,E-06071Badajoz,SpainAbstractAnumericalmethodtosolvethefractionaldiffusionequation,whichcouldalsobeeasilyextendedtomanyotherfractionaldynamicsequations,isconsidered.Thesefractionalequationshavebeenproposedinordertodescribeanomaloustransportcharacterizedbynon-MarkoviankineticsandthebreakdownofFick’slaw.Inthispaperwecombinetheforwardtimecenteredspace(FTCS)method,wellknownforthenumericalintegrationofordinarydiffusionequations,withtheGr¨unwald-Letnikovdefinitionofthefractionalderivativeoperatortoobtainanexplicitfrac-tionalFTCSschemeforsolvingthefractionaldiffusionequation.Theresultingmethodisamenabletoastabilityanalysis`alavonNeumann.Weshowthattheanalyticalstabilityboundsareinexcellentagreementwithnumericaltests.Com-parisonbetweenexactanalyticalsolutionsandnumericalpredictionsaremade.Keywords:Fractionaldiffusionequation,vonNeumannstabilityanalysis,parabolicintegro-differentialequationsPACS:02.70.Bf,05.40.+j,02.50.-r∗Correspondingauthor.Emailaddress:[e-mail](S.B.Yuste).URL:(S.B.Yuste).PreprintsubmittedtoElsevierScience1February20081IntroductionFractionaldifferentialequationshavebeenahighlyspecializedandisolatedfieldofmathematicsformanyyears[1].However,inthelastdecadetherehasbeenincreasinginterestinthedescriptionofphysicalandchemicalpro-cessesbymeansofequationsinvolvingfractionalderivativesandintegrals.Thismathematicaltechniquehasabroadpotentialrangeofapplications[2]:relaxationinpolymersystems,dynamicsofproteinmoleculesandthediffu-sionofcontaminantsincomplexgeologicalformations[3,4,5]aresomeofthemostrecentlysuggested[6].Fractionalkineticequationshaveprovedparticularlyusefulinthecontextofanomalousslowdiffusion(subdiffusion)[7].Anomalousdiffusionischaracter-izedbyanasymptoticbehaviorofthemeansquaredisplacementoftheformDx2(t)E∼2KγΓ(1+γ)tγ,(1)whereγistheanomalousdiffusionexponent.Theprocessisusuallyreferredtoassubdiffusivewhen0γ1.Ordinary(orBrownian)diffusioncorre-spondstoγ=1withK1=D(thediffusioncoefficient).Fromacontinuous(macroscopic)pointofview,thediffusionprocessisdescribedbythediffu-sionequationut(x,t)=Duxx(x,t),whereu(x,t)representstheprobabilitydensityoffindingaparticleatxattimet,andwhereuηζ...isthepartialderiva-tivewithrespecttothevariablesη,ζ...Fromamicroscopicpointofview,thecontinuousdescriptionisknowntobeconnectedwithaMarkovprocessinwhichthemicroscopicparticles(randomwalkers)performstochasticjumpsoffinitemeanandfinitevariance.Intheseconditionsthecentrallimittheo-remholdsforthesumofthesejumpsandEinstein’slawforthemeansquare2displacementensues[Eq.(1)withγ=1].Ontheotherhand,ifanunderlyingnon-Markovianmicroscopicprocessisassumedinwhichrandomwalkersperformjumpsattimeschosenfromadistributionwithanalgebraiclong-timetailt−γ−1,thenthediffusionprocessisanomalous[7,8].InthesecircumstancesthecentrallimittheorembreaksdownandonemustapplythegeneralizedL´evy-Gnedenkostatistics[7,9]whichformthebasisofEq.(1).Itturnsoutthattheprobabilitydensityfunctionu(x,t)thatdescribestheseanomalousdiffusiveparticlesfollowsthefractionaldiffusionequation[7,10,11,12]:∂∂tu(x,t)=Kγ0D1−γt∂2∂x2u(x,t)(2)where0D1−γtisthefractionalderivativedefinedthroughtheRiemann-Liouvilleoperator(seeSec.2).Fractionalsubdiffusion-advectionequations,andfrac-tionalFokker-Planckequationshavealsobeenproposed[13,14,15,16]andevensubdiffusion-limitedreactionshavebeendiscussedwithinthisframework[17].Inthemathematicalliterature,theseequationsareusuallyreferredtoasparabolicintegro-differentialequationswithweaklysingularkernels[18].Thesecurrentapplicationsoffractionaldifferentialequationsandmanyothersthatmaywellbedevisedinthenearfuturemakeitimperativetosearchformethodsofsolution.Someexactanalyticalsolutionsforafewcases,althoughimportant,havebeenobtainedbymeansoftheMellintransform[11,12]andthemethodofimages[19].Thepowerfulmethodofseparationofvariablescanalsobeappliedtofractionalequationsinthesamewayasfortheusualdiffusionequations(anexampleisgiveninSec.4).Anotherroutetosolvingfractionalequationsisthroughtheintegrationoftheproductofthesolu-tionofthecorrespondingnon-fractionalequation(theBrowniancounterpart3obtainedbysettingγ→1)andaone-sidedL´evystabledensity[7,20,21].However,asalsofortheBrowniancase,theavailabilityofnumericalmethodsforsolving(2)wouldbemostdesirable,especiallyforthosecaseswherenoanalyticalsolutionisavailable.OnepossibilitywasdiscussedrecentlybyR.Gorenfloetal.[22,23,24]whopresentedaschemetobuilddiscretemodelsofrandomwalkssuitablefortheMonteCarlosimulationofrandomvariableswithaprobabilitydensitygovernedbyfractionaldiffusionequations.AnothermorestandardapproachistobuilddifferenceschemesofthetypeusedforsolvingVolterratypeintegro-differentialequations[18].Inthisline,someim-plicit(backwardEulerandCrank-Nicholson)methodshavebeenproposed[18,25,26,27,28,29,30].InthispaperweshallusetheforwardEulerdifferenceformulaforthetimederivative∂u/∂tinEq.(2)tobuildanexplicitmethodthatwewillcallthefractionalForward