殷鳳蘭 尹曉燕

摘 要： 對(duì)于一類帶有多個(gè)點(diǎn)源的分?jǐn)?shù)階維擴(kuò)散方程問(wèn)題，應(yīng)用有限差分法給出了一個(gè)數(shù)值求解格式，并在已知點(diǎn)源個(gè)數(shù)及其位置的前提下，應(yīng)用最佳攝動(dòng)量正則化算法對(duì)源強(qiáng)進(jìn)行了數(shù)值反演，并討論了擾動(dòng)數(shù)據(jù)和微分階數(shù)的不同選取對(duì)反演算法的影響.

關(guān)鍵詞： 分?jǐn)?shù)階擴(kuò)散方程 多點(diǎn)源 最佳攝動(dòng)量正則化算法

整數(shù)階擴(kuò)散方程是根據(jù)質(zhì)量守恒定律和費(fèi)克梯度擴(kuò)散定律得出的，但是實(shí)際情況下許多擴(kuò)散現(xiàn)象并不滿足經(jīng)典的費(fèi)克梯度擴(kuò)散定律，通常稱之為“反?！睌U(kuò)散.分?jǐn)?shù)階導(dǎo)數(shù)與反常擴(kuò)散過(guò)程本質(zhì)上有一定的相似性，所以分?jǐn)?shù)階擴(kuò)散方程應(yīng)運(yùn)而生，其研究已受到人們?cè)絹?lái)越多的關(guān)注.對(duì)于分?jǐn)?shù)階微分方程的求解，由于解析解一般都含有特殊函數(shù)或者復(fù)雜級(jí)數(shù)的形式，計(jì)算起來(lái)比較困難，因此越來(lái)越多的研究者討論分?jǐn)?shù)階微分方程的數(shù)值解，而對(duì)于分?jǐn)?shù)階擴(kuò)散方程的源項(xiàng)反演問(wèn)題卻少有研究.本文研究多點(diǎn)源分?jǐn)?shù)階對(duì)流-擴(kuò)散方程的差分求解方法，進(jìn)而應(yīng)用最佳攝動(dòng)量正則化算法對(duì)分?jǐn)?shù)階擴(kuò)散方程中的多點(diǎn)源源強(qiáng)進(jìn)行了數(shù)值反演.

考慮一個(gè)帶多個(gè)點(diǎn)源的時(shí)間分?jǐn)?shù)階對(duì)流擴(kuò)散方程的源強(qiáng)反演問(wèn)題：

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

其中u=u（x，t）為t時(shí)刻x處污染物的濃度分布，D為擴(kuò)散系數(shù)，D≥0，v為對(duì)流系數(shù)，v≥0，f（x，t）=■Q■δ（x-x■），x■為污染源的坐標(biāo)位置，Q■為污染物排放強(qiáng)度（單位時(shí)間排放量），M為點(diǎn)污染源的個(gè)數(shù)，δ為狄拉克函數(shù).■是Caputo意義下的分?jǐn)?shù)階導(dǎo)數(shù)，定義為

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

式中Γ表示Γ函數(shù)，Γ（x）=e■t■，且0<α≤1.

初邊值條件給定為：

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

假設(shè)知道污染源的個(gè)數(shù)，以及每個(gè)源的位置坐標(biāo)，而需要確定各個(gè)污染源的強(qiáng)度值.這時(shí)，需要額外的關(guān)于污染物濃度分布的附加數(shù)據(jù)，聯(lián)合模型（1）—（2）形成一個(gè)源強(qiáng)度識(shí)別的反問(wèn)題.設(shè)在出流端t=T處觀測(cè)到多處的濃度值為：

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

這樣，由附加數(shù)據(jù)（3）式聯(lián)合模型（1）—（2）式即構(gòu)成一維多點(diǎn)源時(shí)間分?jǐn)?shù)階擴(kuò)散模型的源強(qiáng)識(shí)別反問(wèn)題.

1.正問(wèn)題的求解

