一个时间分数阶扩散方程的源强识别反问题

殷凤兰 尹晓燕

摘 要： 对于一类带有多个点源的分数阶维扩散方程问题，应用有限差分法给出了一个数值求解格式，并在已知点源个数及其位置的前提下，应用最佳摄动量正则化算法对源强进行了数值反演，并讨论了扰动数据和微分阶数的不同选取对反演算法的影响.

关键词： 分数阶扩散方程 多点源 最佳摄动量正则化算法

整数阶扩散方程是根据质量守恒定律和费克梯度扩散定律得出的，但是实际情况下许多扩散现象并不满足经典的费克梯度扩散定律，通常称之为“反常”扩散.分数阶导数与反常扩散过程本质上有一定的相似性，所以分数阶扩散方程应运而生，其研究已受到人们越来越多的关注.对于分数阶微分方程的求解，由于解析解一般都含有特殊函数或者复杂级数的形式，计算起来比较困难，因此越来越多的研究者讨论分数阶微分方程的数值解，而对于分数阶扩散方程的源项反演问题却少有研究.本文研究多点源分数阶对流-扩散方程的差分求解方法，进而应用最佳摄动量正则化算法对分数阶扩散方程中的多点源源强进行了数值反演.

考虑一个带多个点源的时间分数阶对流扩散方程的源强反演问题：

■-D■+v■=f（x，t） 0

其中u=u（x，t）为t时刻x处污染物的浓度分布，D为扩散系数，D≥0，v为对流系数，v≥0，f（x，t）=■Q■δ（x-x■），x■为污染源的坐标位置，Q■为污染物排放强度（单位时间排放量），M为点污染源的个数，δ为狄拉克函数.■是Caputo意义下的分数阶导数，定义为

？鄣■■u（x，t）=■？蘩■■■■，0<α<1■，α=1

式中Γ表示Γ函数，Γ（x）=e■t■，且0<α≤1.

初邊值条件给定为：

u（x，0）=g（x），u（0，t）=u（l，t）=0（2）

假设知道污染源的个数，以及每个源的位置坐标，而需要确定各个污染源的强度值.这时，需要额外的关于污染物浓度分布的附加数据，联合模型（1）—（2）形成一个源强度识别的反问题.设在出流端t=T处观测到多处的浓度值为：

u（x■，T）=■■，k=1，2，…，K（3）

这样，由附加数据（3）式联合模型（1）—（2）式即构成一维多点源时间分数阶扩散模型的源强识别反问题.

1.正问题的求解

下面应用隐式差分方法给出问题的数值解.设x■=ih，h>0；i=0，1，2，…，M；t■=nτ，τ>0，n=0，1，2，…，N，其中h和τ分别是空间和时间步长，且Mh=l，Nτ=T.在方程（1）中采用如下有限差分近似：

？鄣■■u（x■，t■）=■■■？蘩■■■+O（τ）（4）

=■{u（x■，t■）-u（x■，t■）+■[u（x■，t■）-u（x■，t■）]}[（k+1）■-k■]+O（τ）

■|■=■+O（h■）（5）

■|■=■+O（h）（6）

将式（4）—（6）代入方程（1），可得

■{u（x■，t■）-u（x■，t■）+■[u（x■，t■）-u（x■，t■）]}[（k+1）■-k■]-D■+v■=f■（7）

其中i=1，2，…，M-1，f■=■Q■δ（ih-x■）.设p=■，q=■，u■■是u（x■，t■）的数值近似，即得

-pu■■+（1+2p-q）u■■-（q-p）u■■=u■■-■（u■■-u■■）[（k+1）■-k■]+τ■Γ（2-α）f■（8）

u■■=0，u■■=0，u■■=g（x■）

因此，可得隐式差分格式：

当n=0时，

-（q+p）u■■+（1+2p+q）u■■-pu■■=u■■+τ■Γ（2-α）f■（9）

当n>0时，

-（q+p）u■■+（1+2p+q）u■■-pu■■=（2-2■）u■■-■u■■[2（k+1）■-（k+2）■-k■]+u■■[（n+1）■-n■]+τ■Γ（2-α）f■（10）

