(完整版)扩散问题的偏微分方程模型,数学建模

第七节 扩散问题的偏微分方程模型
物质的扩散问题，在石油开采、环境污染、疾病流行、化学反应、新闻传播、煤矿瓦斯爆炸、农田墒情、水利工程、生态问题、房屋基建、神经传导、药物在人体内分布以及超导、液晶、燃烧等诸多自然科学与工程技术领域，十分普遍地存在着. 显然，对这些问题的研究是十分必要的，其中的数学含量极大. 事实上，凡与反应扩散有关的现象，大都能由线性或非线性抛物型偏微分方程作为数学模型来定量或定性地加以解决.
ＭＣＭ的试题来自实际，是“真问题⊕数学建模⊕计算机处理”的“三合一”准科研性质的一种竞赛，对上述这种有普遍意义和数学含量高，必须用计算机处理才能得到数值解的扩散问题，当然成为试题的重要来源，例如，ＡＭＣＭ-90Ａ，就是这类试题；ＡＭＣＭ-90Ａ要研究治疗帕金森症的多巴胺（dopamine ）在人脑中的分布，此药液注射后在脑子里经历的是扩散衰减过程，可以由线性抛物型方程这一数学模型来刻划. ＡＭＣＭ-90Ａ要研究单层住宅混凝土地板中的温度变化，也属扩散（热传导）问题，其数学模型与ＡＭＣＭ-90Ａ一样，也是线性抛物型方程.
本文交代扩散问题建模的思路以及如何推导出相应的抛物型方程，如何利用积分变换求解、如何确定方程与解的表达式中的参数等关键数学过程，且以ＡＭＣＭ-90Ａ题为例，显示一个较细致的分析、建模、求解过程.
§1 抛物型方程的导出
设(,,,)u x y z t 是t 时刻点(,,)x y z 处一种物质的浓度. 任取一个闭曲面S ，它所围的区域是Ω，由于扩散，从t 到t t +∆时刻这段时间内，通过S 流入Ω的质量为
2
221(cos cos cos )dSd t t
t
S
u u u M a b c t x y z
αβγ+∆∂∂∂=++∂∂∂⎰
⎰⎰. 由高斯公式得
2222
221222()d d d d t t
t
u u u M a b c x y z t x y z +∆Ω
∂∂∂=++∂∂∂⎰
⎰⎰⎰. （1） 其中，2
2
2
,,a b c 分别是沿,,x y z 方向的扩散系数. 由于衰减（例如吸收、代谢等），Ω内的质量减少为
2
2d d d d t t
t
M k u x y z t +∆Ω
=⎰
⎰⎰⎰， （2） 其中2
k 是衰减系数.
由物质不灭定律，在Ω内由于扩散与衰减的合作用，积存于Ω内的质量为12M M -.
换一种角度看，Ω内由于深度之变化引起的质量增加为
3[(,,,)(,,,)]d d d d d d d . (3)t t
t
M u x y z t t u x y z t x y z
u
x y z t t Ω
+∆Ω
=+∆-∂=∂⎰⎰⎰⎰
⎰⎰⎰
显然312M M M =-，即
222
2222222d d d d ()d d d d .t t
t
t t
t
u
x y z t t u u u a b c k u x y z t x y z
+∆Ω
+∆Ω
∂∂∂∂∂=++-∂∂∂⎰
⎰⎰⎰⎰
⎰⎰⎰
由,,t t ∆Ω之任意性得
222
2222
222u u u u a b c k u t x y z
∂∂∂∂=++-∂∂∂∂ （4） 方程（４）是常系数线性抛物型方程，它就是有衰减的扩散过程的数学模型，对于具体问题，尚需与相应的定解条件（初始条件与边界条件等）匹配才能求得确定情况下的解.
§2 Ｄirac 函数
物理学家Dirac 为了物理模型之需要，硬是引入了一个当时颇遭微词的，使得数学与物理学传统密切关系出现裂痕的“怪”函数：
0,0,() ()1.,0,x x x dx x δδ+∞
-∞
≠⎧==⎨∞=⎩⎰ （5）
它的背景是清晰的，以一条无穷长的杆子为例，沿杆建立了一维坐标系，点的坐标为x ，杆的线密度是()x ρ，在(,]x -∞段，杆子质量为()m x ，则有
d ()
(), ()d ().d x m x x x x m x x
ρρ-∞==⎰. （6）
设此无穷长的杆子总质量为１，质量集中在0x x =点，则应有
001,,()0,,
x x m x x x >⎧=⎨<⎩ 或写成 0()()m x H x x =-，
其中()H x 为
1,0,
()0,0,
x H x x >⎧=⎨
<⎩ 如果沿用（6）中的算法，则在质量集中分布的这种情形有
00,,
(),0.
x x x x ρ≠⎧=⎨
∞=⎩
且
0()d ()x
x x H x x ρ-∞
=-⎰
，
于是得
()d 1.x x ρ+∞
-∞
=⎰
. （7）
但是，从传统数学观点看，若一个函数除某点处处为零，则不论哪种意义下的积分，都必定为零，（7）式岂能成立！但是，δ函数对于物理学而言是如此之有用，以致物理学家正当地拒绝放弃它. 尽管当时数学家们大都嘲笑这种函数，但Ｐ.Ａ.Ｍ.Ｄirac 及其追随者们在物理领域却收获颇丰，Dirac 于1933年获诺贝尔物理奖. 当然Dirac 也意识到()x δ不是一个通常的函数，至于找一种什么办法来阐明()x δ这一符号的合法性，那就是数学家的任务了. 1940年，法国数学家许瓦兹（L.Schwartz ）严格证明了应用()x δ的正确性，把δ函数置于坚实的
数学基础上；1950年，L. Schwartz 获数学界最高奖Fields 奖.
δ函数的重要性质有：
1）0()d 1x x x δ+∞
-∞
-=⎰. （8）
2）00()()d ()x x f x x f x δ+∞
-∞
-=⎰
. （9）
其中()(,)f x C ∈-∞+∞，即0()x x δ-摘出了()f x 在0x x =的值.
3）
00()
()dH x x x x dx
δ-=-. （10）
4）()x δ的导数是存在的，不过要到积分号下去理解：
00()()(),x x f x dx f x δ+∞
-∞''-=-⎰ （11）
()()00()()(1)().n n n x x f x dx f x δ+∞
-∞
-=-⎰
（12）
事实上，由于0()x x δ-在,+∞-∞处为零，则形式地用分部积分公式
000()()
()()d ()()d ,
x x f x x x f x x
x x f x x δδδ+∞
+∞
-∞
-∞
+∞-∞
'---'=-⎰⎰
其中，()(,)n
f x C ∈-∞+∞，于是有（11）与（12）公式.
5）对于()(,)x C ϕ∈-∞+∞，有
000()()()()x x x x x x ϕδϕδ-=-. （13）
6）1
()() (0)||
bx x b b δδ=≠. （14）
7）000000(,,)()()()x x y y z z x x y y z z δδδδ---=---. （15）
8）付立叶变换
00[()].i x y y e λδ--= （16）
[()] 1.x δ= （17）
11221122[()()][()][()].C x x C x x C x x C x x δδδδ-+-=-+- （18）
9）拉普拉斯变换
00[(),[() 1.x x x e x δδδ--== （19）
11221122[()]()][()[()].C x x C x x C x x C x x δδδδ-+-=-+- （20） 从上面的定义与性质看出，Delta 函数()x δ与一般可微函数还是有重大区别的，我们说
它是“广义函数. ”
§3 Cauchy 问题的解
设扩散源在点000(,,)x y z 处，则此扩散问题满足Cauchy 问题
222
2222
222
000, (21)(,,,0)()()(). (22)
u u u u a b c k u t
x y z u x y z M x x y y z z δδδ⎧∂∂∂∂=++-⎪∂∂∂∂⎨⎪=---⎩
对（21）（22）进行付立叶变换，且令
123ˆ(,,), (,)[(,,,)]u
t u x y z t λλλλλ==， 由于
22222
2123222ˆˆˆ[], [], [],u u u u
u u x y z
λλλ∂∂∂=-=-=-∂∂∂ 102030000()[(,,,0)]
[()][()][()] ,
i x y z u x y z M x x y y z z Me λλλδδδ-++=---= 故得常微分方程Cauchy 问题
1020302222222
123()ˆ()0,ˆ(0,).i x y z du a b c k u
dt
u Me
λλλλλλλ-++⎧++++=⎪⎨⎪=⎩ 得唯一解
2222222123102030()()
ˆ(,)a b c k t i x y z u
t Me λλλλ
λλλ-+++-++=. （23）
对（23）求逆变换
1
-，由于
2
122
2
1
4[]a x
a e
λ-
--=， 2
110
2
2
1
240[]()i x e a
a e
x x λλ--
--=
-， 故得
1
2222000222
ˆ(,,,)[]()()()exp 444
u x y z t u x x y y z z k t a t b t c t -=
⎧⎫---=
----⎨⎬⎩⎭
2222000222
()()().444x x y y z z k t a t b t c t ⎧⎫---=----⎨⎬⎩⎭
（24） 如果认为经过了相当长时间后，扩散已经终止，物质分布处于平衡状态，则方程（4）中的
0u
t
∂=∂，于是有线性椭圆型方程的边值问题 222
2222222
0, (,,)(,,)(,,).
D u u u a b c k u x y z D x
y z u x y z x y z ϕ∂⎧∂∂∂++-=∈⎪∂∂∂⎨⎪=⎩
也可以用付立叶变换求解. 当然，根据实际情况，还可以考虑第二边条件(,,)D
u
x y z n ∂∂=ψ∂或第三边条件[](,,)D u
u x y z n
α
βρ∂∂+=∂等，其中D ∂是区域D 的边界，n 是外法线方向，,αβ是实常数.
§4 参数估计
在Cauchy 问题（21）（22）的解（23）中，有四个未知的参数,,,a b c k ，它们分别是扩散与衰减过程中的扩散系数与衰减系数的算术平方根. 至于点源的质量与位置
000,(,,)M x y z 是已知的.设观测取样为：11112222(,,,), (,,,),,(,,,),n n n n x y z m x y z m x y z m
取样时刻为1t =（不然设00, t t t τ=是取样时间，则（21）变成2200t xx yy U t a U t b U =++
2200zz t c U t k U -，对τ而言，取样时间为１,而方程形状与（21）一致），把在(,,)i i i x y z 点
观测到的物质密度i m 与公式（24）都取对数，令1t =，则
2222
000222()()()ln (,,,1)ln []444x x y y z z u x y z abc k a b c ---=--+++. （25） 令222000222()()()111
,,,,,,444x x y y z z X Y Z a b c αβγ---=
===-=-=-
2
ln ln abc k ε=--，则（25）写成 ln (,,,1)W u x y z X Y Z αβγε==+++，
（26） 而我们已观测得(,,,)1,2,,i i i i X Y Z W i n =的数据，用三元回归分析方法求出,,,αβγε的
估计值如下：
ˆˆˆˆ()W X Y Z εαβγ=-++， （27）
其中
1111
1111, , , ,n n n n
k i i i k k k k W W X X Y Y Z Z n n n n ========∑∑∑∑
ˆˆˆ,,α
βγ满足方程组 11121310
2122232031323330ˆˆˆ,,ˆˆˆ,,ˆˆˆ,.
l l l l l l l l l l l l αβγαβγαβγ⎧++=⎪⎪++=⎨⎪++=⎪⎩ 其中
102011
3012
2
211223311
1
121311
23211
()(), ()(),
()(),
(), (), (),()(), ()(),
()(), n n
k k k k k k n
k k k n
n n
k k k k k k n
n
k k k k k k n
k k k l W X W W l Y Y W W l Z Z W W l X X l Y Y l Z Z l X X Y Y l X X Z Z l Y Y Z Z l ==========--=--=--=-=-=-=--=--=--∑∑∑∑∑∑∑∑∑1231133223, , .
l l l l l ===
由ˆˆˆ,,α
βγ可求得2
2
2
,,a b c 的估计值，即222111
ˆˆˆ, , ˆa b c
αβγ
=-=-=-. 又由于 2
ln k abc ε=+- （28） 由（27）式可得ˆε
，再把ˆˆˆ,,a b c 代入（28）得 2
ˆˆˆˆˆln k
abc ε=+- （29）
至此得到参数2222,,,a b c k 的估计值2222ˆˆˆˆ,,,a b c k ，把它们代入（24）分别替代2222,,,a b c k ，
则得不含未知参数的解(,,,)u x y z t 的近似表达式.
§5 竞赛试题分析
AMCM-90A 不可用本文的思路与方法加以解决；该试题由东华盛顿大学数学系Ｙves Nievergelt 提供，要求研究药物在脑中的分布，题文称：
“研究脑功能失调的人员欲测试新的药物的效果，例如治疗帕金森症往脑部注射多巴胺（Dopamine ）的效果，为了精确估计药物影响到的脑部区域，它们必须估计注射后药物在脑内空间分布区域的大小和形状.
“研究数据包括50个圆柱体组织样本的每个样本药物含量的测定值（如图6-1），每个圆柱体长0.76mm ，直径0.66mm ，这些互相平行的圆柱体样本的中心位于网格距为1m m ×0.76×m m ×1mm 的格点上，所以圆柱体互相间在底面上接触，侧面互不接触. 注射是在最高计数的那个圆柱体的中心附近进行的. 自然在圆柱体之间以及由圆柱体样本的覆盖的区域外也有药物.
“试估计受到药物影响的区域中药物的分布. ”
“一个单位表示一个闪烁微粒的计数，或多巴胺的4.753×10-18克分子量，例如表6-1指出位于后排当中那个圆柱体的含药量是28353个单位. ”
后方垂直截面
164
442 1320 414 188 480 7022 14411 5158 352 2091 23027 28353 13138 681 789 21260 20921 11731 727 213 1303
3765
1715
453
前方垂直截面
163 324 432 243
166 712 1055 6098 1048 232 2137 15531 19742 4785 330 444 11431 14960 3182 301 294
2061
1036 258
188
图6-1
数学模型只是实际问题的近似，要建立数学模型，一般首先要对所研究的实际问题进行必要和允许的简化与假设，而且，不同的简化与假设，又可能导致不同的数学模型，例如[2]
是抛物型方程模型，而[3]则是椭圆方程模型.
假设：
（1）注射前大脑中的多巴胺含量可以忽略不计.
（2）大脑中多巴胺注射液经历着扩散与衰减的过程，且沿,,x y z 三个方向的扩散系数分别是常数，衰减使质量之减少与深度成正比.
（3）注射点在后排中央那个圆柱中心，即注射点的坐标000(,,)x y z 已知，注射量有医疗记录可查，是已知的.
（4）注射瞬间完成，可视为点源delta 函数. （5）取样也是瞬间完成，取样时间已知为1t =.
（6）样本区域与整个大脑相比可以忽略，样本组织远离脑之边界，不受大脑边界面的影响.
在以上假设之下，显然可以用本文前面讲过的思路来建模，于是得AMCM-90A 的数学模型为Cauchy 问题（21）（22），解的表达式为（24），且用三元回归分析来估出参数,,,a b c k ，于是可以求得任意位置任意时刻药物的深度.
如果所给数据认为是在平衡状态测得的，药物注射进脑后，从高深度处向低深度处扩散，与扩散同时，一部分药物进入脑细胞被吸收固定，扩散系数与吸收系数都是常数，但过一段时间，所有药物都被脑细胞所固定，达到了平衡态. 在这种假设下，[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 ，于是有
v
k v hv t
∂=∆-∂. （30）
又
w
hv t
∂=∂，即吸收速度与游离的浓度成正比，代入（30）得 ()v k w
w t h t t
∂∂∂=∆-
∂∂∂. （31） 对（31）关于t 从０到+∞积分得
000t t t k v w w h
+∞+∞+∞
====∆-. （32）
由于最后无游离药物，故(,,,)0v x y z +∞=，又开始时(0)t =无被吸收的药物，故(,,,0)0, (,,,0)0w x y z w x y z =∆=；平衡状态在t =+∞时达到，这时(,,)u x y z = (,,,)w x y z +∞，于是由（32）得
(,,,0)k
u u v x y z h
-∆+=， （33）
其中(,,,0)v x y z 是开始时的浓度分布，近似于注射点的点源脉冲函数. 把此注射点取为坐标
原点(0,0,0)，则(,,,0)(,,),v x y z L x y z L δ=是注射量，于是2
k h σ⎛⎫= ⎪⎝
⎭记
2(,,)u u L x y z σδ-∆+=， （34）
作付立叶变换得
22222222ˆˆ(),ˆ,
1()
s u u L L
u
s σξησξη+++==+++ 再作反变换得
u σ-=-， （35）
其中C 是可计算常数.
如果考虑各向不同性，设,,x y z 方向上扩散系数分别为222
,,a b c ，注射点在
000(,,)x y z ，则
222222000222()()()u u u a b c u L x x y y z z x y z δδδ⎛⎫
∂∂∂-+++=--- ⎪∂∂∂⎝⎭
， 于是解为
(,,)u x y z =
exp 1⎧⎪-⎨⎪⎩，（36） （36）中的D 可计算常数.
用前面类似的方法可以进行参数估计.
在建模过程中，点源函数的使用显然与实况有差别；尤其是认为扩散系数与吸收系数都是常数，对于人脑这种有复杂结构的区域，这种假设与实际不会完全符合；夜间与白天（睡与醒）对这些系数有无影响？脑中各点这些系数是否有变？除时间位置应考虑外，可能还与药液浓度有关. 如此看来，脑内药液分布的数学模型很可能不是常系数线性偏微分方程，而是函数系数的线性微分方程甚至是非线性偏微分方程. 这时，其解不再能用封闭公式来表达，求解过程会变得极为复杂，所以也可以考虑是否试用其他数学模型来解，例如在平衡态的假设下，用回归分析方法建立药液的模拟分布(,,)u f x y z =.
对一个实际问题，其数学模型未必唯一，各模型间孰优孰劣，没有一般的判别法，须经实践来检验.
参 考 文 献
[1]叶其孝，大学生数学建模竞赛辅导教材，湖南教育出版社，1993.
[2]Christopher, R. Malone, Gian Pauletto, James, I. Zoellick, Distribution of Dopamine in the Brain, The Journal of Under graduate Mathematics, and its Applications, vol. 12(1991), Special Issue: The 1991 Mathematical Contest in Modeling, pp. 211-223.
[3]孙晓东，荆秦，梁俊，脑中药物分布的数学模型，数学的实践与认识，1991年No. 4，63-69. [4]中国科学院数理统计组，常用数理统计方法，科学出版社，19784.。