楊海婷， 張學(xué)瑩

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

近年來，無網(wǎng)格方法在計(jì)算偏微分方程數(shù)值解方面發(fā)展迅速，當(dāng)遇到網(wǎng)格劃分困難、網(wǎng)格突變以及網(wǎng)格移動(dòng)等問題時(shí)，無網(wǎng)格方法具有明顯的優(yōu)勢(shì)。

全局配點(diǎn)型 RBFs方法［1－2］是一種真正意義上的無網(wǎng)格方法，它成功地克服了上述困難，但是采用它求解偏微分方程時(shí)隨著插值節(jié)點(diǎn)數(shù)的增加，系數(shù)矩陣的條件數(shù)增大，容易導(dǎo)致系數(shù)陣奇異，滿秩甚至病態(tài)，這給處理大規(guī)模的問題帶來不便。為了改進(jìn)這些問題，文獻(xiàn)［3－5］采用局部徑向基函數(shù)節(jié)點(diǎn)插值法減少了矩陣的病態(tài)性。之后，CHEN C S等人采用LMQ思想，在近似特別解(Method of Approximate Particular Solution，MAPS)［6］的基礎(chǔ)上提出了局部近似特別解法(簡(jiǎn)稱 LMAPS)［7］。LMAPS方法具有MAPS和RBFs的優(yōu)點(diǎn)，它不僅規(guī)避了全局矩陣的病態(tài)和滿秩，而且可以用較少的節(jié)點(diǎn)數(shù)獲得更高的精度。目前，LMAPS方法已經(jīng)取得了很大的發(fā)展并廣泛應(yīng)用于物理領(lǐng)域。

文中介紹了MQ函數(shù)中形參c的優(yōu)化方法，分別選取MQ函數(shù)和Matern函數(shù)作為L(zhǎng)MAPS方法中的徑向基函數(shù)，并將這兩種函數(shù)應(yīng)用到規(guī)則區(qū)域和不規(guī)則區(qū)域上的偏微分方程進(jìn)行數(shù)值求解，同時(shí)給出相應(yīng)的誤差比較分析。

1 LMAPS方法的基本理論

1.1 控制方程

考慮偏微分方程

其中，λ，a，b，c是給定的常數(shù)，Δ是二階線性Laplace算子，B 是邊界算子，f(x，y)，g(x，y)和 u0是給定的函數(shù)。

1.2 LM APS方法的基本理論

這里‖·‖是歐幾里得范數(shù)，φ是徑向基函數(shù)。

此處

并且

由于局部支撐域內(nèi)各插值節(jié)點(diǎn)相異，從而矩陣Φns可逆，故有

對(duì)于?xp∈Ω，將方程(1)在空間上離散得到

其中

并且

通過在適當(dāng)位置添加零元素將向量Ψns延拓成n維向量 Ψn［7］，則有

1.3 MQ函數(shù)和M atern函數(shù)

局部近似特別解方法是一種基于RBFs的無網(wǎng)格方法，用它求解偏微分方程得到的數(shù)值解的精度很大程度上取決于基函數(shù)。MQ函數(shù)具有收斂速度快、插值精度高的特點(diǎn)，而相比MQ函數(shù)，Matern函數(shù)插值的LMAPS方法則省去了調(diào)節(jié)形參c的時(shí)間，在降低計(jì)算量的基礎(chǔ)上提高了計(jì)算效率。文中選取MQ函數(shù)和Matern函數(shù)作為L(zhǎng)MAPS方法中的基函數(shù)，并根據(jù)數(shù)值結(jié)果進(jìn)行誤差比較分析。這里采用MQ函數(shù)的一般形式:

其中，x=(x，y)表示在物理空間上支撐域內(nèi)點(diǎn)的坐標(biāo)，c是形狀參數(shù)，對(duì)方程(5)求積分［8-9］可得

其中，r= ‖x－xj‖。

如果選取Matern函數(shù)作為基函數(shù)有

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

