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

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

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

3、)处一种物质的浓度.任取一个闭曲面S,它所围的区域是Q,由于扩散,M1从t到t+&时刻这段时间内，通过S流入C的质量为tT2Fu一2Fu二2Fu=(acosq一bcos-ccos)dSdt.tS二xcy二z由高斯公式得-2.2.2M1/2二U2二U2二U(1)!.!.!.(abc-)dxdydzdt,1x:y二z其中，a2,b2,c2分别是沿x,y,z方向的扩散系数.由于衰减(例如吸收、代谢等)，Q内的质量减少为t2,M2=I0kudxdydzdt,(2)其中k2是衰减系数.由物质不灭定律，在建内由于扩散与衰减的合作用，积存于C内的质量为M1-M2.换一种角度看，。内由于深度之变化引起

4、的质量增加为M3：iiiu(x,y,z,tt)-u(x,y,z,t)dxdydzQttt：udxdydzdt.(3)t二t显然M3=MiM2,即文案大全标准实用dzdt._2_2_2t/21u2；-u2二u2=(a-2b-2c-2-ku)dxdydzdt.ti!:x.:yjz由AttC之任意性得.2.2.2.U2二u2u21u2二a2b2c-2-ku2222:t二x二y二z方程(4)是常系数线性抛物型方程，它就是有衰减的扩散过程的数学模型，尚需与相应的定解条件(初始条件与边界条件等)匹配才能求得确定情况下的解对于具体问题,§2Dirac函数物理学家Dirac为了物理模型之需要，硬是引

5、入了一个当时颇遭微词的,理学传统密切关系出现裂痕的“怪”函数：使得数学与物；、(x)dx=1.(5)它的背景是清晰的，以一条无穷长的杆子为例，沿杆建立了一维坐标系，点的坐标为的线密度是P(x),在(*,x段，杆子质量为m(x),则有dm(x)x二.(x),(x)dx=m(x).dx-二设此无穷长的杆子总质量为1,质量集中在x=x0点，则应有1,xXo,、一、m(x)=i或与成m(x)=H(xXo),0,x:Xo,其中H(x)为(6)1,x0,H(x)=',0,x二0,如果沿用(6)中的算法，则在质量集中分布的这种情形有、O,X=X0,：(x)二0二,x=0.x且JP(x)dx=H(x-

6、x0),于是得.二：(x)dx=1.但是，从传统数学观点看，若一个函数除某点处处为零，则不论哪种意义下的积分，都必定为零，(7)式岂能成立！但是，6函数对于物理学而言是如此之有用，以致物理学家正当地拒绝放弃它.尽管当时数学家们大都嘲笑这种函数，但P.A.M.Dirac及其追随者们在物理领域却收获颇丰，Dirac于1933年获诺贝尔物理奖.当然Dirac也意识到6(x)不是一个通常的函数，至于找一种什么办法来阐明6(x)这一符号的合法性，那就是数学家的任务了1940年，法国数学家许瓦兹(L.Schwartz)严格证明了应用6(x)的正确性，把6函数置于文案大全标准实用坚实的数学基础上；1950年

