張 蘭
(西安航空職業技術學院,陜西 西安 710089)
隨著數學理論和計算技術的快速發展,一些污染擴散問題可以通過數學中的反問題進行數值模擬[1-2],如可以通過反問題的計算來確定污染源等[3-8]。然而客觀世界中的污染存在著漫擴散現象,無法由經典的整數微積分模型來描述。分數階擴散模型是由學者Maneldbrot[9]提出的,其具有良好的記憶性,可以用來描述污染物的擴散,模擬實際現象。因此,尋找行之有效的分數階擴散方程的反演方法對于治理污染問題有非常大的意義和研究價值,也可以更好地解決具有漫擴散現象的其他問題。
目前,對分數階擴散方程的研究主要集中在正問題方面,對反演問題的研究相對較少。谷文娟等[10]利用最佳攝動量正則化方法對分數階擴散方程的擴散系數進行了反演計算;王俊剛[11]運用幾種Tikhonov正則化方法求解了時間分數階擴散方程的反初值和反源項問題;Zhang等[12]將反問題轉為優化問題,對分數階擴散方程的源項反演問題進行了數值模擬。在分數階擴散方程的反演問題中,正則化方法是使用較多的方法,通過正則化方法將反問題轉化為泛函極值問題,可以很好地降低反問題的求解難度。而正則化方法中,解的精度又依賴于參數的合理選取[13-18],因此選擇合適的正則化參數是分數階擴散方程反演問題的重要一環。
近年來,隨著智能算法具有快速發展和易收斂于全局最優解的特點,越來越多的智能算法被應用于優化問題。量子粒子群優化(quantum-behaved particle swarm optimization, QPSO)算法是2004年由Sun等[19]在粒子群的基礎上提出的一種智能優化算法,算法中粒子在量子空間中搜索最優點,該算法由于具有較快的收斂速度和良好的全局收斂性被廣泛應用于工程中,但量子粒子群算法應用于分數階擴散方程的反演問題鮮有報道。本文針對分數階擴散方程反初值問題,采用Tikhonov正則化計算,將量子粒子群優化算法應用于求解Tikhonov正則化參數,數值仿真實驗表明該算法是有效的。
粒子群優化(particle swarm optimization, PSO)算法是1995年由美國學者Kennedy等[20]提出的一種智能優化算法,該算法模擬鳥群的覓食行為,通過不斷更新粒子的最優個體和最優群體來尋找最優解。該算法簡單,收斂速度快,易于實現,是一種通用性很強的優化算法,能夠用于求解各種復雜的優化問題。Sun在量子機制下提出的量子粒子群優化算法不用計算目標函數的梯度,只需要知道目標函數的值即可。通過實驗測試發現,量子粒子群優化算法能更好地收斂全局最優解。
在量子粒子群優化算法中,每一個粒子會向局部吸引點p聚集,粒子的狀態采用波函數描述,假設第i個粒子的位置是Xi,其個體最優位置為pbesti,群體最優位置為gbesti,群體數為M,引入群體平均最優位置mbest,則粒子的更新公式如下:
pi(t+1)=αpbesti(t)+(1-α)gbesti(t)
(1)
(2)
當u≥0.5時,有
Xi(t+1)=pi(t+1)-b|mbest(t)-Xi(t)|ln(1/u(t+1))
(3)
當u<0.5時,有
Xi(t+1)=pi(t+1)+b|mbest(t)-Xi(t)|ln(1/u(t+1))
(4)
式中:α和u為(0,1)的隨機數;pi(t+1)為(t+1)時刻第i個粒子的最佳位置;b為收縮/擴張系數,b=1-0.5t/maxiter,maxiter為最大迭代次數。
本文取如下一維時間分數階方程:
x∈Ω,t∈(0,T),0<α<1
(5)

u(x,t)=0,x∈?Ω,t∈(0,T)
u(x,0)=f(x),x∈Ω
(6)
(7)
式中:?Ω為區域Ω邊界;Γ為Gamma函數;u(x,0)為t=0時初值;f(x)為源項;τ為特定時刻;uτ(x,t)為特定時刻的值;uτ(x)為α=1時特定時刻的值。
分數階方程正問題如式(5),已知初值u(x,0)=f(x),求解u(x,T)=g(x),其中u(x,T)為終值時刻的觀測數據;反問題則為已知觀測數據gδ(x),求解初值f(x):=u(x,0),其中δ為引入的絕對誤差水平。則正問題式(5)所對應的反初值問題如下:
(8)
其中gδ(x)滿足如下不等式:
‖gδ(x)-g(x)‖≤δ
(9)
式中:‖·‖為歐幾里得范數。反問題具有不適定性,對于分數階反問題相比正問題的求解難度更大,學者Tikhonov[21]提出了解決反問題的正則化方法,通過泛函將反問題轉化為適定的優化問題進行求解,該方法一經出現便很快被應用于各個領域的反問題求解。
對于式(8)所示的反問題,采用Tikhonov正則化方法進行求解,記正問題式(5)的解為L(f(x))=u(x,T),L為對稱一致橢圓算子,則反問題轉化為如下正則化問題:
(10)
式中:λ為正則化系數。對于上述問題式(10),Tikhonov正則化解fλ可以通過下面的線性方程組來求解:
(LTL+λI)fλ=LTgδ
(11)
式中:I為單位矩陣。對于固定的λ(λ>0),式(11)有唯一解,即正則化解fλ可以表示為:
fλ=(LTL+λI)-2LTgδ
(12)
正則化解fλ的穩定性和近似程度與參數λ的選取有關,λ越大,解越穩定;λ越小,解越近似于原問題。因此,參數的選取成為正則化方法的重要一環。參數選擇的主要方法有L準則法、廣義偏差法、Morozov偏差理論法等,這些方法都是通過轉化為非線性方程后再進行最優參數的求解,但需要先選取初值,初值的選取也會影響到結果的精確性,因此本文結合智能優化算法中的量子粒子群算法,在Morozov偏差原理[18]的基礎上進行參數選取。首先將問題轉化為優化問題:
λ=argλmin|‖Lfλ-gδ‖2-δ|
(13)
然后采用QPSO算法來處理式(13)所示的極小值問題,從而獲得最優的正則化參數。

下面給出一個算例[11], 取Ω=(0,1),α(x)=x2+1,c(x)=-(x-1),取真解為f(x)=xα(1-x)1-αexαsin(7πx)。用QPSO算法求解分數階擴散方程的正則化參數,并編寫相應的程序對分數階反問題進行數值計算。取時間方向的N=100、空間方向的M=100計算正問題,取M=50、N=100計算反問題。在Morozov偏差原理下,分別取分數階的次數α為0.3和0.7,誤差ε為0.005,0.050和0.010。運用QPSO算法計算正則化參數λ,當α=0.3時,得到λ分別為1.2×10-8,2.1×10-8,1.1×10-7;當α=0.7時,得到λ分別為6.5×10-10,2.3×10-9,6.2×10-9。通過上述方法求解得到的f(x)如圖1所示。從圖中可以看出,精確解與正則化解吻合較好,因此本文算法是有效的。
正則化方法是求解反問題常用的策略,可將反問題轉化為優化問題來求解。本文將量子粒子群優化算法應用于一維分數階擴散方程的反演,并對反演方法進行了分析,數值算例結果說明采用量子粒子群算法求解分數階擴散方程可以得到高精度的數值解。本文只研究了一維分數階擴散方程的反初值問題,其他類似的問題以及二維分數階擴散方程的反演問題,將是下一步研究的內容。