陳林芳，張學(xué)瑩

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

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

1 時(shí)空偏微分方程的局部MAPS形式

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

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

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

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

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

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

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

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

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

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

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

通過在適當(dāng)位置添加零元素將式(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)省略號(hào)位置添加零元素，成功地用替換局部的，使得方程從局部形式變成全局形式，構(gòu)成一個(gè)全局的稀疏矩陣，大大提高了計(jì)算的效率，成功降低了求解大規(guī)模問題時(shí)的困難。

1.3 用LMAPS法離散時(shí)空偏微分方程

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

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

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

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

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

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

1.4 Matern徑向基函數(shù)

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

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

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

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

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

2 數(shù)值算例

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

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

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

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

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

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

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

表1 當(dāng)c=80，ns=5，Δt=10－4，N=441時(shí)的4級(jí)Rung-Kutta方法的誤差比較

表2 當(dāng)c=80，ns=5，Δt=10－4時(shí)LMAPS方法和局部Kansa方法的誤差比較

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

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

圖1 當(dāng)t=0.1，N=961時(shí)采用LMAPS方法和Kansa方法求得的誤差

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

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

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

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

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

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

3 結(jié)束語

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

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