扩散问题的偏微分方程模型,数学建模

1、 #(22) #(3) (22) (3) #第七节扩散问题的偏微分方程模型物质的扩散问题，在石油开采、环境污染、疾病流行、化学反应、新闻传播、煤矿瓦斯爆炸、农田墒情、水利工程、生态问题、房屋基建、神经传导、药物在人体内分布以及超导、液晶、燃烧等诸多自然科学与工程技术领域，十分普遍地存在着.显然，对这些问题的研究是十分必要的，其中的数学含量极大.事实上，凡与反应扩散有关的现象，大都能由线性或非线性抛物型偏微分方程作为数学模型来定量或定性地加以解决.MCM的试题来自实际，是“真问题数学建模计算机处理”的“三合一”准科研性质的一种竞赛，对上述这种有普遍意义和数学含量高，必须用计算机处理才能得到数值解

2、的扩散问题，当然成为试题的重要来源，例如，AMCM-90A,就是这类试题；AMCM-90A要研究治疗帕金森症的多巴胺(dopamine)在人脑中的分布，此药液注射后在脑子里经历的是扩散衰减过程，可以由线性抛物型方程这一数学模型来刻划.AMCM-90A要研究单层住宅混凝土地板中的温度变化，也属扩散(热传导)问题，其数学模型与AMCM-90A样，也是线性抛物型方程.本文交代扩散问题建模的思路以及如何推导出相应的抛物型方程，如何利用积分变换求解、如何确定方程与解的表达式中的参数等关键数学过程，且以AMCM-90A题为例，显示一个较细致的分析、建模、求解过程.1抛物型方程的导出设u(x,y,z,t)是

3、t时刻点(x,y,z)处一种物质的浓度.任取一个闭曲面S，它所围的区域是Q，由于扩散，从t到t+At时刻这段时间内，通过S流入Q的质量为M=F,AtJ!(a2竺cos,b2cos+c2竺cosy)dSdt.1t0 x0y0z由高斯公式得m,Atm1tQS/02%702%02%、1)2)(a2,b2,c2)dxdydzdt0 x20y20z2其中，a2,b2,c2分别是沿x,y,z方向的扩散系数.由于衰减(例如吸收、代谢等)，Q内的质量减少为M=Jt,At.Jk2udxdydzdt,2t其中k2是衰减系数.由物质不灭定律，在Q内由于扩散与衰减的合作用，积存于Q内的质量为M!-M2.换一种角度看，

4、Q内由于3深度之变化引起的质量增加为u(x,y,z,t,At)一u(x,y,z,t)dxdydzQ=Jt,AtJJJdydzdt.t0tQ显然M3=叫-M2，即T山JJJ竺dxdydzdtt&,B2uBz2-k2u)dxdydzdt.J“JJ(a2伽b2伽c2t0 x2By2,由At,t,之任意性得BuB2uB2uB2ua2+b2+c2一k2u(4)BtBx2By2Bz2方程(4)是常系数线性抛物型方程，它就是有衰减的扩散过程的数学模型，对于具体问题,尚需与相应的定解条件(初始条件与边界条件等)匹配才能求得确定情况下的解.2Dirac函数物理学家Dirac为了物理模型之需要，硬是引入了一个当时

5、颇遭微词的，使得数学与物理学传统密切关系出现裂痕的“怪”函数：u0,x丰0,feTOC o 1-5 h z6(x)x,00,x0,H(x)0,x0,如果沿用(6)中的算法，则在质量集中分布的这种情形有0,x丰x,P(x)o且Jxp(x)dxH(x-x),-g0于是得g,x0.7)Jgp(x)dx1.-g但是，从传统数学观点看，若一个函数除某点处处为零，则不论哪种意义下的积分，都必定为零，(7)式岂能成立！但是，6函数对于物理学而言是如此之有用，以致物理学家正当地拒绝放弃它.尽管当时数学家们大都嘲笑这种函数，但P.A.M.Dirac及其追随者们在物理领域却收获颇丰，Dirac于1933年获诺贝尔

6、物理奖.当然Dirac也意识到6(x)不是一个通常的函数，至于找一种什么办法来阐明6(x)这一符号的合法性，那就是数学家的任务了.1940年,法国数学家许瓦兹(L.Schwartz)严格证明了应用6(x)的正确性，把6函数置于坚实的0(x-x)f(x)0数学基础上；1950年，L.Schwartz获数学界最高奖Fields奖.函数的重要性质有：f+,(x-x)dx1-,0J+，(x-x)f(x)dxf(x).-,00其中f(x)GC(-，,+，)，即(x-x)摘出了f(x)在xx的值.00dH(x-x)显(x一x).dx0(x)的导数是存在的，不过要到积分号下去理解：+,(x-x)f(x)dx

