雙調(diào)和方程的有限積分方法

李書(shū)偉，徐定華，余　躍

(浙江理工大學(xué)理學(xué)院，杭州 310018)

?

雙調(diào)和方程的有限積分方法

李書(shū)偉，徐定華，余躍

(浙江理工大學(xué)理學(xué)院，杭州 310018)

利用有限積分法求解平面矩形區(qū)域雙調(diào)和方程邊值問(wèn)題。首先，對(duì)雙調(diào)和方程以及邊界條件分別進(jìn)行積分,得到一帶有任意函數(shù)的線性常微分方程組；其次，將積分產(chǎn)生的任意函數(shù)分別進(jìn)行插值估計(jì)，進(jìn)而轉(zhuǎn)化成為一可求解的線性代數(shù)方程組；最后，利用正則化方法求解奇異線性方程組，獲得近似解誤差估計(jì)。通過(guò)Matlab進(jìn)行數(shù)值模擬實(shí)驗(yàn)獲得數(shù)值結(jié)果，并進(jìn)行誤差分析。數(shù)值結(jié)果表明，與有限差分法、有限元法以及廣義有限差分法相比較，有限積分法具有更高精度。

雙調(diào)和方程；有限積分法；正則化方法；誤差估計(jì)；數(shù)值解

0　引　言

雙調(diào)和方程主要來(lái)源于流體力學(xué)和彈性力學(xué)，在工程技術(shù)方面應(yīng)用極為廣泛，相關(guān)理論和方法的研究一直是人們關(guān)注的熱點(diǎn)，有限差分法[1-2]和有限元法[3]是常用的求解雙調(diào)和方程的兩種數(shù)值解法。有限差分法的優(yōu)點(diǎn)是方法簡(jiǎn)單、計(jì)算量小，但是對(duì)邊界區(qū)域規(guī)則性的要求非常高，難以處理復(fù)雜網(wǎng)格。有限元法雖然可以針對(duì)任何邊界區(qū)域，但對(duì)于區(qū)域網(wǎng)格的劃分有嚴(yán)格的要求。有限體積法[4]和邊界元法[5]也是求解雙調(diào)和方程非常重要的數(shù)值解法。有限體積法雖然自身包含幾何信息、易處理復(fù)雜網(wǎng)格，但是計(jì)算比較復(fù)雜、不易提高精度。邊界元法最大的優(yōu)點(diǎn)是，可以使問(wèn)題的空間維數(shù)降低，從而使計(jì)算工作量及其所需計(jì)算機(jī)容量大大減小，但是該方法需要求出問(wèn)題的基本解，推導(dǎo)過(guò)程比較復(fù)雜。無(wú)網(wǎng)格法[6-8]是近年來(lái)新興的一種數(shù)值方法，它是針對(duì)任何邊界區(qū)域，并且只需在計(jì)算區(qū)域上選取若干節(jié)點(diǎn)及其徑向基函數(shù)，不需要?jiǎng)澐志W(wǎng)格。

目前，偏微分方程邊值問(wèn)題的求解方法大多為有限差分法、有限元法、有限體積法和邊界元法。迄今為止，對(duì)應(yīng)用有限積分法求解偏微分方程的研究工作還不多，在該方法理論上的研究很少。本文研究了雙調(diào)和方程邊值問(wèn)題的有限積分求解方法、數(shù)值算法、數(shù)值模擬以及該方法的收斂性分析，并通過(guò)模擬平面薄板受一均勻載荷時(shí)，研究靜態(tài)平衡下薄板與載荷接觸面上的位移分布與誤差分析。本文針對(duì)研究雙調(diào)和方程邊值問(wèn)題的求解方法，擬采用Wen等[9]提出的有限積分法對(duì)該問(wèn)題進(jìn)行數(shù)值求解，考慮了該方法的理論基礎(chǔ)——積分的離散與函數(shù)的插值表示，結(jié)合文獻(xiàn)[9]給出的利用有限積分法求解一維偏微分方程邊值問(wèn)題的解題思路，以求出雙調(diào)和方程的數(shù)值解，并從理論上對(duì)該方法的收斂性做了分析。

1　預(yù)備知識(shí)

1.1積分離散形式