2.最佳摄动量正则化算法

最佳摄动量正则化算法是一种基于Tikhonov正则化、求解参数反演问题的迭代方法.对于源强反问题（1）—（2），设污染源的位置已知，即x■，x■，…，x■为，需要反演相应的源强度Q■，Q■，…，Q■.记R=[Q■，Q■，…，Q■.]，在R已知的前提下，借助上述求解正问题的方法，求出正问题（1）—（2）的解u（x，t；R），进而得到其在时刻t=T和x■处的值，记为u（x■，t；R）.这样求解反演问题（1）—（2）的一种优化方法就是使得计算值与观测值在某种误差意义下最小，通常化为下述极小问题的求解

min||u（x■，T；R）-■■||■+μ||R||■■，（11）

其中μ为正则参数.根据摄动算法，极小问题的求解又转化为对于给定的R■，求解最佳摄动量δR■，进而由下式确定出R■的一种迭代算法

R■=R■+δR■，n=0，1，2，…（12）

其中，δR■是下述泛函的极小解

F（δR■）=||u（x■，T；R■+δR■）-■■||■■+μ||δR■||■■.（13）

将u（l，t■；R■+δR■）在R■处作泰勒展开，略去高阶项，近似可得

F（δR■）=||u（x■，T；R■）-■■+？塄■u（l，t■；R■）·δR■||■■+μ||δR■||■■.

对t进行离散0=x■

F（δR■）=■[u（x■，T；R■）]+？塄■u（x■，T；R■）·δR■-■■]■+μδR■■δR■.

再根据最小二乘法的思想，求解minF（δR■）相当于求解规范方程

μI+G■GδR■=G■（η-ξ），（14）

其中

G=（g■）■，g■=■；

ξ=（u（x■，T；R），u（x■，T；R），…，u（x■，T；R））■；η=（■■，■■，…，■■）■，

γ称为数值微分步长.以下给出反演R的算法步骤：

步骤1.给定初始猜测向量R■和數值微分步长γ，求向量η，ξ及矩阵G；

步骤2.按照通常的最佳摄动量算法进行计算，由下式求出δR■

δR■=（αI+G■G）■G■（η-ξ）；（15）

步骤3.对于给定的精度eps，判断是否满足||δR■||≤eps.若是，则R■即为所求，算法终止；否则，由（12）式得到R■，再转到步骤1继续进行.

3.数值模拟

问题（1）—（2）中，取初始函数取为g（x）=10x，取M=100，N=100，扩散系数D=0.001，v=0，q=4，Q■=1，Q■=3，Q■=5，Q■=7，x■=0.3，x■=0.7，x■=0.8，x■=0.9，即源强真值为R=（1，3，5，7）.设微分阶数α=0.8.其中δ表示扰动水平，反演误差用Err=||R■-R||■表示，n表示迭代次数，R■=（Q■■，Q■■，Q■■，Q■■）表示反演解.反演结果如下：

利用反演得到的源强进一步模拟t=10时的污染物浓度分布，并与t=10时实际的污染物浓度分布进行对比，见图1、图2.

考察不同的微分阶数α对算法的影响，计算结果见表2：

通过上述数值模拟可以看出，此算法稳定性较好，即使附加数据有较大扰动所得反演结果仍然比较接近真实值，微分阶数对反演算法的影响也不大.

参考文献：

[1]Tong S J，Zheng W，Chen B Z.Analysis of the pollution consequences on leakage and seepage flow of poisonous liquid[J].Ind ustrial Safety and Environmental Protection，2006，32（10）：56-58.

[2]Chen W.A speculative study of 2/3-order fractional Laplacian modeling of turbulence：some thoughts and conjectures[J].C haos，2006，16：02316.

[3]Wang Sheng，Ma Zhengfei，Yao Huqing.Fourier-Bessel series algorithm in fractal diffusion model for porous material[J].Chinese J.C omputat Phys，2008，25（3）：289-295.

[4]苏超伟.偏微分方程逆问题的数值解法及其应用[M].西安：西北工业大学出版社，1995.