扩散方程

2024-06-22

扩散方程（精选八篇）

扩散方程篇1

水是人类赖以生存的资源, 保护水资源就是保护人类自身。随着经济社会的发展, 河流的水质污染状况越来越严重。水体热污染、副营养化等状况已经严重威胁到生态平衡。由于在河流流域内的不同地区的主要污染物不同, 以及污染状况往往是由多种污染物造成, 因此对污染的治理的方法也相对复杂。同时河流湖泊对生产、生活的各个方面也发挥这重要的作用。例如:内河航运、旅游休闲、调节气候等。现今, 随着科技的迅速发展, 如何准确的判定出一个地区的主要污染物也越来越重要。

本文基于对污染物指标的归一化处理, 通过borad排序法对水质状况进行了评价。然后利用matlab结合反映扩散理论建立出数学模型。通过引用长江流域内的观测数据证明了模型的合理性。

1污染指标的选取与样本数据

河流的污染状况往往是由多种因素所决定的, 根据水质检测[1]的主要指标可以将水质状况分为六种不同的类型, 同时结合国家标准 (GB 3838-2002) [2]的规定, 对于地表水质的评价指标一共有24项指标, 但是其中对水质影响最大的一共有四项指标:PH值、溶解氧、高锰酸盐、氨氮。具体标准如表1所示。

其中第I、II、III类水为可饮用水, 溶解氧 (DO) 的指标与高锰酸盐 (CODMn) 和氨氮 (NH3-N) 呈反向变动。

2水质的综合评价模型

在进行水质的综合评价时需要对其余因素进行量纲化处理, 使得数据具有可比性。通过对表1的分析, 由于第I、II、III类水为可饮用水, 因此可以将第III类水的对于四种主要污染物的指标作为临界值进行比较。

(1) 对于高锰酸盐 (CODMn) 与氨氮 (NH3-N) 的处理。因为在第III类水中, 高锰酸盐 (CODMn) 对应的指数为6。同时假设在河流流域内第i个观测站第j个时刻时高锰酸盐 (CODMn) 的实际观测值为δ'1。对于高锰酸盐 (CODMn) 通过与第III类水质的标准进行比较去量纲, 可以得到δ1=δ'1/6。同理氨氮 (NH3-H) 即为δ2=δ'2/1.0。

(2) 对于溶解氧 (DO) 的处理。一般情况下污染值越大, 则污染越严重。但是对于溶解氧 (DO) , 则污染物的含量与污染状况则呈反比例变动的关系。因此若假设溶解氧的实际观测值为δ3", 根据反比例关系可以得到δ'3=1/δ3"。然后与标准值比较去量纲化得到δ3=5/δ'3。

(3) 对于PH值的处理。由表1可知, 对于非劣V类水PH值趋于6-9之间。同时根据实验分析, 通常认为当PH值越接近于7时河流被污染的可能性则越小。若假设PH的实际观测值为δ'4则可以对PH值进行如下处理:

(4) 四种指标的综合评价。根据不同的指标所对应的分数等级, 通过利用非模糊数学判断矩阵法对四种不同的观测物设置权重αk (其中k=1, 2, 3, 4) 。然后结合得分系数矩阵γkij进行加权求和, 则可以得到污染指标εik的数学模型即为:

其中:i——不同的观测点, j——不同的观测时间

模型的检验, 通过引用长江水质的观测数据, 对2003年到2005年间的数据进行分析。通过对高锰酸盐、氨氮、溶解氧、PH值的权重αk进行计算并且进行排序, 可以得到在长江流域17个观测点处28个月中的污染指标εik综合得分排序入表2所示。

根据表中数据可得:水质状况最好的是在湖北丹江口胡家领地区, 其次是江苏南京林山;水质状况最差的是在江西南昌滁槎其次是四川乐山岷江大桥和四川泸州沱江二桥。因此与当年的评测结果基本相符。

3污染源位置的确定

3.1基于扩散方程[4]的数学模型

根据参考文献, 可以知道扩散方程表征了流动系统的质量传递规律, 通过对其进行求解可以得到物质的浓度分布。所以根据扩散方程:

因此每个测量点的水流量P根据实际测出的污染物的浓度T, 可以得出污染物的量为P×T。若第一个观测站点的污染量忽略不计, 则可以得出每个测量点自己产生的污染物的量的数学模型为:

其中:Pi——表示第i个观测点的水流量, Ti——表示第i个观测点所观测的污染物的量。

3.2模型的检验

代入2003年到2005年间长江干流的数据, 取降解系数λ=0.2利用matlab可以得到CODMn和NH3-N两种污染物的量如表3所示。

综上所诉, 高锰酸盐 (CODMn) 的主要污染湖南岳阳、江苏南京、江西九江等地区;而氨氮 (NH3-N) 的主要污染地在湖南岳阳、江西九江等地区。因此在治理污染时, 应该对不同的地区采取的不同的治污措施以优化环境的治理效率。

4结论

本文通过对水质资源的污染研究状况给出了测评方法, 并且带入长江流域内的观测数据进行验证。

在四种主要的污染指标影响下, 水质状况最好的地区是在会被丹江口胡家领地区, 其次是江苏南京林山。而水质最差的地区是江西南昌滁槎其次是四川乐山岷江大桥和泸州沱江二桥。基于扩散方程可知。湖南岳阳、江苏南京的主要污染物为高锰酸盐 (CODMn) 。而在湖南岳阳、江西九江地区的主要污染物也包括氨氮 (NH3-N) 。

通过对长江流域观测数据的分析, 得到的结论是可靠的, 与实际发生的状况差距并不大。因此, 模型具有合理性。

参考文献

[1]蒯圣龙主编.水污染与水质监测[M].合肥市:合肥工业大学出版社, 2010 (09) .

[2]安乡县水利局, 中华人民共和国国家标准GB3838-2002[S].http://sl.anxiang.gov.cn/art/2013/5/29/art_2058_42608.html, 2013, 5 (29) .

反应扩散方程解的渐近性态篇2

反应扩散方程解的渐近性态

研究了一类反应扩散方程初始边值问题.利用微分不等式理论得到了问题解的.渐近性态.

作者：莫嘉琪唐荣荣 Mo Jiaqi Tang Rongrong 作者单位：安徽师范大学数学系,芜湖,241000;湖州师范学院数学系,湖州,313000;上海高校计算科学院E-研究院上海交通大学研究所,上海,240 刊名：数学物理学报 ISTIC PKU英文刊名：ACTA MATHEMATICA SCIENTIA 年，卷(期)： 27(2) 分类号：O175.29 关键词：反应扩散方程渐近性态微分不等式

扩散方程篇3

本文讨论如下非自制反应扩散方程:

非自治指数吸引子的存在性, 其中ΩRn是具有光滑边界的有界区域.f是Lipschitz连续函数, 并且满足

和

g (x, t) ∈Lb2 (R, L2 (Ω) ) (i.e., 平移有界) , 并且满足

其中M>0, Ci>0, i=1, 2, 3, 4.

2 预备知识

这部分, 我们给出相应的基本概念和非自治指数吸引子的存在性定理.

(ⅰ) 对t≥τ, τ∈R, M (t) 是正不变的:U (t, τ) MM;

定义2[7]设E和E1是Banach空间, 且E1紧嵌入到E.设B是E1中的一有界子集, 对于给定的正常数ε∈[0, 1]和δ, K>0, 定义非线性算子S:E→E的一个类Sδ, ε, K (B) 如下:

其中σ-邻域赋予E1中的拓扑结构.

(ⅱ) S可分解为如下形式

其中S0和S1满足下列情形

和

结合引理3, 对于任意的τ∈R, Ug映H10 (Ω) 中的有界集为有界集, 即得如下推论.

为了后边的证明, 我们需要把Ug (t, τ) uτ=u (t) 分解为和的形式

其中D (t, τ) uτ=v (t) 和Kg (t, τ) uτ=w (t) 分别为如下两个方程的解

和

引理5[6]对于任意τ∈R, 方程 (1) 的解满足如下估计:存在常数k0, 使得对任意t≥τ,

其中Q1 (·) 是在[0, ∞) 上的增函数, 并且k0依赖于λ1, 并且它们都依赖于τ.

在文[7]中, 作者首先估计了离散过程指数吸引子存在性的一些抽象结果, 并且在文[7]中的第4部分给出了由一个离散的指数吸引子转换为连续的指数吸引子的一些技巧.为了方便引用, 我们总结如下:

引理6 (指数吸引子的存在性) 设E和E1是两个Banach空间, 并使得E1紧嵌入到E, B是E1中的一个有界集U (t, τ) (τ∈R, t≥τ) 是在Banach空间E中的过程.进一步假设U (t, τ) 满足

