一個時間分數階擴散方程的源強識別反問題

殷鳳蘭 尹曉燕

摘 要： 對于一類帶有多個點源的分數階維擴散方程問題，應用有限差分法給出了一個數值求解格式，并在已知點源個數及其位置的前提下，應用最佳攝動量正則化算法對源強進行了數值反演，并討論了擾動數據和微分階數的不同選取對反演算法的影響.

關鍵詞： 分數階擴散方程 多點源 最佳攝動量正則化算法

整數階擴散方程是根據質量守恒定律和費克梯度擴散定律得出的，但是實際情況下許多擴散現象并不滿足經典的費克梯度擴散定律，通常稱之為“反常”擴散.分數階導數與反常擴散過程本質上有一定的相似性，所以分數階擴散方程應運而生，其研究已受到人們越來越多的關注.對于分數階微分方程的求解，由于解析解一般都含有特殊函數或者復雜級數的形式，計算起來比較困難，因此越來越多的研究者討論分數階微分方程的數值解，而對于分數階擴散方程的源項反演問題卻少有研究.本文研究多點源分數階對流-擴散方程的差分求解方法，進而應用最佳攝動量正則化算法對分數階擴散方程中的多點源源強進行了數值反演.

考慮一個帶多個點源的時間分數階對流擴散方程的源強反演問題：

■-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.