1.4 形參c的優(yōu)化方法

由于形參c的選取直接影響基于MQ函數(shù)插值的LMAPS方法的近似精度，傳統(tǒng)意義下，通過不斷嘗試不同的形參c［10］或者一些特殊的方法，比如說，Contour Pade’算法［11］來規(guī)避系數(shù)陣的病態(tài)或滿秩從而得到一個(gè)合適的形參。但是這些方法都增加了計(jì)算成本，降低了計(jì)算效率并且給處理大規(guī)模問題帶來諸多不便。

為了尋找最優(yōu)形參c，并且平衡用MQ函數(shù)作為基函數(shù)的LMAPS方法求解PDEs得到的數(shù)值解的精確性和穩(wěn)定性，這里采用以分析交叉驗(yàn)證(Leave－One－Out Cross Validation，LOOCV)［12］方法為理論基礎(chǔ)的Rippa LOOCV算法來解決離散數(shù)據(jù)插值問題中形參c的優(yōu)化以及局部支撐域內(nèi)ns的選取問題［12-13］。

2 數(shù)值算例

為了比較MQ函數(shù)和Matern函數(shù)在LMAPS方法中的近似精度，采用基于這兩種函數(shù)的LMAPS方法求解不規(guī)則區(qū)域和規(guī)則區(qū)域上的偏微分方程，并給出誤差分析。求解過程中選取kd－tree算法進(jìn)行局部區(qū)域相鄰節(jié)點(diǎn)的搜索，并分別通過Rippa LOOCV算法以及不斷調(diào)試形參c來提高計(jì)算的近似精度和計(jì)算效率。最大絕對(duì)誤差(MAE)、最大相對(duì)誤差(MRE)、均方根誤差(RMSE)定義如下:

其中，n是總節(jié)點(diǎn)數(shù)，k是時(shí)間層，^u(xj，tk)和 u(xj，tk)分別表示k時(shí)間層上的數(shù)值解和精確解。

這里選取常見于靜電學(xué)、機(jī)械工程等物理領(lǐng)域的Possion方程以及被用來解決激波、淺水波、交通流動(dòng)力學(xué)等物理問題的Burgers′方程作為數(shù)值算例來驗(yàn)證上述兩種函數(shù)的有效性。

算例1 考慮不規(guī)則區(qū)域下不規(guī)則點(diǎn)分布的二維Possion方程其邊界是一個(gè)不規(guī)則區(qū)域，計(jì)算區(qū)域如圖1所示。

圖1 n=4 000，nb=200時(shí)不規(guī)則區(qū)域下節(jié)點(diǎn)不均勻分布的Possion問題的計(jì)算區(qū)域Fig.1 Computational domain and distribution of the interior and boundary points with n=4 000，nb=200

圖2 是總節(jié)點(diǎn)數(shù)n=4 000，邊界點(diǎn)數(shù)nb=200時(shí)，優(yōu)化形參c后的MQ函數(shù)和Matern函數(shù)作為基函數(shù)得到的相對(duì)誤差對(duì)比圖?？梢钥闯?，這兩種函數(shù)的計(jì)算精度都比較滿意，在整個(gè)不規(guī)則計(jì)算區(qū)域上得到的誤差基本一致，并且都是平滑下降的。相比MQ函數(shù)，Matern函數(shù)在不規(guī)則邊界各角處的誤差出現(xiàn)了輕微的波動(dòng)現(xiàn)象。

圖2 n=4 000，nb=200時(shí)MQ和Matern函數(shù)得到的相對(duì)誤差比較Fig.2 Profile of RAE using MQ and Matern with n=4 000，nb=200