7、，L.Schwartz获数学界最高奖Fields奖.6函数的重要性质有：1) f6(x-x0)dx=1.(8)2) 6(xx0)f(x)dx=f(x0).(9)其中f(x)wC(g，+w),即6(xx0)摘出了f(x)在x=x0的值.c、dH(x-x0),、3) =6(xXo).(10)dx4) 6(x)的导数是存在的，不过要到积分号下去理解：(6'(x-Xo)f(x)dx=-f'(Xo),(11)-cd(x-Xo)f(x)dx=(1)nf(Xo).(12)事实上，由于S(x-x0)在代0,-00处为零，则形式地用分部积分公式(X-Xo)f(x):一；(x一Xo)f(x)dx=

8、；(x-x0)f(x)dx,其中，f(x)wCn(-OCI,z),于是有(11)与(12)公式.5)对于中(x)wC(-°o,有(x)、(x-x0)=(Xo)、(x-Xo).(13)一1|b|(14)7)S(xXo,yy0,zzo)=6(xxo)6(yyo)6(zzo).(15)8)付立叶变换F6(yy0)=e灰.(16)F5(x)=1.(17)FC(x-x1)+C26(x-x2)=C1F6(xx1)+C2F6(x-x2).(18)9)拉普拉斯变换F6(x-Xo)=eiF6(x)=1.(19)FC(x-X1)+C2WX-X2)=CF6(xX1)+C2F6(x-X2).(20)6)(b

9、x)=(x)(b；0)从上面的定义与性质看出，Delta函数6(x)与一般可微函数还是有重大区别的，我们说它是“广义函数.”§3Cauchy问题的解设扩散源在点(Xo。，。)处，则此扩散问题满足Cauchy问题卜-u二aft(21)(22)2C2u,(2C2u,2*u,2-2b-2c-2-ku，:xcy；zu(x,y,z,0)=M'(xx。)'(y-y。)、(z-z。).文案大全标准实用对(21)(22)进行付立叶变换，且令，%)，i?(t,Q=Fu(x,y,z,t),由于.2.2.2U,2U2U2qF2二一'1RF2二一'2U,F2二一'3U

10、,jxfyjzFu(x,y,z,0)=MF、.(x-x0)F、(y-y0)FU(z-Zo)=Me(-y0%z0)故得常微分方程Cauchy问题*(a2；b2c2k2)U?=0,U?(0,)=Me*ix°'2y0酒).得唯一解U(t九)"Me'/兑也2冗*2,也2*(?a%丫0岁34)对(23)求逆变换f，由于(23)故得Fe4a.2_Leqx°e4a2i'22.ax=e，5：-2,、2=>e(x-Xo),5U(x,y,z,t)=FR一MLexp_(x-x0)2_(y-y0)2_(z-z0)22t(2,H)3.a2b2c2p.4a2t4

11、b2t4c2t:Mexp-3-"-0-k2t.8二tabc；t4at4bt4ct如果认为经过了相当长时间后，扩散已经终止，物质分布处于平衡状态，则方程(24)4)二u中的上=0,于是有线性椭圆型方程的边值问题；:t-221U.2xb2.2二U-2-y.220U,2八c-2-kU=0,£-1z(x,y,z)Du(x,y,z)D=/(x,y,z).也可以用付立叶变换求解当然，根据实际情况，还可以考虑第二边条件-(x,y,z)或第三边条件：%P是实常数.+Pu=P(x,y,z)等，其中£D是区域D的边界,n是外法线方向,文案大全标准实用§4参数估计在Cauch

12、y问题(21)(22)的解(23)中，有四个未知的参数a,b,c,k,它们分别是扩散与衰减过程中的扩散系数与衰减系数的算术平方根.至于点源的质量与位置M,(X0,y0,Z0)是已知的.设观测取样为：(xi,yi,zi,mi),(X2,y2,Z2,m2),|,(Xn,yn,Zn,mn),取样时刻为t=1(不然设t=toT,to是取样时间，则(21)变成Ut=toa2Uxx+tob2Uyy+,22.,、10cUzz10kU,对T而言，取样时间为1,而方程形状与(21)一致)，把在(Xi,yi,Zi)点观测到的物质密度mj与公式(24)都取对数，令t=1,则2221/、1M1,r(x-x0)(y-y

13、0)(zy),2n(25)lnu(x,y,z,1)=lnInabc-0-；k.(2、二)34a24b24c2,、2,、2,、2x=_,Y=lsyLz=l£2zl_4;=，444abc,M(2；)32-Inabc-k，则(25)与成W=Inu(x,y,z,1)=二X而我们已观测得估计值如下：(Xi,Y,Zi,Wi)i=1,2JH,n的数据，用三元回归分析方法求出?=W_(MX十够十.?Z),(26)s,p,y,名的(27)其中1/W=一、Wk,nkm-1X一Xi,Y,nkd一1nZZi,nkmfl11<?,+|12+l13?l10,4#,+l22肾+I23'?=l20,l

14、31区,*323*33?=l30.其中文案大全标准实用l10l30l11n八(Wk-又)(Wk.W),k1nC(Zk-Z)(Wk.W),k1nn-2=、(Xk-X),I22="n_I20八(Yk-Y)(Wk.W),k4n22(Yk-Y),133c(Zk-Z)k=1nI12八(Xk-X)(Yk-Y),k1n_123c(Yk-Y)(Zk-Z),k1nI13(Xk-X)(Zk-Z),k421-l12,131=l13,l32-l23-，一c一1c1c1又由于由鱼？,？可求得a2,b2,c2的估计值，即？=,萨=丁，然.?,2一，M(28)(29)k=s+lnabcln1=r-,（2二）3由（

15、27）式可得？，再把？,b,a代入（28）得i?=?lnabl?-ln-M_.（2.二）至此得到参数a2,b2,c2,k2的估计值，仔,，!？，把它们代入（24）分别替代a2,b2,c则得不含未知参数的解u（x,y,z,t）的近似表达式.§5竞赛试题分析AMCM-90A不可用本文的思路与方法加以解决；该试题由东华盛顿大学数学系丫vesNievergelt提供，要求研究药物在脑中的分布，题文称：“研究脑功能失调的人员欲测试新的药物的效果，例如治疗帕金森症往脑部注射多巴胺（Dopamine）的效果，为了精确估计药物影响到的脑部区域，它们必须估计注射后药物在脑内空间分布区域的大小和形状.“

16、研究数据包括50个圆柱体组织样本的每个样本药物含量的测定值（如图6-1）,每个圆柱体长0.76mm,直径0.66mm,这些互相平行的圆柱体样本的中心位于网格距为1mm<0.76xmm<1mmi勺格点上，所以圆柱体互相间在底面上接触，侧面互不接触.注射是在最高计数的那个圆柱体的中心附近进行的.自然在圆柱体之间以及由圆柱体样本的覆盖的区域外也有药物.“试估计受到药物影响的区域中药物的分布.”“一个单位表示一个闪烁微粒的计数，或多巴胺的4.753X10-18克分子量，例如表6-1指出位于后排当中那个圆柱体的含药量是28353个单位.”后方垂直截面164442132041418848070

17、221441151583522091230272835313138681789212602092111731727文案大全标准实用213130337651715453刖力垂直截回1633244322431667121055609810482322137155311974247853304441143114960318230129420611036258188图6-1数学模型只是实际问题的近似，要建立数学模型，一般首先要对所研究的实际问题进行必要和允许的简化与假设，而且，不同的简化与假设，又可能导致不同的数学模型，例如2是抛物型方程模型，而3则是椭圆方程模型.假设：(1)注射前大脑中的多巴胺含量

18、可以忽略不计(2)大脑中多巴胺注射液经历着扩散与衰减的过程，且沿x,y,z三个方向的扩散系数分别是常数，衰减使质量之减少与深度成正比(3)注射点在后排中央那个圆柱中心，即注射点的坐标(x0,y0,%)已知，注射量有医疗记录可查，是已知的.(4)注射瞬间完成，可视为点源delta函数.(5)取样也是瞬间完成，取样时间已知为t=1.(6)样本区域与整个大脑相比可以忽略，样本组织远离脑之边界，不受大脑边界面的影响.在以上假设之下，显然可以用本文前面讲过的思路来建模，于是得AMCM-90的数学模型为Cauchy问题(21)(22),解的表达式为(24),且用三元回归分析来估出参数a,b,c,k,于是可

19、以求得任意位置任意时刻药物的深度如果所给数据认为是在平衡状态测得的，药物注射进脑后，从高深度处向低深度处扩散，与扩散同时，一部分药物进入脑细胞被吸收固定，扩散系数与吸收系数都是常数，但过一段时间，所有药物都被脑细胞所固定，达到了平衡态.在这种假设下，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,于是有(30)(31)A一二kvhv.ft又竺=hv,即吸

20、收速度与游离的浓度成正比，代入(30)得Ftvk：:w=(w)-.fthft：t对(31)关于t从。到积分得文案大全标准实用5k.5-hr-vt_ot&wt_owt-0-(32)一h一一由于最后无游离药物，故v(x,y,z,+厘)=0,又开始时(t=0)无被吸收的药物，故w(x,y,z,0)=0,Aw(x,y,z,0)=0；平衡状态在t=十°0时达到，这时u(x,y,z)=w(x,y,z,十无)，于是由(32)得k.-Au+u=v(x,y,z,0),(33)h其中v(x,y,z,0)是开始时的浓度分布，近似于注射点的点源脉冲函数.把此注射点取为坐标原点(0,0,0),则v(x

21、,y,z,0)=L6(x,y,z),L是注射量，于是f记订2=上|h2.-aAu+u=L6(x,y,z),(34)作付立叶变换得二2(22s2)?i?=L,L1(22s2),再作反变换得(35)u=LC1exp-cx2y2z2x2y2z2其中C是可计算常数如果考虑各向不同性，设x,y,z方向上扩散系数分别为a2,b2,c：注射点在(xcy0,z0),则2-2-2、-a2®+b2®+c2W+u=L6(x-%)6(y-丫0)82-4),：£xy也)于是解为u(x,y,z)=DL121212.看刈)b2(y-y0)c2(Z-Z0)Lexp1-.21.21.2x-xo)+(y-y0)