和

则存在一个指数吸引子τ→εU (τ) Qδ (B) , τ∈R.且满足如下性质

(ⅰ) 吸引子εU (τ) B, (τ∈R) 和它的分形维数是有限的.i.e.,

其中常数C1不依赖于τ;

(ⅱ) 关于U的集族εU是正不变的, i.e.,

(ⅲ) 集族有如下 (一致) 指数吸引子的形式

其中常数C2和α不依赖于t, τ;

(ⅳ) 函数t→εU (t) 是一致H9lder连续的

其中常数C3和k不依赖于t, s.

3 主要结果及证明

证明结论之前, 我们需要一下引理.

引理7设u (t) 是问题 (1) 所对应初值uτ∈V的一个解, 则对于任意的ε>0和τ∈R, 存在正常数Cε和Kε使得

其中u1 (t) , w1 (t) , 满足如下估计

其中常数Cε和Kε不依赖于τ.

其定理的证明类似于文[9]中的引理2.5.这里证明从略.

证明首先, 由引理3, 我们可知存在一个仅依赖于B的常数TB, 使得

用-Δu (t) 和 (1) 式做内积, 我们得到

其中

结合 (2) , (3) 有

将 (20) , (21) 带入 (19) , 整理得

根据Gronwall不等式, 并且在[τ+TB, t]上积分可得到 (16) .

为了方便期间, 我们记

引理9设引理6成立, 对任意固定的时间T>0, 过程Ug (t, τ) 在下面的情形下是H9lder连续的:存在正常数CT, B和c′, 使得对于任意zi∈Qδ (B) (i=1, 2, s∈[0, T]是H9lder连续的

和

证明首先证明 (22) .从引理3和方程 (1) 可知u (t) =Ug (t, τ) uτ (t≥τ) 满足:对于任意的uτ∈Q1 (B) ,

其中MB是和Q1 (B) 有关的常数.对于任意s∈[0, T]和z2∈Q1 (B) , 我们有

现在对于 (23) 式, 根据文[7]中的证明, 我们有:对于t≥T和s∈[0, T],

结合 (22) 和 (24) 式, 我们能得到 (28) .

引理10[8]设 (2—4) 成立, 并且g∈Lb2 (R, L2 (Ω) ) .则对于任意的τ∈R和初始值uτ∈H, 则问题 (1) 有唯一解.进一步, 设{uτi, gi} (uτi∈H, gi∈Lb2 (R, L2 (Ω) ) , i=1, 2) 是两个初始条件, 并且用ui表示问题 (1) 的解, 则如下估计成立:τ≤t≤T≤T+τ, 成立

引理11在引理6的假设下, 存在一个T=T (B) , 使得过程Ug (t, τ) (t≥τ) , 满足如下性质

(ⅰ) 对于任意的τ∈R和t≥0

(ⅱ) 对于任意的z1, z2∈Q1 (B) 和τ∈R:Ug (T+τ, τ) 可分解为如下形式

其中S0T和S1T满足下列条件

和

证明根据引理3和引理9, 结合文[10]中的引理5.6, 对于每一个初始值uτ∈Q1 (B) , 我们有如下的分解:用S (t, τ) uτ表示问题 (1) 的其次方程的解, 并设

对于初始值uτi∈Q1 (B) 和对应的解ui (t) =Ug (t, τ) uτi (i=1, 2) , 设uτ=uτ1-uτ2, 我们分解u1 (t) -u2 (t) 为如下和的形式

下面我们分3步证明.

第1步:从 (29) 我们能得到

其中常数C1和k1仅依赖于特征值λ1.因此, 当T1=t-τ足够大时, 我们有

对于非线性项, 我们有

其中θ∈[0, 1].

第3步:由于引理8, 对于D (A) 中的有界集Q1 (B) , 存在一个仅依赖于B的常数T2, 使得对于任意的τ∈R成立

定理12 (非自治指数吸引子) 设ΩRn是有界光滑区域, 并且f满足 (2—3) , 外力项g (x, t) ∈Lb2 (R, L2 (Ω) ) , 则对应的过程Ug (t, τ) 在V中有一个非自治的指数吸引子{εg (t) }t∈R, 满足如下性质:

(ⅰ) εg (t) 在H中是紧的, 对于任意t∈R, 它的分形维数是有限的.i.e., dimF (εg (t) , H) ≤C1, t∈R, 其中常数C1不依赖于t;

证明根据引理9和引理10, 我们知对于集合B和满足 (4) 的每个外力项g (x, t) , 存在一个非自治的指数吸引子{εg (τ) }t∈R满足 (12—15) .则我们可判定{εg (τ) }t∈R满足引理6中的所有条件.证毕.

参考文献

[1]A.V.Babin, M.I.Vishik.Attractors of Evolution equations[M].North-Holland:Amsterdam, 1992.

[2]V.V.Chepyzhov, M.I.Vishik.Attractors for equations of mathematical physics[M].Providence:American Mathematical Society, 2002.

[3]B.C.Robinson.Infinite-dimensional dynamical systems:An Introduction to Dissipative Parabolic PDEs and the Theory of Global Attractors[M].Cambridge University Press, Cambridge, 2001.

[4]A.Eden, C.Foias, B.Nicolaenko, R.Temam.Exponential attractors for dissipative evolution equations[M].John-Wiley, New York, 1994.

[5]A.Miranville.Exponential attractors for nonautonomous evolution equations[J].Appl.Math.Letters, 1998, 11 (2) :19-22.

[6]M.Efendiev, A.Miranville, and S.V.Zelik.Infinite dimensional exponential attractors for a non-autonomous reaction-diffusion system[J].Math.Nachr, 2003, 248 (2) :72-96.

[7]M.Efendiev, A.Miranville, and S.V.Zelik.Exponential attractors and finite-dimensional reduction for non-autonomous dynamical systems[J].Proc.Royal.Soc.Edinburgh, 2005, 135A:703-730.

[8]R.Czaja, M.Efendiev.Pullback exponential attractors for non-autonomous equations Part II:Applications to reaction-diffusion systems[J].J.Math.Anal.Appl., 2011, 381:766-780.

[9]S.V.Zelik.Asymptotic regularity of solutions of a nonautonomous damped wave equation with a critical growth exponent[J].Comm.Pure Appl.Anal., 2004, 3:921-934.

扩散方程篇4

一类时滞反应扩散方程的稳态解和行波解

本文讨论了一类造血生物模型在Dirichlet边值条件下稳态解的全局吸引性.并利用上,下解技术和单调迭代方法讨论了行波解的存在性.

作者：赵书芬黄建华 ZHAO Shu-fen HUANG Jian-hua 作者单位：国防科技大学,理学院数学与系统科学系,湖南,长沙,410073 刊名：生物数学学报 ISTIC PKU英文刊名：JOURNAL OF BIOMATHEMATICS 年，卷(期)：2008 23(2) 分类号：O175 关键词：上、下解稳态解全局吸引性行波解

扩散方程篇5

图像分割就是把图像分解成具有某种特性的若干区域并提取出感兴趣目标的过程。伴随图像分析的发展, 图像分割方法层出不穷, 从简单的阈值技术到基于数学模型的最优方法, 从二维图像到三维图像[1], 并在实际中有了大量的应用。虽然图像分割方法已经有了很大的发展, 但由于它的复杂性, 仍有很多问题没有很好的解决, 如图像分割具有较强的针对性, 某种图像分割方法只能解决一些特定的应用问题;有时为了改善图像分割性能, 甚至要组合使用几种分割方法等。

自1987年Kass等人[2]提出主动轮廓模型以来, 基于曲线演化 (curve evolution) 的形变模型已被广泛应用于图像分割。形变模型的图像分割方法具有能够有效结合图像本身的低层次视觉属性与待分割目标先验知识的灵活开放的框架, 可获得分割区域的完整表达。目前常见的轮廓演化模型有两种, 即基于参数的模型和基于几何特性的模型[9]。参数形变轮廓模型 (snake模型) 是一种能量函数最小化的形变轮廓曲线, 常常能得到很好的整体效果。但由于模型本身的缺陷, 使得snake模型对初始位置过于敏感, 不能处理拓扑结构改变并且易于陷入局部极值, 使其稳定性难以满足复杂图像分割的要求。基于几何活动轮廓线 (geometric active contours) 的水平集 ( level set) 图像分割方法是处理封闭运动界面随时间将移动的界面作为零水平集 (zero level set) 嵌入高一维的水平集函数中。这样, 由闭超曲面的演化方程可以得到水平集函数的演化方程, 而嵌入的闭超曲面总是其零水平集, 最终只要确定零水平集就可以确定移动界面演化的效果。水平集方法在一定程度上克服了snake模型的缺点, 如对初始轮廓的选择没有特殊的要求, 就可以很好地处理拓扑结构改变等。且水平集方法主要是基于图像的梯度分割图像, 能够直接产生闭合的曲线或曲面, 并对噪声和伪边界有很强的鲁棒性;但是它对初始边界位置十分敏感, 有时还要求人工选择合适的参数。另外水平集方法巨大的计算量是一个无法回避的问题, 当初始轮廓线不断演化、变长, 算法的计算量将是无法忍受的。

考虑到水平集方法巨大的计算量, 本文提出了基于格子波尔兹曼的对流扩散方程活动轮廓分割模型。格子波尔兹曼模型 (Lattice Boltzmann Model) 是一种新的求解偏微分方程的数值工具。它是由陈十一和钱跃兹在1991～1992年间提出的[3,4], 可以用于求解Navier-Stokes等偏微分方程。LBM是一种新的数值工具, 与传统数值方法不同的是, 它不涉及对偏微分方程的直接求解, 而是采用自底向上的方法, 通过模拟粒子的移动间接的实现对偏微分方程的模拟。利用该模型求解对流扩散方程, 取代水平集方程实现图像分割, 可极大地提高算法的速度, 并能够获得一个闭合的分割曲线。

1 LBM 原理

格子波尔兹曼模型 (Lattice Boltzmann Model) 是一种新的求解偏微分方程的数值工具。目前已经被广泛的用来模拟动态流体的流动[5]等方面。格子波尔兹曼模型起源于模拟流体运动规律的格子气元胞自动机, 并有如下改进:

a) 元胞空间从三角形网格改进为方形网格, 简化了表示元胞空间的数据结构, 更易于编程实现;