下面考慮一維情形[9]?？紤]定義在[a,b]上的函數(shù)u(x)，定義域離散為a=x0

(1)

則一次積分可離散為：

k=0,1,…,N

(2)

式(2)用矩陣表示為：

U(1)=A(1)u

(3)

其中：

u=[u0,u1,…，uN]T.

同理，二次積分的離散形式為：

U(2)=A(2)u=(A(1))2u.

更一般地，m次積分的離散結(jié)果為:

U(m)=A(m)u=(A(1))mu.

a=x0

c=y0

為方便表述，把網(wǎng)格節(jié)點(diǎn)標(biāo)記為：

(x0,ys),(x1,ys),…,(xN1,ys),s=0,1,…,N2，

稱這樣的順序?yàn)闃?biāo)準(zhǔn)順序。

記Ux(x,y)和Uy(x,y)分別為對(duì)x和y的積分，則：

s=0,1,2,…,N2，

用矩陣表示為：

具體可寫為：

或簡(jiǎn)記為：

其中A(1)是由方程(3)求得的(N1+1)×(N1+1)維積分矩陣。

關(guān)于x的m次積分的離散形式為：

同理，關(guān)于y的積分的離散形式為：

為了將其轉(zhuǎn)化為標(biāo)準(zhǔn)順序，可進(jìn)行一個(gè)初等行交換，其對(duì)應(yīng)初等矩陣T=(Tmn)，滿足：

其中：i=0,1,…,N1;j=0,1,…,N2

1.2函數(shù)的插值表示

下面考慮線性插值[9]。對(duì)于1.1節(jié)給出的函數(shù)u(x)進(jìn)行線性插值，即使用梯形規(guī)則，得：

(4)

則：

用矩陣表示為：

U(1)=A(1)u.

對(duì)于二階線性插值有：

U(2)=A(2)u=(A(1))2u.

對(duì)于m階線性插值有：

U(m)=A(m)u=(A(1))mu.

=R(x)α+P(x)β

(5)

其中：Ri(x)是徑向基函數(shù)，αi是Ri(x)的系數(shù)，Pj(x)是多項(xiàng)式基函數(shù)，βj是Pj(x)的系數(shù)，徑向插值的一個(gè)附加條件為：

(6)

聯(lián)立方程(5)和(6)得：

(7)

則方程(7)的一階積分插值為：

因此，其一階的有限積分矩陣為：

m階積分的有限積分矩陣為：

A(m)=Am.

2　雙調(diào)和方程的有限積分法

考慮如下雙調(diào)和方程邊值問(wèn)題：

(8)

對(duì)式(8)中方程的x和y分別進(jìn)行四次積分有：

x2f1(y)+xf2(y)+f3(y)+y3g0(x)+

y2g1(x)+yg2(x)+g3(x)=

dξ4dξ3dξ2dξ1dη4dη3dη2dη1

(9)

根據(jù)有限積分，方程(9)可離散為：

+X3Φyρ(0)+X2Φyρ(1)+XΦyρ(2)+Φyρ(3)

+Y3ΦxΛ(0)+Y2ΦxΛ(1)+YΦxΛ(2)+ΦxΛ(3)

(10)

則方程(10)可以根據(jù)邊界條件確定唯一解，并且其導(dǎo)數(shù)邊界條件可以表示為：

(11)

3　誤差分析

考慮方程(10)，令：

則得到如下線性方程組：

簡(jiǎn)記為

(12)

(13)

為避免因系數(shù)矩陣K的奇異性而導(dǎo)致錯(cuò)誤的數(shù)值解，采用正則化方法解方程組(13)，并對(duì)有限積分法計(jì)算結(jié)果做出誤差估計(jì)。

證明可參見(jiàn)文獻(xiàn)[12]中定理3.1。

特別的，當(dāng)δ→0時(shí)，則有如下誤差估計(jì)：

這里‖·‖表示L2-范數(shù)。

對(duì)K進(jìn)行奇異值分解得：

則，有正則解：

令Rα=(αI+KTK)-1KT，由文獻(xiàn)[13]可知：

由此可得：

根據(jù)引理1得：

又因?yàn)椋?/p>

則：

特別地，令δ→0且α=hμ-3/E(μ>3)，則誤差估計(jì)為：

4　數(shù)值試驗(yàn)

