局部化MAPS法求解時空偏微分方程

陳林芳，張學瑩

(河海大學理學院，南京 210098)

一種有效的數(shù)值方法對于求解微分方程是非常關(guān)鍵的，它不僅要能滿足數(shù)值解精度方面的要求，而且還要保證計算效率。目前發(fā)展起來的數(shù)值方法主要有2大類:基于網(wǎng)格數(shù)值解法和無網(wǎng)格方法?；诰W(wǎng)格數(shù)值解法主要有有限差分法(FDM)、有限體積法(FVM)和有限元法(FEM)等;無網(wǎng)格方法主要有光滑粒子流體動力學方法(SPH)［1］、Galerkin 方法(EFG)［2］、重構(gòu)核質(zhì)點方法(RKPM)［3］、微積分方法(DQ)［4］和近似特別解方法(MAPS)［5－7］等。無網(wǎng)格方法在求解積分方程或偏微分方程組時只需要用區(qū)域內(nèi)離散點上的函數(shù)值去近似空間中某點的函數(shù)值和各階導數(shù)值，而這些離散點間無需用網(wǎng)格進行連接，計算格式簡單，能部分或徹底地消除網(wǎng)格劃分帶來的困難，所以無網(wǎng)格方法逐漸受到國內(nèi)外學者的關(guān)注。在眾多的無網(wǎng)格方法中，基于徑向基函數(shù)(radial basis functions，RBFs)的無網(wǎng)格方法更引人注目。最初，Kansa將RBFs引入到求解偏微分方程中，即所謂的Kansa方法［8－9］;1990年Dehghan證明了基于徑向基函數(shù)的Kansa［10］方法，并應用該方法處理了偏微分方程中的多元數(shù)據(jù)，對離散數(shù)據(jù)進行了逼近，最后給出了偏微分方程的數(shù)值解;G.R.Liu提出了徑向基點插值法(RPIM)，有效地將徑向基函數(shù)的方法局部化;C.S.Chen等采用LMQ的思想，在近似特別解(MAPS)的基礎(chǔ)上，對區(qū)域進行局部劃分，得到局部近似特別解(LMAPS)。本文主要用LMAPS法解時空偏微分方程。該方法與Kansa方法有很大相似性，可通過數(shù)值實例進行誤差比較分析。

1 時空偏微分方程的局部MAPS形式

1.1 導數(shù)的局部近似

大規(guī)模的稠密矩陣的病態(tài)矩陣問題已經(jīng)限制了以徑向基函數(shù)為基礎(chǔ)的無網(wǎng)格方法的發(fā)展，而局部近似特別解法是一種比較有效的方法，形參改變時波動較少，并且能夠避免病態(tài)矩陣問題。

這里，Ωp是局部支撐域，并且

其中:‖·‖是歐幾里得范數(shù);φ(‖x－xj‖)是一個徑向基函數(shù)。這里選用MQ-RBF是因為它在大多數(shù)應用中能提供較高的精度，已經(jīng)被研究者廣泛應用。規(guī)范化的MQ-RBFs能夠?qū)懗梢韵滦问?

其中r=‖x－xj‖。把方程(1)重新寫成以下矩陣形式:

此處 Φns=［Φ(‖xp－ x1‖)，Φ(‖xp－ x2‖)，…，Φ(‖xp－ xns‖)］，ans=［a1，a2，…，.ans］T。因為?Ωp，每個節(jié)點都有類似于式(1)的方程，所以得到以下的矩陣式:

這里，Ψns=［Φxk(‖xp－x1‖)，…，Φxk(‖xp－xns‖)］，xk表示第 k層坐標方向，Φxk表示 Φ 關(guān)于 xk的一階偏微商，可計算求得:

同理，u(xp)的二階微商可近似成以下形式:

其中，Ξns=［φ(‖xp－x1‖)，…，φ(‖xp－xns‖)］。將方程(6)代入式(7)和(9)可得:

其中:Ans= Ψnsθns－1;Bns= Ξnsθns－1。

1.2 MAPS方法局部形式推廣到全局形式

通過在適當位置添加零元素將式(11)中的矩陣Bns從局部形式推廣到全局形式。

假設(shè)N=100，ns=5，Ωp={xp，1，xp，2，xp，3，xp，4，xp，5}={x20，x23，x27，x54，x90}，5=，則式(11)可寫成如下形式:

則上式可寫成:Δu(xp)=B100(xp)。