b) 元胞的状态空间从离散的空间改进为连续的空间, 消除了格子气元胞自动机的布尔性质;

c) 碰撞项改进为具有单一松弛时间的BGK (Bhatnagar-Gross-Krook) 估计[6], 从而显著地降低了计算量。

格子波尔兹曼自动机的一种常用的元胞空间如图1所示。每个元胞与其东、西、南、北以及东南、西南、东北、西北共8个方向的邻居相连。每个粒子按这8个方向 (i=1, 2, …, 8) 向元胞邻居移动或者停留在原位置 (i=0) 静止不动, 并在各个网格结点内发生碰撞。

移动阶段可以表示为:

$f_{i}^{*} (r + Δ h \cdot C_{i}, t) = f_{i} (r, t), 1 \leq i \leq 8 (1)$

碰撞阶段可以表示为:

$\begin{array}{l} f_{i} (r, t + Δ t) = f_{i}^{*} (r, t) + \\ \frac{1}{ζ} [f_{i}^{(0)} (r, t) - f_{i}^{*} (r, t)], 1 \leq i \leq 8 (2) \end{array}$

综合式 (1) 和式 (2) , 格子波尔兹曼自动机的演化规则可以表示为:

$\begin{array}{l} f_{i} (r + Δ h \cdot C_{i}, t + Δ t) = f_{i} (r, t) + \\ \frac{1}{ζ} [f_{i}^{(0)} (r, t) - f_{i} (r, t)] (3) \end{array}$

其中:fi (r, t) 为粒子密度分布函数;f $_{i}^{(0)}$ (r, t) 为局部平衡分布函数, 表示元胞内粒子由于碰撞达到平衡态时的状态;ζ为松弛时间, 表示元胞内粒子密度趋向于平衡分布函数的时间, 表征模型趋向于平衡态的快慢;Δh为移动距离;Ci为方向矢量;Δh·Ci表示某个方向的移动距离;Δt为碰撞时间。相对于传统的数值解法, LBM具有稳定性好, 编程简单, 计算效率高的特点, 另外LBM是一个天然的离散系统, 因此是应用于图像处理的极好工具。

2 对流扩散方程模型

对流扩散方程是另一种重要的偏微分方程, 它可以描述质量、热量的运输过程以及反应扩散过程等众多物理现象。对流扩散方程可以表示为如下形式:

$\frac{\partial}{\partial t} ρ = a \nabla^{2} ρ + b | \nabla ρ | (4)$

该方程来源于流体力学动力学模型。其物理意义是, 流体在流动过程中, 各点的粒子密度随时间的变化不仅仅受到因粒子密度不均匀而引起的扩散现象的影响, 还要加上因流体流动而引起的对流项的影响。由对流项引起的粒子密度的变化直接与粒子密度的梯度|∇ρ|相关。基于这样的物理意义, 为了实现利用元胞自动机求解对流扩散方程的目的, 一个简单的方法就是在各向异性扩散模型的演化方程中, 人为的加上一个对流项Δtλ|∇ρ|, 其中λ为一个常数。这样元胞自动机的演化方程可以写为:

$\begin{array}{l} f_{i} (r + Δ h C_{i}, t + Δ t) = g_{i} (r) {f_{i} (r, t) + \\ \frac{1}{ζ} [f_{i}^{(0)} (r, t) - f_{i} (r, t)] + Δ t λ γ} + [1 - g_{i} (r)] \times \\ {f_{i} (r + Δ h C_{i}, t) + \frac{1}{ζ} [f_{i}^{(0)} (r + Δ h C_{i}, t) - \\ f_{i} (r + Δ h C_{i}, t)]} (5) \end{array}$

其中, gi (r) 为某个方向粒子通过概率,

$γ = | \nabla ρ | + \frac{v}{9 λ} (C_{i} \cdot \nabla) ρ$

若定义 $f_{i}^{0} (r) = \frac{ρ (r)}{9}$ 。将式 (5) 左边泰勒展开, 代入粒子分布密度函数, 等式两边同取O (ε2) , 再分别对i求和, 便可以得到该模型的宏观方程:

$\begin{array}{l} \frac{\partial}{\partial t} ρ = - \frac{λ}{9} | \nabla ρ | \sum_{i} g_{i} (r) + \\ \frac{1}{2} v^{2} Δ t \sum_{i} g_{i} (r) (C_{i} \cdot \nabla)^{2} f_{i}^{(0)} (6) \end{array}$

其中取局部平衡函数:

即得:

$\frac{\partial}{\partial t} ρ = \frac{Μ Δ h^{2}}{2 Δ t} \sum_{i} g_{i} (r) (C_{i} \cdot \nabla)^{2} ρ - \frac{λ}{9} | \nabla ρ | \sum_{i} g_{i} (r) (7)$

从式 (7) 可以看出, 这里提出的元胞自动机模型的宏观方程与对流扩散方程是一致的。

图像分割中经常用到的水平集方程, 从实质上讲就是一种对流扩散方程。水平集方程:

$\frac{\partial}{\partial_{t}} ρ = b k | \nabla ϕ | - V \cdot \nabla ϕ (8)$

3 实验分析

式 (8) 的水平集方程, 它可替换为对流扩散方程[7]。但问题是, 无论是对流扩散方程还是水平集方程用传统的数值方法都很难求解, 在式 (8) 中, k为距离场的曲率, 当我们利用式 (8) 实现图像分割时, k的值会随着距离场的变化而变化, 因此每一步迭代中都必须更新k, 这大大增加了算法的计算量, 尤其是当初始轮廓线变得很长的时候, 这种计算量有时是不可忍受的。

根据以上分析, 这里我们利用LBM求解对流扩散方程, 利用对流扩散方程取代水平集方程实现图像分割。实验表明我们设计的算法可以得到闭合的分割曲线, 且无需在每次迭代过程中求距离场的曲率, 另外, 由于LBM的并行计算性能, 可以极大地提高算法的速度。

3.1 基于LBM对流扩散模型的图像分割算法

该算法与水平集算法相似。首先确定初始轮廓, 依据这个初始轮廓, 建立一个距离函数。利用LBM方法模拟该距离函数的演化, 从而实现对图像的分割。详细描述如下:

1) 首先必须对模型初始化。建立热量场 (初始轮廓) , 将图像中的像素点看作是LBM模型中的元胞, 将原始图像看作是元胞自动机模型的初始值, 粒子密度分布函数和局部平衡分布函数初始值设置如下:

$f_{i} (r, 0) = f_{i}^{(0)} (r, 0) = {\begin{cases} Μ \frac{ρ (r, 0)}{9}, i \neq 0 \\ 1 - \frac{8 Μ ρ (r, 0)}{9}, i = 0 \end{cases} (9)$

2) 迭代过程。粒子在热量场中进行扩散, 在各个元胞内部发生碰撞, 从而对元胞的状态产生影响, 即