本節(jié)用有限積分法求解一個(gè)雙調(diào)和方程，其數(shù)值結(jié)果與用有限差分法和有限元方法得出的數(shù)值結(jié)果相比較，通過(guò)實(shí)例驗(yàn)證方法是否有效、適用。

算例1考慮定解問(wèn)題(8)，令：

h(x,y)=8[3x2(1-x)2+3y2(1-y)2+(6x2-6x+1)(6y2-6y+1)],f(x,y)=g(x,y)=0，本文用有限積分法解該問(wèn)題，并將計(jì)算得的數(shù)值結(jié)果與文獻(xiàn)[1]中的十三點(diǎn)格式的有限差分法和廣義有限差分法以及文獻(xiàn)[3]中的有限元方法得出的數(shù)值結(jié)果進(jìn)行對(duì)比。數(shù)值結(jié)果及精確解見(jiàn)表1，數(shù)值誤差見(jiàn)表2。

表1　數(shù)值結(jié)果

表2　誤差結(jié)果

由表1、表2易見(jiàn)該方法具有如下優(yōu)點(diǎn)：該方法是一種邊界型無(wú)網(wǎng)格法，具有不需要?jiǎng)澐志W(wǎng)格的優(yōu)勢(shì)；與其它數(shù)值解法相比較，該方法具有更高的計(jì)算精度；該方法的精度與節(jié)點(diǎn)數(shù)量的多少無(wú)關(guān)。因此有限積分法是處理雙調(diào)和邊值問(wèn)題的一種有效數(shù)值解法。

算例2在彈性力學(xué)[14]里，對(duì)于一受一般載荷的平面薄板，由薄板的小撓度彎曲理論可得該薄板彈性曲面的微分方程為：

Δ2u(x,y)=q/D,(x,y)∈Ω，

作為算例，本文考慮一受均勻載荷q=10 N的平面薄板，板邊長(zhǎng)度a=1 m，寬度b=1 m，厚度δ=0.01 m，彈性模量E=1.372×108Pa，泊松比μ= 0.15，取其邊界條件為固支邊界條件，即：

u(x,y)=0,(x,y)∈Γ，

將計(jì)算所得的數(shù)值結(jié)果與納維葉精確解[13]：

做比較，數(shù)值誤差結(jié)果見(jiàn)表3。

表3　數(shù)值誤差結(jié)果

由表3易見(jiàn)，對(duì)于薄板小撓度彎曲問(wèn)題，有限積分法計(jì)算所得的數(shù)值解很好地接近納維葉法計(jì)算所得的精確解，由此可見(jiàn)，有限積分法對(duì)于求解薄板小撓度彎曲問(wèn)題是有效的數(shù)值算法。

5　結(jié)　論