由上述方程可見，在矩陣B100(xp)省略號位置添加零元素，成功地用替換局部的，使得方程從局部形式變成全局形式，構(gòu)成一個全局的稀疏矩陣，大大提高了計算的效率，成功降低了求解大規(guī)模問題時的困難。

1.3 用LMAPS法離散時空偏微分方程

其中:λ是一個給定的常數(shù);Δ是二階線性拉普拉斯算子;B是邊界微分算子。

簡而言之，H(x，u，uxuy)能用以下方程表示

這里a，b，c是給定系數(shù)函數(shù)。由方程(10)和(11)、離散方程(12)和(15)可得:

其中:n是時間步長的指數(shù);ns是在支撐域Ωi內(nèi)插值點xi的個數(shù)，i=1，…，N。

方程(16)、(17)中的系數(shù)aij，cij，bij可以通過方程(6)的系數(shù)矩陣求得。為了使得到的結(jié)果更加精確，通過在適當位置添加零元素可以將方程(16)、(17)中的系數(shù)矩陣從局部形式推廣到全局形式，時空偏微分方程寫成以下形式:

這里，A，B和C是N×N的稀疏矩陣，F(xiàn)(x)=(f(x1)，f(x2)，…，f(xn))T。通過已知的條件求解該稀疏矩陣，可求得每時間層結(jié)點的近似解。

1.4 Matern徑向基函數(shù)

本文采用Matern徑向基函數(shù):

其中:ΔΦ=r2K2(r);γ是歐拉常數(shù);K0(r)和K1(r)分別是0階和1階第二類貝塞爾函數(shù)［11］(Bessel function)。

1.5 時間層上龍格—庫塔(Rung-kutta)法離散形式

方程(19)是由第n時間層值推導第n+1時間層值，用4階龍格—庫塔方法計算:

這里的空間算子L是通過方程(16)右邊部分給出。u(1)，u(2)，u(3)是預測值，起到中間變量的作用。因為在每一個時間層(不包含初始層)內(nèi)部點的函數(shù)值是未知的，全部的空間導數(shù)能夠通過每個坐標方向上的網(wǎng)格點離散近似得到，較低時間層上的值能夠通過初始值求得，每一時間層上數(shù)值都可以通過相應的初始條件和邊界條件一步步求得。

2 數(shù)值算例

用MAPS方法離散算例1并與Kansa方法進行誤差比較，在時間層上采用Rung-kutta法進行求解。在算例2中用2種徑向基函數(shù)的LMAPS法求解N－S方程組。這2個算例比較最大絕對誤差(MAE)、最大相對誤差(RAE)、均方根誤差(RMSE)，其定義如下:

其中:N是測試點;?(xj，tn)是近似解;u(xj，tn)是精確解。為了提高計算效率和計算精度，本文采用Kdtree算法搜索局部區(qū)域的相鄰節(jié)點。

算例1 在單位正方形［－0.5，0.5］×［－0.5，0.5］內(nèi)考慮關(guān)于式(1)的一個瞬態(tài)擴散方程［12］:

初始條件:P(x，y，0)=1，(x，y)∈Ω∪?Ω。

狄利克雷邊界:P(x，y，t)=0，(x，y)∈?Ω 且 t＞0。

這里取λ=1。該方程描述了包括導電熱的流量和在靜止的流體擴散方面的許多物理學過程，也被認為類似于一滴染料在一杯水中傳播的一種擴散過程。本研究想表明在內(nèi)部時間t從0到1有一個跳躍條件，使該問題變得更具挑戰(zhàn)性。

現(xiàn)在設(shè)定時間步長Δt=10－4，這里選定形參c=80，問題(22)的精確解為:

表1 當c=80，ns=5，Δt=10－4，N=441時的4級Rung-Kutta方法的誤差比較

表2 當c=80，ns=5，Δt=10－4時LMAPS方法和局部Kansa方法的誤差比較

由表1可知:在c，ns，Δt，N相同的條件下應用四級Rung-Kutta方法求數(shù)值解時，隨著時間層的增多，數(shù)值解的誤差越來越小，數(shù)值解的精度越來越高。

從表2可知:① 當c，ns，Δt，t取值一致時，節(jié)點取得越細，劃分得越密，所獲得的誤差越小，數(shù)值解的精度就越高，但所消耗的時間也越多;② 在c，N，t，Δt，ns相同的情況下，LMAPS方法比局部Kansa方法產(chǎn)生的誤差略小一些，LMAPS方法求得的數(shù)值精度高，所消耗的時間略微少一點。