表1為nb=200，n不同時(shí)，MQ和Matern函數(shù)得到的誤差結(jié)果。一般而言，n越多，網(wǎng)格劃分越密，所得的數(shù)值解精度越高，但實(shí)際并非如此。與MQ函數(shù)相反，隨著n的增大，Matern函數(shù)取得的誤差不斷增大，當(dāng)增大到某個(gè)特定值時(shí)，MQ函數(shù)得到的誤差產(chǎn)生了較大的波動(dòng)并出現(xiàn)了增大趨勢(shì)。這說明誤差不僅與總節(jié)點(diǎn)數(shù)目n有關(guān)，還與節(jié)點(diǎn)分布有著密切的關(guān)系。

總的來說，在處理不規(guī)則點(diǎn)不均勻分布問題時(shí)，優(yōu)化形參c后的MQ函數(shù)得到的數(shù)值結(jié)果相比LMAPS－Matern方法具有更高的計(jì)算精度和計(jì)算效率。但是在總節(jié)點(diǎn)數(shù)n增大到一定程度條件下，此時(shí)從數(shù)值解的穩(wěn)定性上來看，Matern函數(shù)相比MQ函數(shù)更有優(yōu)勢(shì)。

表1 MQ和Matern函數(shù)在不同n時(shí)數(shù)值結(jié)果的誤差對(duì)比Tab.1 Errors and com putational time versus nusing MQ and Matern

算例2 考慮二維不定常Burgers′方程組:

其中，Ω ∪ ?Ω = ［0，1］2，0 ＜ t＜ + ∞，方程的精確解［14］ 為

圖3，4 是在 n=441，Re=100，Δt=10－3的條件下，t分別取t=0，t=1，t=2時(shí)的數(shù)值解圖像。由圖可見，MQ函數(shù)和Matern函數(shù)求得的數(shù)值結(jié)果基本一致，隨著時(shí)間的增大，數(shù)值解的梯度變大，解的不連續(xù)性增強(qiáng)。此外，與v相反，得到的u的數(shù)值解的不連續(xù)區(qū)域在時(shí)間t的推動(dòng)下逐漸向坐標(biāo)中心移動(dòng)。

圖3 n=441，Re=100，c=50時(shí)MQ函數(shù)在不同時(shí)刻得到的數(shù)值結(jié)果Fig.3 uand v distributions obtained by MQ with n=441，Re=100，c=50 at different time

圖4 n=441，Re=100時(shí)M atern函數(shù)在不同時(shí)刻得到的數(shù)值結(jié)果Fig.4 uand v distributions obtained by Matern with n=441，Re=100 at different time

表 2，3 是 Δt=10－4，t=0.1 時(shí)，兩種函數(shù)得到的誤差。數(shù)值結(jié)果表明它們都達(dá)到了相似的近似精度，并且隨著Re的增大，誤差均明顯增大，當(dāng)Re增大到某一特定值時(shí)，相比MQ函數(shù)，Matern函數(shù)得到了更精確的結(jié)果。同時(shí)隨著總節(jié)點(diǎn)數(shù)的不斷增多， Matern函數(shù)的優(yōu)勢(shì)也更加明顯。

表2 t=0.1時(shí)MQ函數(shù)的誤差結(jié)果Tab.2 Obtained error by MQ for example 2 at t=0.1

表3 t=0.1時(shí)Matern函數(shù)的誤差結(jié)果Tab.3 Obtained error by Matern for example 2 at t=0.1

結(jié)合圖3，4，在圖5中可以看出，隨著Re的增大，在數(shù)值解梯度較大的區(qū)域，這兩種函數(shù)取得的絕對(duì)誤差相比光滑的區(qū)域要大得多，并且誤差圖像出現(xiàn)了不同層次的波動(dòng)。當(dāng)Re達(dá)到100時(shí)，Burgers′方程接近半雙曲型，此時(shí)用這兩種函數(shù)作為徑向基函數(shù)的LMAPS方法都不能得到穩(wěn)定的數(shù)值解;但是相比MQ函數(shù)，Matern函數(shù)得到了較高的精度以及較小的波動(dòng)?？偟膩碚f，無論從近似精度還是計(jì)算效率上來看，Matern函數(shù)對(duì)于求解規(guī)則區(qū)域上多個(gè)變量的半雙曲型偏微分方程更有優(yōu)勢(shì)。