本文利用有限積分法求解了平面矩形區(qū)域的雙調(diào)和方程的邊值問(wèn)題，與其他方法相比，該方法優(yōu)點(diǎn)是無(wú)需劃分網(wǎng)格并且成功地避免了有限差分方法中的截?cái)嗾`差問(wèn)題，其計(jì)算結(jié)果比其他數(shù)值解法具有更高的精度。

雖然本文研究的是矩形區(qū)域的線性經(jīng)典方程邊值問(wèn)題，但是有限積分法同樣適用于非規(guī)則區(qū)域的非線性分?jǐn)?shù)階方程邊值問(wèn)題。對(duì)于非規(guī)則區(qū)域問(wèn)題,只需要按照區(qū)域范圍對(duì)函數(shù)每個(gè)積分點(diǎn)的積分區(qū)間稍作調(diào)整便可求解；對(duì)于非線性方程邊值問(wèn)題，只需按照文中所示的計(jì)算過(guò)程，將其轉(zhuǎn)化成非線性代數(shù)方程組，通過(guò)迭代解法便可求得數(shù)值解；對(duì)于分?jǐn)?shù)階方程邊值問(wèn)題，由文獻(xiàn)[9]和文獻(xiàn)[15]可知文中所述的有限積分方法同樣適用，只是分?jǐn)?shù)階方程積分所得的代數(shù)方程組有所差異。

[1] 李永海. 解雙調(diào)和方程的一種混合廣義差分法[J]. 吉林大學(xué)自然科學(xué)學(xué)報(bào), 1993(3): 19-30.

[2] CHEN G, LI Z, LIN P. A fast finite difference method for biharmonic equations on irregular domains and its application to an incompressible Stokes flow[J]. Advances in Computational Mathematics, 2008, 29(2): 113-133.

[3] WANG C M, WANG J P. An efficient numerical scheme for the biharmonic equation by weak Galerkin finite element methods on polygonal or polyhedral meshes[J]. Computers and Mathematics with Applications, 2014, 68: 2314-2330.

[4]WANG T. A mixed finite volume element method based on rectangular mesh for biharmonic equations [J]. Journal of Computational and Applied Mathematics, 2004, 172: 117-130.

[5] 溫希重, 李榮華. 求解平面雙調(diào)和方程邊界元法的誤差分析[J]. 高等學(xué)校計(jì)算數(shù)學(xué)學(xué)報(bào), 1988, 4: 298-310.

[6] LI M, JIANG T S, HON Y C. A meshless method based on RBFs method for nonhomogeneous backward heat conduction problem [J]. Engineering Analysis with Boundary Elements, 2010, 34: 785-792.

[7] LI M, CHEN C S, HON Y C. A meshless method for solving nonhomogeneous Cauchy problems[J]. Engineering Analysis with Boundary Elements, 2011, 35: 499-506.

[8] HON Y C, YANG Z. Meshless collocation method by Delta-shaped basis functions for default barrier model[J]. Engineering analysis with boundary elements, 2009, 33(7): 951-958.

[9] WEN P H, HON Y C, LI M, et al. Finite integration method for partial differential equations[J]. Applied Mathematical Modelling, 2013, 37(24): 10092-10106.

[10] DUAN Y, ZHENG Y M, CENG P P. Convergence estimate of the RBF-based meshless method for initial-boundary value problem of wave equations[J]. Engineering Analysis with Boundary Elements, 2012, 36(3): 303-309.

[11] HACKBUSCH W. Elliptic Differential Equations:Theory and Numerical Treatment[M]. Verlag Berlin Heidelberg: Springer, 1992: 103-109.

[12] NARCOWICH F J. Recent developments in error estimates for scattered-data interpolation via radial basis functions [J]. Numerical Algorithms, 2005, 39(1/3): 307-315.

[13] 劉繼軍, 不適定問(wèn)題的正則化方法及應(yīng)用[M]. 北京：科學(xué)出版社, 2005: 46-55.

[14] 徐芝綸. 彈性力學(xué): 下冊(cè)[M]. 四版, 北京：高等教育出版社, 2006: 1-41.

[15] HEYDARI M H, HOOSHMANDASL M R, MAALEK GHAINI F M. An efficient computational method for solving fractional biharmonic equation[J]. Computers and Mathematics with Applications, 2014, 68(3) :269-287.

(責(zé)任編輯： 康鋒)

Finite Integration Method for Biharmonic Equations

LIShuwei,XUDinghua,YUYue

(School of Science, Zhejiang Sci-Tech University, Hangzhou 310018,China)

In this paper, a finite integration method is applied to deal with boundary value problem of biharmonic equations on rectangle domains. Firstly, biharmonic equation and boundary conditions are integrated to gain a system of linear ordinary differential equations equipped with some arbitrary functions. Secondly, interpolation estimation is conducted for arbitrary function generated from integration. The boundary value problem is thus transformed to the system of linear algebraic equations which can be solved. Finally, regularization method is used to solve singular system of linear equations and to gain the error estimate of approximate solution. Numerical simulation experiment is carried out by Matlab software, and error analysis is conducted. Compared with finite difference method, finite element method and generalized finite difference, the results indicate that finite integration method presents higher precision.

biharmonic equations; finite integration method; regularization method; error estimate; numerical solution

10.3969/j.issn.1673-3851.2016.01.022

2015-04-17

國(guó)家自然科學(xué)基金項(xiàng)目(11071221，11471287)

李書(shū)偉(1989-)，男，山東德州人，碩士研究生，主要從事反問(wèn)題理論及應(yīng)用研究。

徐定華，E-mail：dhxu6708@zstu.edu.cn

O242.2, O241.8, O343.1

A

1673- 3851 (2016) 01- 0133- 07 引用頁(yè)碼： 010802