圖1 當t=0.1，N=961時采用LMAPS方法和Kansa方法求得的誤差

算例2 考慮應用LMAPS和local Kansa方法分析方腔內(nèi)的非定常不可壓Navier-Stokes方程組流體流動問題，其控制方程如下:

圖2 采用全局local Kansa方法和全局LMAPS法求解N－S方程組的誤差分布

圖3 當n=441，t=0.1時LMAPS關(guān)于形參c的誤差

表3 基于不同徑向基函數(shù)的LMAPS法求解N－S方程組所得的誤差對比

由圖2可見:這2種方法的全局化求解N－S方程組得到的近似解與精確解吻合得很好，計算精度也非常高，都是行之有效的方法。

由圖3左圖可見:隨著形參c的增大，φ的數(shù)值解精度會變得很高，但當形參c達到一定值時誤差趨于穩(wěn)定。由圖3右圖可見:網(wǎng)格劃分越細密則計算精度越高，且形參c的最優(yōu)值也會逐漸增大。由表3可知:用MQ函數(shù)插值的LMAPS方法求解N－S方程的誤差精度能達到10－6，很明顯Matern函數(shù)插值的LMAPS法求解的精度達不到那么高，但是它省去了選取形參的時間，誤差精度也非常高。

3 結(jié)束語

本文用LMAPS和local Kansa基于RBFs的無網(wǎng)格方法求解時空偏微分方程。算例1證明了在時間層上用4級Rung-kutta法求解是一種行之有效的方法，具有較高精度。在算例2中將基于MQ函數(shù)插值和基于Matern函數(shù)插值的LMAPS法所取得的數(shù)值結(jié)果進行了比較，結(jié)果表明:由基于Matern函數(shù)插值的LMAPS法得到的數(shù)值結(jié)果精度很高，而且不需要像MQ函數(shù)那樣不斷調(diào)試合適的形參c，這樣可以節(jié)約很多時間。

［1］lucy L B.A numerical approach to the testing of the fission hypothesis［J］.Astrophysical Jourral，1997(8):24.

［2］Belyschko T，Cu L.Element-free Galerkin methods［J］.Numerical Methods in Engineering，1994，37:56.

［3］liu W，Jun S，Zhang Y.Reproducing kernel particle methods［J］.International Journal for Numerical Methods in Fluids，1995，20:106－1081.

［4］Bellman R E，Kashef B G，Casti J.Differential quadrature a technique for the rapid solution of nonlinear partial differential equations［J］.Journal of computational physics ，1972，10:40－52.

［5］Yao G M，Chen C S，Joseph Kolibal.A localized approach for the method of approximate particular solutions［J］.Computers and Mathematics with Applications，2011，61:2376－2387.

［6］Chen C S ，F(xiàn)an C M，Wen P H，The method of particular solutions for solving certain partial differential equations［J］.Numerical Methods of Partial Differential Equations，2012，28:506－522.

［7］Vertnik R.Local collocation method for phase-change problems，Master’s thesis［D］.［S.l.］:University of Nova Gorica，2007.

［8］Kansa E J.Multiquadrics-a scattered data approximation scheme with applications to comutational fluid-dynamics.I.Surface approximations and partial derivative estimates［J］.Comput Math Appl，1990，19:127－45.

［9］Kansa E J.Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics.II.Solutions to parabolic，hyperbolic and elliptic partialdifferential equations［J］.Comput Math Appl，2000，39:37－123.

［10］Dehghan M.Finite difference procedures for solving a problem arising in modeling and design of certain optoelectronic devices［J］.Math compact simulation，2006，71:16－30.

［11］Alassar R S，Mavromatis H A.A New Integral Involving the Product of Bessel Functions and Associated Laguerre Polynomials［J］.Acta Applicandae Mathematicae，2008(2):263－267.

［12］Xueying Zhang，Huiqing Zhu，Lei-Hsin Kuo.A comparison study of the LMAPS method and the LDQ method for time-dependent problems［J］.Engineering Analysis with Boundary Elements，2013，37:1408－1415.

［13］Shu C，Yao Q，Yeo K S.Block-marching in time with DQ discretization:an efficient method for time-dependent problems［J］.Comput Methods Appl Mech Engrg，2002，191:4590－4592.