圖5 n=441，t=0.2時(shí)MQ和M atern函數(shù)在不同下的絕對(duì)誤差比較Fig.5 Comparison of the error p rofiles of uobtained by MQ and Matern with n=441，t=0.2

3 結(jié)語(yǔ)

介紹了基于徑向基函數(shù)插值的無網(wǎng)格LMAPS方法，分別選取MQ函數(shù)和Matern函數(shù)作為基函數(shù)，并利用局部配點(diǎn)法構(gòu)造低階線性方程組進(jìn)行插值近似來避免全局矩陣的病態(tài)和滿秩。文中選取不規(guī)則區(qū)域上的Possion問題，采用優(yōu)化形參c后MQ函數(shù)與Matern函數(shù)進(jìn)行誤差比較，同時(shí)選取規(guī)則區(qū)域上的二維Burgers′方程作為算例，并給出相應(yīng)的誤差分析。

數(shù)值實(shí)驗(yàn)表明這兩種函數(shù)都取得了較高的精度，由于避免了形參c的選取，優(yōu)化形參c后MQ函數(shù)相比Matern函數(shù)得到更高的近似精度和計(jì)算效率。而無論是求解不規(guī)則區(qū)域點(diǎn)不均勻分布問題還是高維半雙曲型的偏微分方程，Matern函數(shù)得到的數(shù)值解都有著更強(qiáng)的穩(wěn)定性。

［1］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－145.

［2］Kansa E J.Multiquadrics－a scattered data approximation scheme with applications to comutational fluid－dynamics.II.solutions to parabolic，hyperbolic and elliptic partial differential equations［J］.Comput Math Appl，2000，39:123－137.

［3］Wendland H.Piecewise polynomial，positive definite and compactly supported radial basis functions ofminimal degree［J］.Adv Comput Math，1995，4(1):389－396.

［4］WU Z.Multivariate compactly supported positive definite radial functions［J］.Adv ComputMath，1995，4(1):283－292.

［5］Vertnik R，Sarler B.Meshless local radial basis function collocation method for convective－diffusive solid－liquid phase change problems［J］.International Journal of Numerical Methods for Heat and Fluid Flow，2006，16:617－640.

［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］YAO GM，CHEN CS，Kolibal J.A localized approach for themethod of approximate particular solutions［J］.ComputMath Appl，2011，61:2376－2387.

［8］SHU C，DING H，Yeo K S.Local radial basis function－based differential quadrature method and its application to solve twodimensional incompressible Navier－Stokes equations［J］.Computer Methods in Applied Mechanics and Engineering，2003，192:941－954.

［9］CHEN C S，F(xiàn)AN CM，WEN PH.Themethod of particular solutions for solving elliptic problems with variable coefficients［J］.International Journal of Computational Methods，2011，8:545－559.

［10］Kansa E J，Carlson R E.Improved accuracy of multiquadric interpolation using variable shape parameters［J］.Comput Math Appl，1992，24:99－120.

［11］Bengrt F，WrightG.Stable computation ofmultiquadric interpolants for all values of the shape parameter［J］.ComputMath Appl，2004，47:497－523.

［12］Rippa S.An algorithm for selecting a good value for the parameter c in radial basis function interpolation［J］.Adv ComputMath，1999，11:193－210.

［13］Bengrt F，Julia Z.The Runge phenomenon and spatially variable shape parameters in RBF interpolation［J］.ComputMath Appl，2007，54:379－398.

［14］Young D L，F(xiàn)AN C M，HU S P，et al.The Eulerian－Lagrangian method of fundamental solutions for two－dimensional unsteady Burgers equations［J］.Engineering Analysiswith Boundary Element，2008，32:395－412.

(責(zé)任編輯:楊 勇)