有限体积法求解分数阶cable方程

有限体积法求解分数阶cable方程 有限体积法求解分数阶cable方程 第29卷第2期 2O11年3月 泉州师范学院 JournalofQuanzhouNormalUniversity Vo1.29No.2 Mar.2011 有限体积法求解分数阶cable方程 沈淑君,吴春红,刘青霞. (1.华侨大学科学学院,福建泉州362021;2.厦t-J~A-r学院数理系,福建厦门361024; 3.厦门大学数学科学学院,福建厦门361005) 摘要:cable方程是神经元动力学中最重要的基本方程之一,而用于描述神经纤维活动的分数阶cable方 程能够更好地模拟神经元的动力学行为.文章采用有限体积法离散得到数值逼近格式,求解一维和二维的分数 阶cable方程.并用提出的数值方法求解一维和二维情况的两个数值例子,从而说明数值方法的有效性. 关键词:分数阶cable方程;数值近似;有限体积法 中图分类号:O241.82文献标识码:A文章编号:1009-8224(2011)02-0031-06 2O世纪4O年代,研究人员了解到神经轴突有着被动膜的特性,并且神经轴突可以被分布式RC电路所刻 画.这样,神经轴突的被动膜性质可以用潜艇电缆来模拟.1855年,WilliamThomson 通过潜艇电缆发现信号 的数学理论.而cable方程 -:.,,?, 山一一?u(x,,)一u(x,),?n,t>O,(1)0' 通常用于描绘神经元的电缆性质L1制. cable方程是从联系电动性质与细胞膜之间的现象学中产生的[1-23,并且也可以从物理学中描述离子电扩 散运动的Nemst-Planck方程导出[3.Nernst-Planck方程的模型已被用于模拟电压关联的浓度梯度.Qian和 Sejnowskics~应用一维的Nemst-Planck方程模拟在树突状脊柱中的激感性后联的电位.换言之,在一些简单的 假设下,神经元中离子电扩散运动的Nemst-Planck方程与cable方程是等价的.VanEgeraat和Wikswo等 人[6]应用一个一维Nemst-Planck方程研究长时间范围的受伤轴突的轴突传输.Lindsay等人[g]回顾了从 Maxwell'S方程到cable方程的发展过程.由此,cable方程是神经元动力学中的最重要的基本方程之一. 近年来,研究者已把传统的cable方程推广到用于描述神经纤维活动的带有分数阶导数的cable方程 — ./一+, 一D一?u(x,,)一一U(,),z?0,,>o,(2) C' 并且用这种带有分数阶导数的cable方程来模拟神经元的动力学行为效果更好. Henry和Langlands[1~3讨论了带有时间分数阶导数算子的cable方程,它可用于模拟具有电紧张性质的 刺神经元末端神经的活动.这些方程是由用来模拟树突状脊柱中诱捕性质引起的反常次扩散的带有分数阶算 子的Nernst—Planck方程中导出的.Langlands等人[112引入了分数阶Nernst-Planck方程,并指出分数阶cable 方程用来作为在神经细胞中当分子扩散是反常次扩散时离子的电扩散的宏观模型,并且得到在电流注入的无 穷或半无穷情况下带有特殊函数的级数形式的解. cable方程被用来模拟神经元的性质之后,研究者开始讨论cable方程的求解问题.姜常珍等[切建立阈下 线性状态外加刺激时的连续型电缆方程,提出利用积分变换方法求出点电源刺激下电缆方程的解析解,并用 计算机对无髓神经纤维细胞跨膜电位进行模拟仿真,动作电位产生的可能性.王鹤和姜常珍m利用傅立 收稿日期:2OlO一12—2o 作者简介:沈淑君(1977,),女,福建诏安人,讲师,博士,从事分数阶微分方程数值解研究. 基金项目:国家自然科学基金资助项目(11001090,11026094);福建省自然科学基金资助项目(2010J01011, 2010J05009);华侨大学科研启动基金资助项目(O8BS507) 32泉州师范学院2011年3月 叶变换方法求出了阈下电刺激时连续型电缆方程的解析式,并探讨了在工程上计算膜电压以及估算误差的方 法.Langlands等人[1t]利用分离变量法求得了分数阶cable方程的解析解,但得到解析解的方法具有较大的局 限性,针对不同的边界条件有可能得不到结果. 由于大多数分数阶微分方程的解析解除了简单的初边值条件外难以用显式给出,所以研究者更注重寻求 数值方法.Liu等人Ds-ls研究空间分数阶偏微分方程,将此方程转化为常微分方程系统,并用向后差商公式来 求解.Roop["]在中的一个有界区域内用连续的分片多项式做为三角剖分的基函数,研究Galerkin逼 近.Lin和xu[用有限差分和谱方法求解了时间分数阶扩散方程.Liu等人用随机游走和有限差分法讨论 了IAvy-Feller对流扩散过程.Zhuang等人EZo]在解反常次扩散方程时,提出了一个新的隐式数值方法和两个 改进收敛阶的数值技巧.本文利用有限体积法求解分数阶cable方程. 1数值逼近格式 1.1分数阶cable方程 本文考虑带有如下初边值条件: u(x,0)一?(z),?0;(3) u(x,,)l舶一(,),t>0(4) 的分数阶cable方程(2),其中:0?y1,y2<1,?是拉普拉斯算子,u(x,)示细胞膜电 位,0是一有界区域, 是的边界,而符号一和D分别表示1一y1,1一y2阶的Riemman-Liouville分数阶 导数,其定义为: ,, D],-Pu(x'z)=南蓦J坼 其中:I1(?)是Gamma函数. 下面对二维分数阶cable方程(2D-FCE): 一 (嘉+x,y,t)一x,y,t),n—Eo,L.-]XEo,].(5) 在空间与时间上离散得到逼近格式.一维分数阶cable方程的数值离散形式类似 地可以推得. 1.2用有限体积法(F,,M)离散空间导数 用如下网格离散空间变量: 五一ihx,hz一,一jhy,h一冠? 其中:h,h是空间步长.对偶网格为: - 27/I-1/2=(f+告)^,i==:0,1,2,…,一1. 辨l/2=(+?)^,.『一0,1,2,…,一1_ 对应于每个节点(五,)有一个控制体积,记此控制体积为: Vf.j一{(,)jX川z<<,1/2<Y<蜘1/2),i一1,2,…,M工;一1,2,…,. 且其边界包含四个逆时针方向的向量. 将式(5)左右两边在控制体积内进行积分得 J'j',出一盯(嘉+雾)如一?如(6) 记U(,)一u(x,,).对式(6)左边和右边的第二项分别进行如下近似: j'J'. ?.1l,j'J'"dxdy~.Ui,j(). 而对式(6)右边第一项应用Green公式有: J..『c嘉+孝出一』一考+ 用中点公式和一阶导数近似,有: 第2期沈淑君,等:有限体积法求解分数阶cable方程33 (嘉+雾)出?(考,靠,一a鱼yI(工.)hx+(az~----Ic一宝]x~va,Yj)? (丝吐一丝生)+(坠卜上一丝)一 凡v凡v凡凡. "n凡,,' 用上面的逼近式,对二维分数阶cable方程(5)在空间上的离散构成了如下近似形 式: ?D一(卜.()--2(,)+.())+.D(H)一2()+}卜.,,())一,().(7) 1.3离散时间分数阶导数 采用如下形式离散时间变量: 一O,1'2,…,N,r一. 其中:r是时19步长. 为简便起见,采用下面的近似: D~-rv(7/)dT/=南南 南? 式(8)第一项可以近似为: . d?d一)[(7--n73= 其中:系数一(s+1)一5. 式(8)第二项可以近似为: 南砖叻一而1曼s~O rt,砖竺?南n-1[一c"cb7一 而1[()一()][(n--5)一(n--s--1)]一EE()一口(?)]6,r 由此可得: j=rl出?而r)+而[((]一 而而(_6(.(9) 对空19离散形式式(7)两边从t到t井进行积分,有 Ui,(,井1)?仉j,j(,)+rlEui,一(,州)一2ui.j(,井1)+".升1(,升1)]+ n?(61一)[,产l(,,广)一2uf'J(,,广)+,(,一)]+r2[1.(,)一2ulo(,)+U/q-1,j(1)]+ 一 l一1 )[,,()一2(一)+)]一r3()一r3(一.』().(1o) r2(一 其中:rl一(1+ y1)'一南,一. 令"为u(x,,)的近似值,则可得到二维分数阶cable方程(5)的隐式格式为: ,r一1 一"+rl("掣一2;+urg-1..)+rl?(6,卜.一)(.一2z'+UW-S)+r2(n,+产l.一2苗+)+j;0 ,r1,r1 rz ?(6一)(,,一Zu7+"r州r-s,,)一r…Un+l—r3?(一)w-.(11),=Oj=O 类似地,可推得一维分数阶cable方程的隐式格式为: 泉州师范学院2011年3月 矿一u7+r2(一2+n+.1)+rz?(一)(一2u7-+)一 r3一r3?(一)uT-,.(12) 2数值例子 2.1一维情形 例1考虑如下分数阶cable方程"]: 丁au(x,t)=一)'0?z?1, t>0'] u(O,)一2,t>0,(13) u(x暑20'?Xj,O)一一z,??1.J 在图1中,取y1=72—1,它显示的是标准的反应一扩散过程,并且数值解很快就达 到稳定状态.而在图2 中,取y1一y2—0.5,由于反常扩散的影响,数值解需要花费比图1稍长的时间才能达到稳定状态.在图3和图 4中,当y1,y2取不同的数值时,解的性态与图1和图2有着显着的不同.图3显示当y1<y2时,数值解要经过 较长时间才能达到稳定状态.而在图4中,当y1>72时,经过一段时间的扩散,数 值解逐渐收敛回初始条件 u(x,O)一2一的状态.这些结果与文[14]中得到的解析解的结果一致. Distance 图1当y1一y2—1,t一0.1,1.0,10.O, 100.0时的数值解 图3当y1=0.5,y2—1.0,t一0.1, 1.0,10.0,i00.0时的数值解 DistanceJ=【 图2当一0.5,y2=0.5,f一0.1,1.0, 10.0,100.0时的数值解 图4当y1=1.0,y2=0.5,t一0.1, 1.0,10.0,100.0时的数值解 第2期沈淑君,等:有限体积法求解分数阶cable方程35 2.2二维情形 例2考虑下面的分数阶cable方程: 垒一D{(晏+寻)u(x,y,t)--(,,+厂(z,,,o?z,3,?,>o'1 u(O,y,)一(不,y,,)一0,0?y?不,t>0,}. u(x,0,,)=u(x,,,)一0,0?z?万,t>0,I u(x,y,O)一0,0?x,Y?兀J 其中:厂(x,y,t)一sin()sin()(2++). 式(14)的精确解是:u(x,y,)一t2sin(~)sin(ny).此例中,取y1—0.85,y2—0.9. O 雪0 卫 董o 0 图5当t一1.O,f一0.01,一hy=0.025时的数值解.图6当z一1.0时的精确解 图5显示的是,当t一1.O,y一0.01,=hy一0.025时,用本文所提出的数值方法求解 的结果.图6显示 的是当t:==1.0时的精确解.从中可见,图5的数值解与图6的精确解比较吻合. 参考文献: E1]RALIW.BranchingdendritictreesandmotoneuronmelTlbraneresistivity[J].ExpNeural,1959,1:491—527. [2]RALIw.TheNervousSystem,inHandbookofPhysiology[M].AmericanPhysiologicalSociety,Bethesda,MD,1977. [3]TUCKWELIHCIntroductiontoTheoreticalNeurobiology[M].Cambridge:CambridgeUniversityPress,1988. E4]WANGP,SONGQH,WANGJ,eta1.Researchoncableequationunderanydirectionalelectricalfield[J].JournalofTianjinUniver- sity,2004,37(6):515—519. [53QIANN,SEJNOW.SKITJ.Anelectro-di.usionmodelforcomputingm~'nbranepotentialsandiomcconcentrationsinbranching dendrites,spines,andaxons[J].BiolCybem,1989,62:1-15. [6]VANEGERAATJM,WIKSWOJP.Amodelforaxonalpropagationincorporatingbothradialandaxialiometransport[J].Biophys J,1993,64:1287—1298. [73WANGY.Amodifiedcablefunctionforreperesenttheexcitationofperipheralnervesbytransversefieldinducedbypulsedmagnetic field[c].Turkey,2001Proceedingsofthe23rdAnnualEMBSInternationalConference,2001;896-898. [8]KRZYZANSKIW.AnapplicationofaweaksolutionofthecableequationtotheRallmodelofanervecell[J].MathematicalMedi— cineandBiology,2001,18(4):377-407. [9]LINDSAYKA,ROSENBERGJR,TUCKERG.FromMaxwell'sequationstothecableequationandbeyond[J].ProgressinBio- physicsandMolecularBiology.2004,85(1):71-116. [103HENRYBI,LANGLANDSTA^,LFmctio~lcablemodelsforspinyneuronaldendrites[J].PhysicalReviewLetters,PRL100, 2008,128103. [11]LANGIANDSTAM,HENRYBI,WEARNESUFractionalcableequationmodelsforanomalouselectrodi.usioninnervecells: infinitedomainsolutions[J].JournalofMathematicalBiology,2009,59(6):761—808. [12]JIANGCZ,WANGH,WANGJ,eta1.Thesolvingandsimulationofcableequationunderthestimulationofpointelectricalsource [J].JournalofBiomedicalEngineering,2005,22(1):43—46. ?石4O?O0O0—0-1卫武,I 36泉州师范学院2011年3月 [13]WANGH,JIANGC乙 Resolutionofcontinuouscableequationbyintegraltransformationmethodanditseli'oranalysis[J].Jounlal 0fTianiinNormalUniversity:Natura1ScienceEdition,2004,24(3):56—58. [14]LANGIANDST,HENRYB,WEARNESSolutionsofafractionalcableequaiton:finitecase[R].AppliedMath髓1aticsReport AMR05/35,UniversityofNewSouthWales,2005. [15]LIUF,ANHV,TURNERI.NumericalsolutionofspacefractionalFokker-Planckequation[J].J.CompandApplMath,2004,166: 209—219. [163LIUF,ANHV,TURNERI,eta1.NumericalsimulationforsolutetransportinfractalporousmediaEJ].ANZIAMJ,2004,45(E): 461-473. [173ROoPJP.ComputationalaspectsofFEMapproximationoffractionaladvectiondispersionequationonboundeddomainsinR2口]. JCompAPPlMath,2006,193(1):243—268. [183LINY,XUCFinitedifference/spectralapproximationsforthetime-fractionaldiffusionequation[J].JCompPhys,2007,225(2): 1533-1552. [19]LILTQ,LrUF,TURNERI,eta1.ApproximationoftheLecy-Felleradvection-dispersionprocessbyrandomwalkandfinitediffer- eneemethod[J].JCompPhys,2007,222:57-70. [2o3ZHUANGP,LIUF,ANHV,eta1.Newsolutionandanalyticaltechniquesoftheimplicitnumericalmethodfortheanomalous sumiffusionequation[J].sIAMJNumerAnal,2008,46(2):1079—1095.9. FractionalCableEquationbyFiniteVolumeMethod SHENShu-jun,WUChun-hong2,LIUQing-xia. (1.SchoolofMathematicalSciences,HuaqiaoUniversity,Fujian362021,China; 2.DepartmentofMathematicsandPhysics,XiamenUniversityofTechnology,Fujian361024,China~ 3.SchoolofMathematicalSciences,XiamenUniversity,Fujian361005,China) Abstract:Thecableequationisoneofthemostfundamentalequationsforroodelingneuronaldynamics. The'fractiona1cableequationhasbeenusedtodescribenervefiber.Andithasobtainedbetterresultin simulatingneuron'Sdynamicsbehavior.Inthispaper,wemakethefirstattempttoapplythefinite volumemethodinthesolutionsofthefractionalcableequationinone-dimensionalcaseand two-dimensionalcase.FirstlythenumericalapproximationschemebyfinitevolumemethodiS presente&Then.twonumericalexamplesirMudingone-dimensionalcaseandtwo-dimensionalcaseare presentedbyusingthenumericaltechnique.whichshowstheefficiencyofthenumerica1method. Keywords:fractionalcableequation;numerica1approximation;finitevolumemethod.