$\begin{array}{l} f_{i} (r + Δ h C_{i}, (n + 1) Δ t) = g_{i} (r) {f_{i}^{*} (r, n Δ t) + \\ \frac{1}{ζ} [f_{i}^{(0)} (r, t) - f_{i}^{*} (r, n Δ t)]} + [1 - g_{i} (r)] \cdot \\ {f_{i}^{*} (r + Δ h C_{i}, n Δ t) + \frac{1}{ζ} [f_{i}^{(0)} (r + Δ h C_{i}, n Δ t) - \\ f_{i}^{*} (r + Δ h C_{i}, n Δ t)]} (10) \end{array}$

更新元胞状态 $ρ (r, (n + 1) Δ t) = \sum_{i = 0}^{\infty} f_{i} (r, (n + 1) Δ t)$ , 通过找到满足ρ (r, t) =0的点更新分割轮廓线边缘。

每次碰撞后必须重新计算局部平衡分布函数和元胞自动机模型。

$f_{i}^{(0)} (r, n Δ t) = {\begin{cases} Μ \frac{ρ (r, n Δ t)}{9}, i \neq 0 \\ 1 - \frac{8 Μ ρ (r, n Δ t)}{9}, i = 0 \end{cases} (11)$

如果达到最大迭代次数N, 则退出迭代。

3) 输出分割目标边缘。

3.2 参数选择

算法中, 距离函数定义为网格中各点到初始轮廓线Γ的距离。实验中, 当网格处于分割轮廓线Γ内部时, 距离函数为正;当处于分割轮廓线Γ外部时, 距离函数为负。由于水平集方法中, 分割轮廓线的演化速度是定义在分割轮廓线Γ上的 (从LBM模型的宏观方程可以看到, 该演化速度在LBM演化方程中对应着半透膜粒子通过概率gi (r) ) , 因此必须将演化速度gi (r) 向整个网格空间扩展。扩展的方法是, 如网格中任意一点r′到分割轮廓线Γ上一点r的距离是r′到分割轮廓线Γ的最短距离, 那么定义r′处的gi (r′) 等于gi (r) 。

以上算法中的另一个关键点是, 距离函数的重新初始化问题。由于距离函数始终满足|ρ|=1, 因此我们采用如下的迭代方法重新初始化距离函数[8]:

$ρ_{t} = 1 - | \nabla ρ_{t - 1} | + ρ_{t - 1} (12)$

从式 (12) 可以看出, 当|∇ρt-1|大于1或接近于1的时候, ρt不断减小, 而当|∇ρt-1|小于1的时候, ρt不断增大。只有当|∇ρt-1|等于1的时候, 迭代达到稳定状态。而且这个迭代过程的收敛非常迅速, 不需要花费很长的计算时间。

在试验中, 我们实际选取M=0.001, Δt=0.01。

3.3 评价与分析

为了能够对本文提出的算法进行一个定性的评价, 我们将该算法分别应用到了CT图像、DSA图像以及MRI图像的分割, 如图2～图4所示。

图2 (a) 是原图像, 该图像为128×128的8 bit灰度图。图2 (b) 为利用本算法的分割结果。比较图2 (a) 和图2 (b) 可以看到, 图像中的目标区域 (肝脏) 在图2 (b) 中已经被很好分割出来了。图3给出了利用本算法对数字减影血管造影图像中的血管进行分割的结果, 而图4则给出了利用本算法对脑部核磁共振图像中大脑脑干的分割结果。

从实验结果看, 本文中提出的基于LBM的图像分割算法是可行和有效的。

为评价算法的计算速度, 下面仅给出算法与水平集方法和窄带水平集方法的计算时间。我们分别利用这3种方法分割图5 (a) 中的黑色圆形区域, 而表1则给出了它们的计算时间。整个实验基于WindowsXP平台, 用Matlab编程实现。图5 (a) 为82×82的二值图像, 图5 (b) 为初始分割轮廓线。从表1中可以看到, 利用本文提出的算法对图5 (a) 进行图像分割时共耗时561 s, 而窄带水平集方法和经典的水平集方法的耗时都要大于这个时间, 这说明利用我们的算法, 图像分割的速度可以得到极大的提高。

为了检测LBM模型计算的精确度, 我们利用基于LBM的对流扩散方程求解边界-初始值问题。图6给出了利用LBM的数值解和解析解在时间t=70时的水平线y=70上的值。从图中可以看出, 在边界地区, 由于边界处理的原因数值解与解析解的误差较大, 而在其它的大部分区域数值解与解析解几乎完全重合。这说明, 提出的基于LBM的对流扩散模型能够精确的模拟对流扩散方程。

4 结束语

文中为了进行图像分割, 首先将水平集方程替换为对流扩散方程, 但是这两个方程都十分的难以求解。研究中尝试了利用LBM求解来对流扩散方程。基于这个对流扩散模型, 我们实现了基于LBM的水平集图像分割。实验证明该方法能够精确地求解对流扩散方程, 且该方法是稳定的。图像分割的实验结果表明, 该算法能够很好的实现对包括CT图像、核磁共振图像、数字减影血管造影图像在内的各种图像的图像分割, 并大大提高了处理速度。更为重要的是, 该算法可以做成并行图像处理系统, 这更会使图像分割的速度得到一个飞跃。

参考文献

[1]PARKER B J, FENG D J.Graph-based Mumford-Shah seg-mentation of dynamic PET with application to input functionestimation[J].IEEE Transactions on Nuclear Science, 2005, 52 (1) :79-89.

[2]KASSM, WITKINA, TERZOPOULOS D.Snakes:active cont-ourmodels[J].International Journal of Computer Vision, 1987, 1 (4) :321-331.

[3]CHENS Y, CHENHD, MARTINEZ D, et al.LatticeBoltzmannmodel for sim ulation of magnetohydrod ynmnics[J].PhysicsReviewLetter, 1991, 67 (27) :3776-3379.

[4]QIAN Y H, HUMIERES D D, LALLEMAND P.Lattice BGKmodelsforNavier-Stokes equationl[J].Europhysics Letters, 1992, 17 (6) :479-484.

[5]KOELMAN J M V A.A simple Lattice Boltamann Scheme fornavier-stokes fluid flow[J].Europhysics Letters, 1991, 15 (6) :603-607.

[6]GHATHAGAR P L, GROSS E P, KROOK M.A model forcollision processes in gases.I.Small amplitude processes incharged and neutral one-component systems[J].Physical Re-view, 1954, 94 (3) :511-525.

[7]于志伟, 陶波.图像重建的自适应随机元胞自动机算法[J].计算机学报, 1997, 20 (10) :943-948.

[8]KARAFYLLIDIS I.ANDREADIS P.TZIONAS P, et al.Acellular automaton for the determination of the mean velocityof moving objects and its VLSI implementation[J].PatternRecognition, 1995, 29 (4) :689-699.

[9]任继军, 何明一.一种整合单通道各向异性扩散信息的水平集图像分割方法[J].计算机应用研究, 2008, 25 (1) :300-301.

[10]吴晓红, 罗代升, 王正勇, 等.基于水平集和迭代自组织算法的图像分割[J].计算机工程与应用, 2007, 43 (3) :16-18.

[11]陈樟, 苏伟, 万敏.微通道气体流动的格子波尔兹曼模[J].功能材料与器件学报, 2008, 14 (1) :51-54.

扩散方程篇6

扩散方程是一类重要的抛物型偏微分方程,涉及到传热,传质等传导现象,在诸如能源,机械,动力等诸多领域,都会要用到与时间相关的扩散方程来描述。由于边界元法具有降维,数据准备工作量省,精度高等优点,故采用边界元方法求解此类方程是一种可选的方法。从上世纪60年代将边界元法应用于数值计算以来,人们对抛物型方程的边界元解法的实施主要依据直接边界元解法。Onishi,Zhu等人曾研究过采用带时间变量的基本解,利用间接边界积分方程及其等价的Galerkin变分形式求解二维扩散方程,但他们的讨论偏重于理论,未实施于数值计算。Galerkin方法是基于变分原理基础上的一种把微分方程或积分方程转化为等价的变分方程,然后离散变分方程求原方程数值解的方法。Galerkin方法比配点法更便于进行理论分析,比如解的存在唯一性,近似解的收敛性及误差估计等。由于实施Galerkin边界元解法,有求解四重奇异积分及超越函数积分的麻烦,少有实施用于数值计算的报道。

针对二维扩散方程的Dirichlet初边值问题,具体实施基于单层位势的间接边界积分方程的Galerkin边界元解法。在采用常单元的情况下,推导了计算中涉及到的关于时间和空间的四重奇异积分的所有公式,并较好地处理了超越函数。最后采用Fortran90编制的程序进行了数值试验,说明方法的可行性和有效性