7、-f(x),00J+，(n)(x-x)f(x)dx(-1)nf(n)(x).一，00事实上，由于(x-x)在+，-,处为零，则形式地用分部积分公式+，一+,(x-x)f(x)dx-,一，0=J+8(x-x)f(x)dx,一，0其中，f(x)GCn(-,+，)，于是有(11)与(12)公式.对于申(x)gC(-，,+，)，有p(x)(x-x)=p(x)(x-x).000(bx)=丄(x)(b丰0).|b|(x-x,y-y,z-z)=(x-x)(y-y)(z-z).000000付立叶变换F(yy)=e-叫.0F(x)=1.FC(x-x)+C(x-x)=CF(x-x)+CF(x-x).1122112

8、29)拉普拉斯变换F(x-x)=e-x0,F(x)=1.0FC(x-x)+C(x-x)=CF(x-x)+CF(x-x).11221122从上面的定义与性质看出，Delta函数(x)与一般可微函数还是有重大区别的,它是“广义函数.”3Cauchy问题的解设扩散源在点(x0,y0,z0)处，则此扩散问题满足Cauchy问题dud2urd2ud2u7a2+b2+c2一k2U,dtdx2dy2dz2u(x,y,z,0)M(x一x)(y一y)(z一z).(8)(9)10)11)12)(13)(14)(15)(16)(17)(18)(19)(20)我们说(21)000对(21)(22)进行付立叶变换，且令

9、=(,),u(t,)二Fu(x,y,z,t),123由于，2u,2u,2u、F=2u,F=2u,F=2u,x21,y22，z23Fu(x,y,z,0)=MF(x-x)F(y-y)F(z-z)000=Me-i(Vo+码yo+人z0),故得常微分方程Cauchy问题du+(a2+b22+c22+k2)u=0,dt123u(0,)=Me(x0+2y0+3z0).得唯一解23)u(t,)=Me-(a22+b22+c22+k2)t-i(1x0+2y0+3z0).对(23)求逆变换F-1，由于-12aF-1e4a2=e-ax2，PF-1e-iix0-1e4a2=e-a2(x-x)2，兀0故得u(x,y,z

10、,t)=F-而M(2.8tabc7|4a2t4b2t4c2t如果认为经过了相当长时间后，扩散已经终止，物质分布处于平衡状态，则方程(4)24) (22) #(3) #,u中的芬=0，于是有线性椭圆型方程的边值问题a2宜+b2取+c2包-k2u=0,(x,y,z)eD,x2,y2，z2u(x,y,z)=申(x,y,z).Il,D=(x,y,z)也可以用付立叶变换求解.当然，根据实际情况，还可以考虑第二边条件,u,n,D,u或第二边条件a+u=p(x,y,z)等，其中,D是区域D的边界，n是外法线方向,n，da,是实常数. # (3) #4参数估计在Cauchy问题(21)(22)的解(23)中，

11、有四个未知的参数a,b,c,k,它们分别是扩散与衰减过程中的扩散系数与衰减系数的算术平方根.至于点源的质量与位置M,(x,y,z是已知的.设观测取样为：(X,y,z,m),(x,y,z,m)，,(x,y,z,m),00011112222nnnn取样时刻为t1(不然设ttt,t是取样时间，则(21)变成Uta2U+tb2U+言，00t0 xx0yytc2U-1k2U,对工而言，取样时间为1,而方程形状与(21)一致)，把在(x,y,z)点0zz0i观测到的物质密度m与公式(24)都取对数，令t1，则iM(x-x)2(y-y)2(z-z)225)lnu(x,y,z,1)ln-lnabc-o+o+o

12、+k2.TOC o 1-5 h z(2折)34a24b24c2人Y(xx)2v(yy)27(zz)2101v1444a2b2c2Mln-lnabc-k2，则(25)写成(2航)326)Wlnu(x,y,z,1)二X+0Y+YZ,而我们已观测得(X,Y,Z,W)i1,2,n的数据，用三元回归分析方法求出,0,1,的估计值如下27)其中W1艺nk1W,kXL工x,Y1工y,Z1工nik1nik1Z,nik1 # #(3) # (3) #i,10+111,TOC o 1-5 h z111213101,+10+111,212223201,+10+111.3330I3132一一其中l=工(Wkk11(z

13、-Z)(W-W),30kkk11=工(X-X)2,1=工(Y-Y)2,1k22kk1k11=工(x-X)(Y-Y),kkk11=工(Y-Y)(Z-Z),23kk-X)(w-W),k=工(Y-Y)(Wkl20kk1l13k1ll,2112k1-W),工(Z-Z)2,33kk1(x-X)(z-z),kll,3113ll.3223TOC o 1-5 h z111由B,V可求得a2,b2,c2的估计值，即a2=一=,b2=,c2.又由于p,28)k2=8+InabcInM，(2折)3由（27）式可得，再把a,b,c代入（28）得:M29)k2=s+Inabc一In（2麻）3至此得到参数a2,b2,c2