下面應(yīng)用隱式差分方法給出問(wèn)題的數(shù)值解.設(shè)x■=ih，h>0；i=0，1，2，…，M；t■=nτ，τ>0，n=0，1，2，…，N，其中h和τ分別是空間和時(shí)間步長(zhǎng)，且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■）.設(shè)p=■，q=■，u■■是u（x■，t■）的數(shù)值近似，即得

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

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

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

當(dāng)n=0時(shí)，

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

當(dāng)n>0時(shí)，

-（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.最佳攝動(dòng)量正則化算法

最佳攝動(dòng)量正則化算法是一種基于Tikhonov正則化、求解參數(shù)反演問(wèn)題的迭代方法.對(duì)于源強(qiáng)反問(wèn)題（1）—（2），設(shè)污染源的位置已知，即x■，x■，…，x■為，需要反演相應(yīng)的源強(qiáng)度Q■，Q■，…，Q■.記R=[Q■，Q■，…，Q■.]，在R已知的前提下，借助上述求解正問(wèn)題的方法，求出正問(wèn)題（1）—（2）的解u（x，t；R），進(jìn)而得到其在時(shí)刻t=T和x■處的值，記為u（x■，t；R）.這樣求解反演問(wèn)題（1）—（2）的一種優(yōu)化方法就是使得計(jì)算值與觀測(cè)值在某種誤差意義下最小，通常化為下述極小問(wèn)題的求解

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

其中μ為正則參數(shù).根據(jù)攝動(dòng)算法，極小問(wèn)題的求解又轉(zhuǎn)化為對(duì)于給定的R■，求解最佳攝動(dòng)量δR■，進(jìn)而由下式確定出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■處作泰勒展開(kāi)，略去高階項(xiàng)，近似可得

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

對(duì)t進(jìn)行離散0=x■

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

再根據(jù)最小二乘法的思想，求解minF（δR■）相當(dāng)于求解規(guī)范方程

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

其中

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

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

γ稱為數(shù)值微分步長(zhǎng).以下給出反演R的算法步驟：

步驟1.給定初始猜測(cè)向量R■和數(shù)值微分步長(zhǎng)γ，求向量η，ξ及矩陣G；

步驟2.按照通常的最佳攝動(dòng)量算法進(jìn)行計(jì)算，由下式求出δR■

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

步驟3.對(duì)于給定的精度eps，判斷是否滿足||δR■||≤eps.若是，則R■即為所求，算法終止；否則，由（12）式得到R■，再轉(zhuǎn)到步驟1繼續(xù)進(jìn)行.

3.數(shù)值模擬

問(wèn)題（1）—（2）中，取初始函數(shù)取為g（x）=10x，取M=100，N=100，擴(kuò)散系數(shù)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，即源強(qiáng)真值為R=（1，3，5，7）.設(shè)微分階數(shù)α=0.8.其中δ表示擾動(dòng)水平，反演誤差用Err=||R■-R||■表示，n表示迭代次數(shù)，R■=（Q■■，Q■■，Q■■，Q■■）表示反演解.反演結(jié)果如下：

利用反演得到的源強(qiáng)進(jìn)一步模擬t=10時(shí)的污染物濃度分布，并與t=10時(shí)實(shí)際的污染物濃度分布進(jìn)行對(duì)比，見(jiàn)圖1、圖2.

考察不同的微分階數(shù)α對(duì)算法的影響，計(jì)算結(jié)果見(jiàn)表2：

通過(guò)上述數(shù)值模擬可以看出，此算法穩(wěn)定性較好，即使附加數(shù)據(jù)有較大擾動(dòng)所得反演結(jié)果仍然比較接近真實(shí)值，微分階數(shù)對(duì)反演算法的影響也不大.

參考文獻(xiàn)：

[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]蘇超偉.偏微分方程逆問(wèn)題的數(shù)值解法及其應(yīng)用[M].西安：西北工業(yè)大學(xué)出版社，1995.