1 二维扩散方程的间接边界积分方程的等价变分公式

考虑如下二维扩散方程的初边值问题:

其中Ψ是R2中的单连通开区域,Γ=Ψ是Ψ的边界并充分光滑,Ψ=R2Ψ是Ψ的闭包在R2中的补域,I是时间变量t的取值范围,I∶0

式中,H(t-τ)是海维赛德(Heaviside)函数,表示当t<τ时u*恒为0。方程(1)的求解可通过格林公式归化为求解它的间接边界积分方程中的未知量q(ξ,τ),如下式:

(2)式有如下的等价变分形式,对u∈L2(Ψ),g∈L2(0,T;H1⒊(Γ)),求q∈L2(0,T;H-1⒊(Γ)),使得p∈L2(0,T;H-1⒊(Γ)),满足a(q,p)=f(p),其中,

解的存在唯一性可以用Lax-Milgram定理证明。

把边界Γ离散为N个单元Γi(i=1,2,…,N),把区域Ψ分成M个单元em(m=1,2,…,M),把时间区间(0,T)等分为K个时间步,Δt=T K,tk=kΔt(k=1,2,…,K)。本文采用常单元进行计算,即假设q,p在每个边界单元和每个时间步上都是常数。

变分问题的离散形式可表示为下面的式子,

其中:

2 积分计算

首先解析计算(6)式中对时间t,τ的积分,即计算,

(1)计算内层积分

(2)计算外层积分

所以,

用同样的方法可以计算得到,

其中,Ei(·)是指数积分函数,可以展开成如下的级数形式:

其中,c0是欧拉常数,其近似值c0≈0.5772156649。

容易看出,当j=i-1,i,i+1时,(9),(10)是奇异积分,分两种情况来处理。

(1)当i=j时,Gijkk的奇异部分可以解析计算,不含奇异性的部分用二重高斯积分计算。事实上,因为Γi是直线段,可以设r=x-ξ=α-βl,那么,dsxdsξ=l2(α-β)2,其中,l是Γi的长度,则有,

(2)当i≠j时,对计算中对含有奇异性的部分第一重采用解析积分,第二重采用高斯积分,而对不含奇异性的部分用二重高斯积分计算。下面是解析积分的推导。

设x是积分源点,积分线段是从A到B的有向线段,沿线段的单位方向是的外法向量。x到的距离为d,R是x到ξ上任一点的向量,向量R1,R2和R在m上的投影为S1,S2和S,

则Gijkk的奇异项经过一重解析积分得,

同理,Gijkl所有奇异项也可解析计算,而剩下的积分全部采用高斯积分计算。Fik的计算同样可以先对时间解析积分,然后使用高斯积分。其中,

当Gijkk,Gijkl,Fik全部计算出来之后,就可以建立线性代数方程组,解出所有qjl,然后通过qjl求得任意点x在任意tk时刻的函数值u(x,tk)。

3 指数积分函数的处理

从前面的积分计算中看到,指数积分函数Ei(x)在计算中多次出现,由于Ei(x)中含有一个无穷级数项,在实际数值计算中,无法得到它的精确值,只能求一个近似值来代替它参与运算。这个近似值的精度大小将会很大程度的影响最终结果的精度。

(1)由其形式(11)式可知,当x较小时,Ei(x)收敛速度是比较快的,计算时可设定一个误差限,当通项小于误差限时可停止计算。

(2)当x较大时,Ei(x)的收敛速度非常慢,如果仍然采用上面的方法来处理,n值需要取得比较大,循环的次数会大量增加,从而导致程序效率低下。有文献给出了x较大时Ei(x)的级数展开式:

事实上,(13)式是通过分部积分推导出来的。对作N次分部积分,结果有,

它的截断项是与一般的级数不同,当x比较大时,上式收敛较快,如n=5,余差为

当x>10时,便有可见已经满足计算要求。

4 数值实验

单位半径的圆域问题,具有零初始条件和如图2所示的与时间有关的边界条件,参数c=5。采用本文的方法计算出r=0和r=0.8的结果,和文献[8]中的解析解比较,结果如图3所示,可以看出,数值结果具有相当高的精度。由此可见,本文中处理四重奇异积分和超越函数的方法是可行的,有效的。

摘要：对二维热传导方程的Dirichlet初边值问题,采用带时间变量的基本解,利用基于单层位势的间接边界积分方程及其等价的Galerkin变分形式求解,该方法涉及到与时空相关的四重奇异积分的计算。在采用常单元离散的情况下,推导了具体实施数值计算所需的积分公式,完成了数值实验,验证了该方法的有效性和可行性。

关键词：扩散方程,Galerkin边界元法,间接边界积分方程,奇异积分

参考文献

[1]Zhu J.AnIndirect Boundary Element Methodinthe Solutionof the Diffu-sion Equation[C].Boundary Element VIII Conference,1986:707-714.

[2]ShawR P.An integral equation approach to diffusion[J].Int.J.HeatMass Transfer,17.1974:693-699.

[3]Brebbia C A,Walker S.Boundary Element techniques in engineering[M].Newnes-Butterworths,London.1980.

[4]Costabel M.BoundaryIntegral Operatorsfor the Heat Equation[C].In-tegral Equations and Operator Theory,1990:498-552.

[5]Onishi K.Galerkin method for boundary integral equations in transientheat conduction[C]//C A Brebbia.WL Wendland.G Kuhn.editors.Boundary Elements IX,1987:231-248.

[6]Attaway D C.The BEMfor the diffusion equation:a feasibility study[J]//Morino L,Piva Reds.Boundaryintegral equation methods,theoryand application.Rome,1990:75-81.

[7]祝家麟.椭圆边值问题的边界元分析[M].北京:科学出版社,1991.

[8]C.A布瑞比亚,J CF泰勒斯,L C诺贝尔.边界元的理论和工程应用[M].龙述尧,刘腾喜,蔡松柏,译.国防工业出版社,1988.

扩散方程篇7

考虑如下椭圆方程:

式 (1) 中, Ω∈Rn为有界多角形区域, 在区域Ω上a= (a1, a2, …, an) 满足diva=0, ε为常数, 并且有∈≪|a|·diamΩ, f∈L2Ω。

该方程不仅在流体力学中有着重要的意义, 在石油开采, 地下水渗流等工程问题上也被广泛的应用着, 因此有着重要的现实意义, 然而由于ε≪|a|·diamΩ, 既强对流占优, 若用标准的Galerkin方法求解, 会发现因边界层引发的, 近似解的波动将会影响到整个区域, 而原方程的精确解并不存在这种波动。泡函数法的基本思想便是, 在标准Galerkin方法的基础上, 通过在每个小单元T上引入T内非0, ∂T和T外恒为0的泡函数, 来增加离散空间Vh基的个数, 从而解决近似解的震荡问题[1,2,3]。鉴于泡函数法在解决椭圆型强对流扩散问题时所显示出的优越性, 本文通过构造适当的泡函数空间把泡函数法应用于相应的抛物问题, 从而将椭圆问题的泡函数方法推广到抛物问题, 并给出了收敛性分析和数值算例。

1 问题的提出及算法的实现

考虑如下强对流—扩散方程的初边值问题

式 (2) 中, Ω∈Rn为有界多解形区域, ut为u关于时间t的导数, 在区域Ω上a= (a1, a2, …, an) 满足

diυa=0, ε为常数, 并且有ε≪|a|·diamΩ, g (·, t) ∈

L2 (Ω) 。

在时间维上取时间步长为Δt, 用

$\frac{u (t_{n}) - u (t_{n} - 1)}{Δ t}$ 来逼近 $u_{t} (t_{n}), 1 \leq n \leq \frac{Τ_{E}}{Δ t}$ 。则式 (2) 的 (a) 式在t=tn时可转化为

$\frac{u (t_{n}) - u (t_{n - 1})}{Δ t} - ε Δ u (t_{n}) + a \nabla u (t_{n}) = g (t_{n}) + e^{n}$ (3)

式 (3) 中 $e^{n} = \frac{u (t_{n - 1}) - u (t_{n} - 1)}{Δ t} - u_{t}$ 。

令 $σ = \frac{1}{Δ t}, f (t_{n}) = g (t_{n}) + \frac{u (t_{n - 1})}{Δ t}$ , 省略en, 则在t=tn时刻式 (2) 可转化为如下问题

式 (4) 中, l (u (tn) ) =σu (tn) -εΔu (tn) +aᐁu (tn) 。