14、,k2的估计值a2,b2,C2,k2，把它们代入（24）分别替代a2,b2,c2,k2,则得不含未知参数的解u（x,y,z,t）的近似表达式.5竞赛试题分析AMCM-90A不可用本文的思路与方法加以解决；该试题由东华盛顿大学数学系YvesNievergelt提供，要求研究药物在脑中的分布，题文称：“研究脑功能失调的人员欲测试新的药物的效果，例如治疗帕金森症往脑部注射多巴胺（Dopamine）的效果，为了精确估计药物影响到的脑部区域，它们必须估计注射后药物在脑内空间分布区域的大小和形状.“研究数据包括50个圆柱体组织样本的每个样本药物含量的测定值（如图6-1），每个圆柱体长0.76mm,直径0.

15、66mm,这些互相平行的圆柱体样本的中心位于网格距为lmX0.7XmX1mm的格点上，所以圆柱体互相间在底面上接触，侧面互不接触.注射是在最高计数的那个圆柱体的中心附近进行的.自然在圆柱体之间以及由圆柱体样本的覆盖的区域外也有药物.“试估计受到药物影响的区域中药物的分布.”“一个单位表示一个闪烁微粒的计数，或多巴胺的4.75X10-18克分子量，例如表6-1指出位于后排当中那个圆柱体的含药量是28353个单位.”后方垂直截面16444213204141884807022144115158352209123027283531313868178921260209211173172721313033

16、7651715453前方垂直截面1633244322431667121055609810482322137155311974247853304441143114960318230129420611036258188图6-1数学模型只是实际问题的近似，要建立数学模型，一般首先要对所研究的实际问题进行必要和允许的简化与假设，而且，不同的简化与假设，又可能导致不同的数学模型，例如2是抛物型方程模型，而3则是椭圆方程模型.假设：注射前大脑中的多巴胺含量可以忽略不计.大脑中多巴胺注射液经历着扩散与衰减的过程，且沿x,y,z三个方向的扩散系数分别是常数，衰减使质量之减少与深度成正比.注射点在后排中央那个圆

17、柱中心，即注射点的坐标(x,y,z)已知，注射量有医000疗记录可查，是已知的.注射瞬间完成，可视为点源delta函数.取样也是瞬间完成，取样时间已知为t1.样本区域与整个大脑相比可以忽略，样本组织远离脑之边界，不受大脑边界面的影响.在以上假设之下，显然可以用本文前面讲过的思路来建模，于是得AMCM-90A的数学模型为Cauchy问题(21)(22),解的表达式为(24)，且用三元回归分析来估出参数a,b,c,k，于是可以求得任意位置任意时刻药物的深度.如果所给数据认为是在平衡状态测得的，药物注射进脑后，从高深度处向低深度处扩散，与扩散同时，一部分药物进入脑细胞被吸收固定，扩散系数与吸收系数都

18、是常数，但过一段时间，所有药物都被脑细胞所固定，达到了平衡态.在这种假设下，3给出了下述的分析、建模、求解过程.设v(x,y,z,t)是t时刻在(x,y,z)点处游离的药物浓度，w(x,y,z,t)是t时刻(x,y,z)点处吸收固定的药物浓度，u(x,y,z)是达到平衡态时(x,y,z)点处吸收固定的药物浓度.又设游离药物在各方向上有相同的扩散系数k，吸收系数为h，于是有k,v一hvat又石=hv，即吸收速度与游离的浓度成正比，代入(30)得dvatka/人、aw(,w)-hatat对(31)关于t从0到+8积分得v+8t0+8一wt0+8t030)31)32) # (3) #由于最后无游离药

19、物，故v(x,y,z,+8)0，又开始时(t0)无被吸收的药物，故w(x,y,z,0)0,w(x,y,z,0)=0；平衡状态在t=+8时达到，这时u(x,y,z)=w(x,y,z,+8)，于是由(32)得33)v(x,y,z,0)其中v(x,y,z,0)是开始时的浓度分布，近似于注射点的点源脉冲函数.把此注射点取为坐标(、k原点(0,0,0)，则v(x,y,z,0)L(x,y,z),L是注射量，于是记2二Ih丿34)y2,u+uL(x,y,z)，作付立叶变换得2(g叶其孝，大学生数学建模竞赛辅导教材，湖南教育出版社，1993.Christopher,R.Malone,GianPauletto,James,I.Zoellick,DistributionofDopamineintheBrain,TheJournalofUndergraduateMathematics,anditsApplications,vol.12(1991),SpecialIssue:The1991MathematicalContestinModeling,pp.211-223.孙晓东，荆秦，梁俊，脑中药物分布的数学模型，数学的实践与认识，1991年No.4,63-69.中国科学院数理统计组，常用数理统计方法，科学出版社，19784