事实上, 由于初值u (t0) =u (x, 0) =υ已知, 这实质上已经把一个抛物问题转化成了一个椭圆问题。因此可以用泡函数法来求解该问题, 得到问题式 (2) 在tn时刻的全离散格式:

求u $_{h}^{n}$ ∈V $_{h}^{n}$ 满足

式 (5) 中, $f_{h}^{n} = g (t_{n}) + \frac{u_{h}^{n - 1}}{Δ t}, V_{h}^{n} = V_{L} + V_{B}^{n}$ 为求解空间, VL⊂H $_{0}^{1}$ (Ω) 为前述分片多项式空间V $_{B}^{n}$ 是由VL按第二节提供的方法在tn层上构造的泡函数空间, 由于在不同的时间层上, 泡函数空间是不同的。因此V $_{B}^{n}$ 与n有关, 从而V $_{h}^{n}$ 也与n有关。

由lax-Milgram引理, 问题式 (5) 的解是存在唯一的。

2 误差分析

设问题式 (2) 和式 (5) 的解分别为u (tn) 和u $_{h}^{n}$ , 并记ζn=u (tn) -u $_{h}^{n}$ , 得到‖ζn‖0, 2Ω的估计。

定义1 椭圆投影Ph:H $_{0}^{1}$ (Ω) →VL+H10 (Ωh) 满足

ε (ᐁ (Phu-u) , ᐁχ) + (aᐁ (phu-u) , χ) =0,

∀χ∈VL+H $_{0}^{1}$ (Ωh) 。

其中, H $_{0}^{1}$ (Ωh) =♁TH $_{0}^{1}$ (T) 是由在所有单元边界上均为零的函数构成的H1 (Ω) 的一个子空间, Ωh:=∪{T, T∈Th}但不包括单元边界。由Lax-Milgram引理, 对给定的u, Phu是存在唯一的。

为得到‖Phu-u‖0, 2, Ω和‖Phu-u‖1, 2, Ω的估计, 先引入如下引理。

引理1[4] 存在常数C>0, 对任意有界多角形区域K⊂Rn, α≥0, ε>0, U∈H1 (K) , 均有函数z∈H $_{0}^{1}$ (K) 满足

$\begin{array}{l} a ε^{- \frac{1}{2}} ∥ U - z ∥_{0, 2, k} + ε^{\frac{1}{2}} | U - z |_{1, 2, Κ} \leq \\ C {α^{\frac{1}{2}} ∥ U ∥_{0, 2, \partial Κ} + ε^{\frac{1}{2}} | U |_{1, 2, Κ}} 。 \end{array}$

定理1 令VL=V $_{L}^{k}$ , 整数k≥1, 0<r≤1, 若u∈Hk+r (Ω) ∩H $_{0}^{1}$ (Ω) , 则有

$ε^{\frac{1}{2}} | u - Ρ_{h} u |_{1 ‚ 2 ‚ Ω} \leq$

C $\sum_{t \in Τ_{h}} h_{Τ}^{2 (k - 1)} (ε h_{Τ}^{2 r} + | a | h^{2 r_{Τ}})$

|u| $_{k + r, 2, Τ}^{2}$ $^{\frac{1}{2}}$

$\begin{array}{l} \leq C (ε h^{2 k + 2 r - 2} + \\ | a | h^{2 k + 2 r - 1})^{\frac{1}{2}} | U |_{k + 2 r, 2, Ω} 。 \end{array}$

$ε ∥ u - Ρ_{h} u_{h} ∥_{0, 2, Ω} \leq C (ε h^{2} + | a | h^{3})^{\frac{1}{2}}$

$\sum_{t \in Τ_{h}} h_{Τ}^{2 (k - 1)} (ε h_{Τ}^{2 r} + | a | h_{Τ}^{2 r + 1})$

$\begin{array}{l} | u |_{k + r, 2, Τ}^{2}}^{\frac{1}{2}} \leq C (ε h^{2} + \\ | a | h^{3})^{\frac{1}{2}} C (ε h^{2 k + 2 r - 2} + \\ | a | h^{2 k + 2 r - 1})^{\frac{1}{2}} | U |_{k + 2 r, 2, Ω} 。 \end{array}$

其中|a|:=ess supx∈Ω|a (x) |。

证明取任意的 $\hat{u} \in V_{L} + Η_{0}^{1} (Ω_{h}) ‚ z \in Η_{0}^{1} (Ω_{h})$ , 记u: $= u - \hat{u}$ 。

由diva=0得

(aᐁ (u-Phu) , u-phu) =0。

结合Pu的定义有

ε|u-Phu| $_{1, 2, Ω}^{2}$ =ε (ᐁ (u-Phu) , ᐁ (u-Phu) ) +

(aᐁ (u-Phu) , u-Phu) =ε (ᐁ (u-Phu) ,

ᐁ (U-z) ) + (aᐁ (u-Phu) , U-z) ) ≤

$\begin{array}{l} C ε^{\frac{1}{2}} | u - Ρ_{h} u |_{1, 2, Ω} {\sum_{Τ \in Τ_{h}} (ε^{\frac{1}{2}}) | u - z |_{1, 2, Τ} + \\ ε^{\frac{1}{2}} | a | ∥ U - z ∥_{0, 2, Τ})^{2}}^{\frac{1}{2}} 。 \end{array}$

再由引理 (式) 1

$\begin{array}{l} ε^{\frac{1}{2}} | u - Ρ_{h} u |_{1, 2 Ω} \leq C {\sum_{Τ \in Τ_{h}} (| a |^{\frac{1}{2}} ∥ U ∥_{0, 2, \partial Τ} + \\ ε^{\frac{1}{2}} | U |_{1 ‚ 2 ‚ Τ})^{\frac{1}{2}}}^{\frac{1}{2}} (6) \end{array}$

再由 $\hat{u}$ 的任意性, 取 $\hat{u}$ 为u在V $_{L}^{Κ}$ 上的某种插值 $\hat{u} = Ι_{L} u$ , 满足

将式 (7) 代入式 (6) 即得 $ε^{\frac{1}{2}} | u - p_{h} u h |_{1, 2, Ω}$ 的估计。

再利用对偶论证, 得ε‖u-Phuh‖0, 2, Ω的估计。

假设ϕ∈L2 (Ω) , ψ∈H2 (Ω) ∩H $_{0}^{1}$ (Ω) 是问题

的解, 则有正则性估计ε‖ψ‖2, 2, Ω≤C

‖φ‖0, 2, Ω。由

(u-Phu, ϕ) =-ε (u-Phu, Δψ) - (u-Phu, a·ᐁψ) =

ε (ᐁ (u-phu) , ᐁψ) + (aᐁ (u-Phu) , ψ) 。

以及Phu的定义, 取 $\bar{ψ} = ψ - Ι_{L} ψ$ , 则对∀z∈H $_{0}^{1}$ (Ωh) 利用引理1有

$\begin{array}{l} (u - Ρ_{h} u, ϕ) = ε (\nabla (u - Ρ_{h} u), \nabla (\bar{ψ} - z)) + (a \times \\ \nabla (u - Ρ_{h} u), \bar{Ψ} - Ζ) \leq C ε^{\frac{1}{2}} | u - \end{array}$

$\begin{array}{l} Ρ_{h} u |_{1, 2 Ω} {\sum_{Τ \in Τ_{h}} (ε^{\frac{1}{2}}) | \bar{ψ} - z |_{1, 2, Τ} + \\ ε^{\frac{1}{2}} | a | ∥ \bar{ψ} - z ∥_{0, 2, Τ})^{2}}^{\frac{1}{2}} \leq \end{array}$

$\begin{array}{l} C ϵ^{\frac{1}{2}} | u - Ρ_{h} u |_{1, 2 Ω} {\sum_{Τ \in Τ_{h}} (ϵ^{\frac{1}{2}}) | \bar{ψ} - \\ Ι_{L} ψ |_{1, 2, Τ} + | a |^{\frac{1}{2}} ∥ ψ - \\ i_{L} ψ ∥_{0, 2, \partial Τ})^{2}}^{\frac{1}{2}} \leq C (ε h^{2} + \end{array}$

$| a | h^{3})^{\frac{1}{2}} ∥ ψ ∥ 2, 2, Ω ϵ^{\frac{1}{2}} | u -$

Phu|1, 2, Ω。

因此有

ε‖u-Phu‖0, 2, Ω=0≠0supϕ∈L2 (Ω)

$\begin{array}{l} \frac{ε (u - Ρ_{h} u, ϕ)}{∥ ϕ ∥_{0, 2, Ω}} \leq \\ (ε h^{2} + | a | h^{3})^{\frac{1}{2}} ε^{\frac{1}{2}} | u - \end{array}$

Phu|1, 2, Ω。

再由的 $ε^{\frac{1}{2}} | u - Ρ_{n} u |_{1, 2, Ω}$ 估计, 结论得证。

下面我们来对‖u (tn) -u $_{h}^{n}$ ‖0, 2, Ω进行估计, 即整体L2误差估计。

定理 2令VL=V $_{L}^{k}$ , 整数k≥1, 0<r≤1, 若υ (x) ∈Hk+r (Ω) , ut (x) ∈L2 (0, TE;Hk+r (Ω) ) 则有如下整体误差估计:

$\begin{array}{l} ε ∥ ζ^{n} ∥_{0 ‚ 2 ‚ Ω} \leq C (ε h^{2} + | a | h^{3})^{\frac{1}{2}} (ε h^{2 k + 2 r - 2} + \\ | a | h^{2 k + 2 r - 1})^{\frac{1}{2}} (f_{0}^{t n} | u_{t} |_{k + r, 2, Ω} d s + \\ | υ (x) |_{k + r, 2, Ω}) + ε ∥ u_{0} - υ ∥_{0, 2, Ω + ε} \\ Δ t f_{0}^{t n} ∥ u_{t t} ∥_{0, 2, Ω} d s 。 \end{array}$

证明令u (tn) -u $_{h}^{n}$ =u (tn) -Phu (tn) +Phu (tn) -u $_{h}^{n}$ =瘙窞n+ξn, 则有‖u (tn) -u $_{h}^{n}$ ‖0, 2, Ω≤‖瘙窞n‖0, 2, Ω+‖ξn‖0, 2, Ω。

而由定理1,

$\begin{array}{l} ∥ ■^{n} ∥_{0 ‚ 2 ‚ Ω} \leq C (ε h^{2} + | a | h^{3})^{\frac{1}{2}} (ε h^{2 k + 2 r - 2} + \\ | a |^{2 k + 2 r - 1})^{\frac{1}{2}} | u |_{k + r, 2 r Ω} \leq C (ε h^{2} + \\ | a | h^{3})^{\frac{1}{2}} (ε h^{2 k + 2 r - 2} + | a | h^{2 k + 2 r - 1})^{\frac{1}{2}} \\ (| υ |_{k + r, 2, Ω} + \int_{0}^{t n} | u_{t} |_{k + r, 2, Ω} d s) 。 \end{array}$

因此只需估计‖ξn‖0, 2, Ω。

由式 (3) , 式 (5) 以及Phu (tn) 的定义可得误差方程

$\begin{array}{l} (u_{t} - \partial_{t} u (t_{n}), χ) + \frac{■^{n} - ■^{n - 1}}{Δ t}, χ) + (\frac{ξ^{n} - ξ^{n - 1}}{Δ t}, χ) + \\ ε (\nabla ξ^{n}, \nabla χ) + (a \nabla ξ^{n}, χ) = 0 (9) \end{array}$

式 (9) 中,

$\partial_{t} u (t_{n}) = \frac{u (t_{n}) - u (t_{n - 1})}{Δ t}, \forall χ \in V_{L} + Η_{0}^{1} (Ω_{h})$ 。

令χ=ξn, 则有

‖ξn‖ $_{0, 2, Ω}^{2}$ -‖ξn‖0, 2, Ω‖ξn-1‖0, 2, Ω≤

‖瘙窞n-瘙窞n-1‖0, 2, Ω‖ξn‖0, 2, Ω+Δt‖ (ut-∂tu (tn) ‖0, 2, Ω

‖ξn‖0, 2, Ω。

从而

ε‖ξn‖0, 2, Ω-ε‖ξn-1‖0, 2, Ω≤

ε‖瘙窞n-瘙窞n-1‖0, 2, Ω+εΔt‖ (ut-∂tu (tn) ‖0, 2, Ω。 (10)

把ε‖瘙窞n-瘙窞n-1‖0, 2, Ω=ε‖ (u (tn) -u (tn-1) ) -

Ph (u (tn) -u (tn-1) ) ‖0, 2, Ω≤

$\begin{array}{l} ε \int_{t n - 1}^{t n} ∥ u_{t} - Ρ_{h} u_{t} ∥_{0, 2, Ω} d s \leq C (ε h^{2} + \\ | a | h^{3})^{\frac{1}{2}} (ε h^{2 k + 2 r - 2} + | a | h^{2 k + 2 r - 1})^{\frac{1}{2}} \int_{t_{n - 1}}^{t_{n}} | u_{t} |_{k + r, 2 Ω d s} 。 \end{array}$

和

Δt‖ (ut-∂tu (tn) ‖0, 2, Ω=‖∫ $_{t_{n - 1}}^{t_{n}}$ utt (s) (s-

tn-1) ds‖utt‖0, 2, Ω≤Δt∫ $_{t_{n - 1}}^{t_{n}}$ ‖utt‖0, 2, 52ds,

代入式 (10) 并对n求和, 有

$\begin{array}{l} ε ∥ ξ^{n} ∥_{0, 2, Ω} \leq ε ∥ ξ^{0} ∥_{0, 2, Ω} + ε Δ t \int_{0}^{t_{n}} ∥ u_{t t} ∥_{0, 2, Ω} d s + \\ C (ε h^{2} + | a | h^{3})^{\frac{1}{2}} (ε h^{2 k + 2 r - 2} + \\ | a | h^{2 k + 2 r - 1})^{\frac{1}{2}} \int_{0}^{t_{n}} | u_{t} |_{k + r, 2 Ω} d s 。 \end{array}$

再由

$\begin{array}{l} ε ∥ ξ^{0} ∥_{0, 2, Ω} = ε ∥ u_{0} - Ρ_{h} v ∥_{0, 2, Ω} \leq ε ∥ u_{0} - \\ v ∥_{0, 2, Ω} + ε ∥ Ρ_{h} υ - υ ∥_{0, 2, Ω} \leq \\ ε ∥ u_{0} - υ ∥_{0, 2, Ω} + C (ε h^{2} + | a | h^{3})^{\frac{1}{2}} \\ (ε h^{2 k + 2 r - 2} + | a | h^{2 k + k r - 1})^{\frac{1}{2}} | υ |_{k + r, 2 Ω} 。 \end{array}$

可得ε‖ξn‖0, 2, Ω的估计。利用三角不等式结论得证。

3 数值算例

由理论分析可知, 当区域Ω为一维空间, VL为线性多项式空间时, 有H $_{0}^{1}$ (Ω) ♁VL=H $_{0}^{1}$ (Ω) , 这意味着有限元解u $_{h}^{n}$ 使问题式 (4) 在空间H $_{0}^{1}$ (Ω) 上是精确成立的, 因此有限元解u $_{h}^{n}$ 在空间上是精确逼近真解u (tn) 的, 在时间维上收敛阶为O (Δt) 。

在单位矩形Ω=[0, 1]×[0, 1]上考虑式 (11) 。

设其真解为 $u (x, y, t) = 1 - e - \frac{\sqrt{x y (1 - x) (1 - y) t}}{ε},$

g (x, y, t) 由真解确定。对该问题, 取ε=0.001, 时间步长取t=0.05, 将Ω 做均匀剖分, 步长为h=0.05, 采用双线性元。再分别利用标准有限元和泡函数法求解, 计算结果如表1。

由表1可见, 标准有限元方法产生了数值振荡, 而泡函数法则避免了了数值振荡, 而且数值解的逼近度也有了明显改善。综上可知, 泡函数法可以解决本文提出的抛物型强对流—扩散方程, 能够得到收敛的有限元解, 从而成功地把泡函数法推广到抛物问题。

摘要：传统的 (Galerkin) 方法在解决强对流扩散问题时, 格式稳定性差并伴有强烈的数值震荡现象。为克服上述缺陷, 人们提出了借助于泡函数构造数值方法以得到稳定化的数值解, 并成功地应用于求解椭圆型强对流扩散问题。则把泡函数法推广到抛物型强对流扩散问题, 给出了最优的L2误差分析, 得到了稳定化的数值解。

关键词：泡函数法,抛物型强对流扩散问题,收敛性

参考文献

[1] Breezi F, Russo A.Choosing bubles for advection-diffusion problems.Math Models Methods Appl Sci, 1994;44:571—587

[2] Breezi F, Fanca L P, Hughes T J R, et al.b=∫g, comput.Methods Appl Mech Engrg, 1997;145;329—339

[3]Brezzi F, Hughes TJ R, Marini D, et al.Apriori error analysis of re-sidual-free bubbles for advection-diffusion problem.SIAM J Number Anal, 1990;36:51—63

扩散方程篇8

在分析P-M模型扩散系数的基础上, 讨论了扩散系数对图像去噪的影响, 然后改进了P-M模型中的扩散系数, 实验结果表明, 改进的扩散系数具有更好的去噪效果。

1 P-M模型

由于各向同性的热传导方程在整个图像范围内采用相同的扩散系数, 所以会出现边界模糊的问题。因此, 考虑通过选择扩散系数, 使扩散程度在各个方向上不同来解决这个问题。考虑将图像的梯度值作为边缘检测算子, 使扩散系数为图像梯度的函数, 随梯度值变化而变化, 在梯度小的地方具有较大值, 随着梯度的增大, 扩散系数逐渐减小, 这样, 就能很好平滑非边缘处的噪声, 而又能保留边缘处的特征。基于上述思想, Perona和Mailik提出了非线性的扩散模型, 简称P-M模型[9], 其表达式为

式 (1) 中, u0表示原图像; (x, y) 表示图像的空间坐标;变量t表示处理顺序的参数, 在离散过程中表示迭代次数;div表示散度算子;表示梯度算子;表示图像的梯度模。是非增函数, 称为扩散系数, 也是尺度函数, 相当于一个边缘检测算子, 用于控制扩散速度, 与图像梯度成反比, 在图像变化光滑的区域进行快速各向异性扩散, 在图像边缘区进行低速扩散甚至不扩散, 达到在去噪的同时保持边缘信息的效果, 满足C (0) =1, 且时退化为各向同性热传导方程。由以上分析可以看出, 扩散系数选择的正确与否决定了扩散方程是否有效。Perona和Mailik给出了式 (1) 中两个既能去掉噪声又能增强边缘的扩散系数函数。

式中, K为梯度阈值, 用于判断图像的特征, 它的取值对于图像的去噪效果起到关键性的作用, 当时, 扩散系数受到抑制, 能够保留边缘信息甚至增强, 时, 扩散强度较大, 能够很好的去噪。若K的取值太大, 容易将边缘点误判为噪声点, 图像会过度平滑导致图像模糊, 若K的取值太小, 不能很好的识别出噪声点, 图像会较早停止扩散导致去噪效果较差, 目前很难找到一个能够良好区分噪声和边缘点的梯度阈值K的方法, 根据图像和高斯噪声分布的特性, 如果图像的信噪较小, K的取值较大;相反, 如果图像的信噪比大, K的取值较小。

P-M模型在理论上是逆向扩散过程, 这种逆向扩散虽然能在一定条件下增强图像边缘的强度, 取得较好的去噪效果, 但是该方程在理论上是不适定的, 不适当的参数设置可能导致方程产生完全畸变的结果。在噪声处, 图像的梯度可能非常大, 此时平滑系数较小, 噪声点被保留下来。

2 CLMC模型

针对P-M模型的不足, 1992年, Cattle, Lions, Morel和Coll对P-M模型作出了改进, 提出了用高斯滤波来优化P-M模型 (CLMC模型) [10], 平滑梯度值较大的噪声点, 降低梯度值, 原图像边缘区域被保留下来, 接下来再运用P-M模型进行处理, 就能平滑梯度值较大的噪声点, 得到很好的去噪效果, 解决了P-M模型存在的某些问题, 在理论上证明该模型的解是适定性[10]。CLMC模型表达式为

是高斯函数[10], σ为高斯函数的滤波尺度。当σ→0时, CLMC模型就是P-M模型, 当σ>0时, CLMC模型中的扩散系数通过在空间域中引入正则化项, 不仅克服了P-M模型中梯度易受噪声影响的缺点, 而且是适定的。

在CLMC模型中, σ的选取是一个难点也是该模型的缺点, 若取值过小, 就失去了高斯滤波的作用, 并且导致扩散过程发散;若取值过大, 则容易造成图像模糊, 丢失细节特征。

3 扩散系数的改进

对P-M模型中的扩散系数修改为

式 (5) 中, 0<k≤1。下面给出C1、C2及不同常数k下新的扩散系数C3的函数图像, 取梯度阈值K=10。

从图像上可以看出, 扩散系数C1在图像的整个梯度范围内进行去噪, 随着的增加, 扩散系数逐渐减小;当很大时, 扩散强度才逐渐趋于零, 使图像在整个区域进行去噪, 同时也模糊了图像的边缘部分。而扩散系数C2在较小时, 有较大的扩散能力, 较大时, 几乎没有扩散作用, 保持了图像的边缘, 同时导致边缘处的噪声被保留。新的扩散系数C3在常数k较小时, 同时具备C1和C2的优点, 既能保持图像的边缘, 又能去除噪声。数值试验结果也说明了这个问题。

4 数值计算方法

将修改的扩散系数应用到CLCM模型, 采用有限差分对式 (4) 进行离散[11,12]。

式中, k表示迭代次数;λ表示时间步长与空间步长之比, λ≤0.5。

5 实验分析

5.1 图像处理效果的评价指标

对于去噪效果衡量, 有主观评价方法和客观评价方法。主观评价方法就是目测方法, 在实际应用中不够精确;客观评价方法是通过图像的方差、均值、信噪比和峰值信噪比等判断。

(1) 均方误差 (mean squared error, MSE)

(2) 信噪比 (signal noise rate, SNR)

(3) 峰值信噪比 (peak signal noise rate, PSNR)

式中, u (i, j) 表示原图像;表示去噪后的图像;M和N表示图像的长和宽。MN即表示图像所有像素点, SNR和PSNR越大表示图像失真越少。以上三个评价指标基本上能够反映出图像的去噪效果。

5.2 实验结果

程序采用MATLABR2010a编写, 实验中选取标准测试图像Lena, 加入均值为0, 方差为0.01的高斯白噪声, 取梯度阈值K=10, 步长λ=0.05, 迭代次数为5, 并选取σ=0.5, 3×3领域的高斯核。图2为本文扩散系数C3 (k=0.2) 去噪与P-M方程扩散系数C1和C2在CLMC模型中去噪效果对比图。

从表1中的数据可以看出, 使用新的扩散系数C3进行去噪, 信噪比 (SNR) 和峰值信噪比的值均比C1和C2的值大, 这说明新的扩散系数去噪效果较好。同时常数k的取值也会影响图像的去噪效果, 从表1可以看出, 当常数k的取值较小时图像的去噪效果较好。实验表明当k>1时, 新的扩散系数的去噪效果没有P-M方程所给出的扩散系数去噪效果好。

6 总结

通过研究PDE方程在图像去噪领域的基础模型:P-M模型和CLMC模型, 重点分析了扩散系数对图像去噪的影响, 并且提出了一种新的扩散系数, 应用到CLMC模型, 采用差分的数值求解, 实验结果表明了新扩散系数的有效性。

参考文献

[1] Li S Z.Markov random field modeling in computer version.3rd ed.London:Spring Verlag, 1999:50—60

[2] Donoho D L.Denoising by soft-thresholding.IEEE Transactions on Iformation Theory, 1995;41 (33) :613—627

[3] Mallat S.A wavelet tour of signal processing.2nd ed.San Diego:Elsevier, 1997:486—491

[4] Chen G Y, Bui T D.Multiwavelets denoising using neighboring coefficients.Signal Processing Letters, IEEE, 2003;10 (7) :211—214

[5] Gabor D.Information theory in electron microscocpy.Laboratory Investigation, 1965;45 (14) :801—807

[6] Jain A K.Partial differential equations and finite-difference methods in image processing.Optimization Theory and Applications, 1977;41 (2) :65—91

[7] Koenderink J J.The structure of images.Biological Cybernetics, 1984;50 (5) :363—370

[8] Witkin A D.Scale-space filtering.Proceedings of the International Joint Conference on Artificial Intelligence, ACM Inc, 1983;39 (2) :1019—1021

[9] Perona P, Malik J.Scale-space and edge detection using anisotropic diffusion.IEEE Trans Pat tern Anal Machine Intel, 1990;12 (7) :629 —639

[10] Catte F, Lions P L, Morel J M, et al.Image selective smoothing and edge detection by nonlinear diffusion.SIAM Journal on Numerical Analysis, 1992;29 (1) :182—193

[11] 王大凯, 侯榆青, 彭进业.图像处理的偏微分方程方法.北京:科学出版社, 2008Wang D K, Hou Y Q, Peng J Y.Image processing method of partial differential equations.Beijing:Science Press, 2008

本文来自 360文秘网(www.360wenmi.com)，转载请保留网址和出处

【扩散方程】相关文章：

对流扩散05-22

扩散处理05-24

扩散效应05-25

空间扩散06-14

扩散理论07-08

行为扩散07-18

扩散模拟07-19

信息扩散08-09

扩散性能08-15

扩散体系08-16