梁方程降階計(jì)算的重心插值配點(diǎn)法

徐子康，王兆清，*，孫浩森，李金

（1.山東建筑大學(xué) 工程力學(xué)研究所，山東 濟(jì)南250101；2.山東建筑大學(xué) 學(xué)報(bào)編輯部，山東濟(jì)南250101；3.山東建筑大學(xué)理學(xué)院，山東 濟(jì)南250101）

梁方程降階計(jì)算的重心插值配點(diǎn)法

徐子康1，王兆清1，*，孫浩森2，李金3

（1.山東建筑大學(xué) 工程力學(xué)研究所，山東 濟(jì)南250101；2.山東建筑大學(xué) 學(xué)報(bào)編輯部，山東濟(jì)南250101；3.山東建筑大學(xué)理學(xué)院，山東 濟(jì)南250101）

采用重心插值配點(diǎn)法求解梁方程時，隨著計(jì)算節(jié)點(diǎn)數(shù)量的持續(xù)增加，其計(jì)算精度將逐步下降。通過對降階計(jì)算重心插值配點(diǎn)法的研究，可為數(shù)值求解梁方程提供一種數(shù)值穩(wěn)定性好、計(jì)算精度高的新方法。文章基于重心Lagrange插值及其微分矩陣，推導(dǎo)了梁方程降階計(jì)算重心插值配點(diǎn)法的公式，并通過數(shù)值算例驗(yàn)證其有效性。結(jié)果表明:隨著計(jì)算節(jié)點(diǎn)數(shù)量的持續(xù)增加，降階法的計(jì)算精度仍保持在10-10～10-12范圍內(nèi)；求解兩端簡支的梁方程時，兩步降階法的計(jì)算精度高于一步降階法；直接法計(jì)算矩陣條件數(shù)與節(jié)點(diǎn)數(shù)的7次方是同階的，而一步降階法計(jì)算矩陣條件數(shù)與節(jié)點(diǎn)數(shù)的4次方是同階的，降階法可以有效地降低計(jì)算矩陣的條件數(shù)，提高計(jì)算精度；重心插值配點(diǎn)法采用矩陣—向量形式的計(jì)算公式，便于程序的編寫，提高了計(jì)算效率。

梁方程；降階法；重心Lagrange插值；配點(diǎn)法

0 引言

梁方程是工程技術(shù)領(lǐng)域中常見的四階常微分方程。數(shù)值分析四階常微分方程邊值問題，按照微分方程離散方式的不同可分為Galerkin法［1-2］和配點(diǎn)法［3-4］。Galerkin法是基于加權(quán)積分表達(dá)的弱形式離散方法，其典型代表是有限元法。有限元法是一種低階數(shù)值方法，要得到問題的高精度數(shù)值解，需要劃分稠密的單元，這極大降低了分析問題的效率。有限元需要采用Galerkin法將控制微分方程離散為代數(shù)方程，在形成剛度矩陣過程中，需要進(jìn)行大量的積分計(jì)算，而數(shù)值積分的精度將直接影響最終的計(jì)算結(jié)果。配點(diǎn)法是基于微分方程強(qiáng)形式的離散方法，也就是強(qiáng)迫微分方程在給定的離散節(jié)點(diǎn)上精確成立［5］。配點(diǎn)法是一種真正意義的無網(wǎng)格方法，其具有計(jì)算精度高、計(jì)算公式簡單、不需要數(shù)值積分和程序?qū)嵤┓奖愕忍攸c(diǎn)［6-8］。

利用配點(diǎn)法的諸多優(yōu)點(diǎn)，國內(nèi)外數(shù)值研究者采用不同的配點(diǎn)型方法求解四階邊值問題。Yogesh利用三次B樣條函數(shù)近似未知函數(shù)，將計(jì)算區(qū)間采用等距節(jié)點(diǎn)離散，提出求解四階邊值問題的三次B樣條配點(diǎn)法，其計(jì)算精度相比Siraj提出的五次非多項(xiàng)式樣條配點(diǎn)法高1～3數(shù)量級［9］。Malek采用配點(diǎn)法將四階常微分方程在計(jì)算區(qū)間上離散，利用擬譜法在每個子區(qū)域上近似未知函數(shù)，將提出的擬譜配點(diǎn)法應(yīng)用于求解一維四階線性微分方程、二維雙調(diào)和方程和二維非矩形域問題［10］。Muite使用實(shí)空間微分法、譜空間預(yù)處理法、實(shí)空間積分法和譜空間積分法4種不同的Chebyshev配點(diǎn)法求解四階邊值問題，驗(yàn)證了不同方法下的計(jì)算精度和數(shù)值穩(wěn)定性［11］。雖然學(xué)者基于配點(diǎn)法提出了多種求解梁方程的數(shù)值方法，但所提方法的計(jì)算公式繁瑣，不便于程序的編寫，因而計(jì)算精度和計(jì)算效率都比較低。

Lagrange插值在數(shù)值分析方面具有重要作用，但在等距節(jié)點(diǎn)插值過程中，當(dāng)節(jié)點(diǎn)數(shù)量較大時，Lagrange插值表現(xiàn)出極大的數(shù)值不穩(wěn)定性。采用Chebyshev或Legendre節(jié)點(diǎn)等特殊分布的插值節(jié)點(diǎn)可以顯著提高Lagrange插值的穩(wěn)定性和計(jì)算精度［12］。王兆清等將Lagrange插值公式改寫成重心公式形式，提出了基于重心Lagrange插值的配點(diǎn)型數(shù)值方法—重心插值配點(diǎn)法（BICM）［6］。該方法首先是采用重心Lagrange插值近似未知函數(shù)并建立未知函數(shù)各階導(dǎo)數(shù)在插值節(jié)點(diǎn)處的微分矩陣，然后利用微分矩陣直接離散梁方程和邊界條件，最后求解矩陣形式的離散方程組。BICM采用矩陣—向量形式的計(jì)算公式，使得計(jì)算程序的編寫非常方便，極大地提高了數(shù)值求解梁方程的計(jì)算精度和計(jì)算效率，現(xiàn)已廣泛應(yīng)用于求解微分方程的邊值問題［13-15］、初值 問 題［16］和 初邊 值 問題［17］。

采用重心插值配點(diǎn)法直接求解梁方程，隨著計(jì)算節(jié)點(diǎn)數(shù)的持續(xù)增加，計(jì)算矩陣的條件數(shù)隨之增大，其計(jì)算精度將逐步下降［5］。為了保證梁方程在不同計(jì)算節(jié)點(diǎn)數(shù)下的求解精度，文章在重心Lagrange插值配點(diǎn)法直接計(jì)算（直接法）的基礎(chǔ)上，針對不同的邊界條件，提出求解梁方程的一步降階重心插值配點(diǎn)法（一步降階法）和兩步降階重心插值配點(diǎn)法（兩步降階法）。兩步降階法的主要思想是先求出中間變量，再利用中間變量和未知函數(shù)的關(guān)系解出未知函數(shù)，實(shí)現(xiàn)兩步重心插值配點(diǎn)法的降階計(jì)算。一步降階法將四階梁控制微分方程表達(dá)為未知函數(shù)和中間變量的耦合微分方程組，利用重心插值配點(diǎn)法求解二階耦合微分方程組。

1 重心Lagrange插值及其微分矩陣

考慮定義在區(qū)間 ［a，b］上的節(jié)點(diǎn)a＝x1＜x2＜… ＜xN＝b，函數(shù)w（x）在節(jié)點(diǎn)處的函數(shù)值為wi＝w（xi），i＝1，2，…，N。函數(shù)w（x）的重心Lagrange插值基函數(shù)由式（1）［5］表示為

重心Lagrange插值基函數(shù) ξi（x）和重心插值權(quán)ωi分別由式（2）和（3）表示為

式中:k＝1，2，…，N；i＝1，2，…，N。

函數(shù)w（x）在節(jié)點(diǎn)處的m階導(dǎo)數(shù)可由式（4）表示為

2 梁方程的重心插值配點(diǎn)法計(jì)算公式

2.1 求解梁方程的直接法

Euler-Bernoulli梁控制微分方程由式（5）表示為

式中:EI為梁的抗彎剛度；P為軸線方向上的荷載；f（x）為橫向荷載。

兩端簡支的邊界條件由式（6）表示為

兩端固支的邊界條件由式（7）表示為

由重心Lagrange插值微分矩陣，式（5）的矩陣形式由式（8）表示為

式（8）可由式（9）簡記為

式中:L＝EID（4）+PD（2）；w＝［w1，w2，…，wN］T；f＝（f（x1），f（x2），…，f（xN））T。

邊界條件（6）或（7）采用重心插值公式（1）離散，則邊界條件（6）或（7）的離散矩陣表達(dá)式可由式（10）表示為

將式（9）的主方程和式（10）的離散邊界方程聯(lián)立組合為一個過約束方程組，由式（11）表示為

上述邊界條件的施加方法稱為附加法，降階法的邊界條件也采用該方法作類似處理。采用最小二乘法求解方程組（11），得到未知函數(shù)在計(jì)算節(jié)點(diǎn)上的函數(shù)值。

2.2 簡支邊界條件下求解梁方程的兩步降階法

引進(jìn)中間變量，中間變量由式（12）表示為

則式（5）、（6）表示的簡支梁邊值問題可分解為2個Poisson方程邊值問題，分別由式（13）和（14）表示為

式（13）和（14）可采用重心插值配點(diǎn)法順次求解，構(gòu)成梁方程的兩步降階計(jì)算。兩步降階法相當(dāng)于將原梁方程化為兩個一維Poisson方程求解。由于邊界條件必須是未知函數(shù)和中間變量的約束，兩步降階法僅適用于求解兩端簡支邊界條件下的梁方程。對于一般邊界條件下的梁方程，需采用一步降階法求解。

2.3 一般邊界條件下求解梁方程的一步降階法

金屬平衡管理是衡量冶煉企業(yè)工藝技術(shù)和管理水平的重要標(biāo)志，通過計(jì)算金屬回收率，反映金屬元素的回收利用狀況和損失方向[2]。

聯(lián)立式（5）和（12），所得的方程組由式（15）表示為

采用重心插值微分矩陣離散方程（15），離散矩陣形式由式（16）表示為

采用附加法施加邊界條件（10），得到的過約束代數(shù)方程組由式（17）表示為

一步降階法將梁方程表達(dá)為未知函數(shù)和中間變量的耦合微分方程組，按照矩陣乘法的運(yùn)算法則，通過增加分塊零矩陣非常方便地施加各類邊界條件，因而該方法不僅適用于求解兩端固支邊界條件下的梁方程，也可應(yīng)用于求解任意邊界條件（如簡支和自由邊界條件）下的梁方程。

3 數(shù)值算例

數(shù)值算例計(jì)算程序利用MATLAB編寫，采用與第二類Chebyshev節(jié)點(diǎn)同分布形式的計(jì)算節(jié)點(diǎn)。直接法和一步降階法的計(jì)算矩陣分別為L1和L2，相對誤差由式（18）定義為

3.1 兩端簡支的梁方程

均布荷載作用下的簡支梁，梁方程的邊值問題由式（19）［5］表示為

數(shù)值計(jì)算時，取梁的抗彎剛度EI＝106、l＝2、 q＝1000，一律采用國際制單位。直接法和兩種降階法計(jì)算的數(shù)據(jù)結(jié)果見表1。不同節(jié)點(diǎn)數(shù)計(jì)算的相對誤差對數(shù)變化曲線如圖1所示。

表1 直接法和降階法計(jì)算矩陣條件數(shù)與計(jì)算的相對誤差Er

圖1 不同節(jié)點(diǎn)數(shù)計(jì)算的相對誤差對數(shù)變化曲線圖

觀察圖1可以發(fā)現(xiàn)，兩步降階法的計(jì)算精度最高，當(dāng)節(jié)點(diǎn)數(shù)量在41個之內(nèi)，3種方法的計(jì)算結(jié)果都保持在高精度范圍內(nèi)，隨著節(jié)點(diǎn)數(shù)量的繼續(xù)增加，直接法的計(jì)算精度不斷下降，當(dāng)節(jié)點(diǎn)數(shù)量增加到201時，直接法的計(jì)算結(jié)果沒有意義而兩種降階法仍然能夠得到高精度解。直接法和一步降階法計(jì)算矩陣條件數(shù)與節(jié)點(diǎn)乘方的關(guān)系曲線如圖2所示。

圖2 直接法和一步降階法計(jì)算矩陣條件數(shù)與節(jié)點(diǎn)乘方的關(guān)系曲線圖

由圖2可以直觀的看出，隨著節(jié)點(diǎn)數(shù)量的增加，直接法和一步降階法計(jì)算矩陣的條件數(shù)不斷增加，直接法計(jì)算矩陣條件數(shù)增長速度最快且增長曲線與N7基本吻合，一步降階法計(jì)算矩陣條件數(shù)增長曲線與N4基本重合。計(jì)算矩陣的條件數(shù)越大，說明矩陣越病態(tài)，計(jì)算的精度就越低。圖2中直接法計(jì)算矩陣迅速增加的條件數(shù)導(dǎo)致圖1中直接法計(jì)算精度的持續(xù)下降，而一步降階法可以極大地降低計(jì)算矩陣條件數(shù)從而得到更高的精度和更好的數(shù)值穩(wěn)定性。

3.2 兩端固支的梁方程

兩端固支的梁方程邊值問題由式（21）表示為

式（21）的解析解由式（22）表示為

梁方程（21）的邊界條件為兩端固支，采用直接法和一步降階法進(jìn)行求解。不同節(jié)點(diǎn)數(shù)計(jì)算的相對誤差對數(shù)變化曲線如圖3所示。直接法和一步降階法計(jì)算的數(shù)據(jù)結(jié)果見表2。

圖3 不同數(shù)量節(jié)點(diǎn)計(jì)算的相對誤差對數(shù)變化曲線圖

表2 直接法和一步降階法計(jì)算矩陣條件數(shù)與計(jì)算的相對誤差Er

觀察圖3發(fā)現(xiàn)，采用11個計(jì)算節(jié)點(diǎn)，直接法和一步降階法即可得到問題的高精度解，當(dāng)計(jì)算節(jié)點(diǎn)大于71時，直接法的計(jì)算精度迅速下降，而一步降階法仍保持在高精度范圍內(nèi)。直接法和一步降階法計(jì)算矩陣條件數(shù)與節(jié)點(diǎn)乘方的關(guān)系曲線如圖4所示。

圖4 直接法和一步降階法計(jì)算矩陣條件數(shù)與節(jié)點(diǎn)乘方的關(guān)系曲線圖

圖4驗(yàn)證了直接法計(jì)算矩陣條件數(shù)與N7是同階的，而一步降階法計(jì)算矩陣條件數(shù)與N4是同階的，進(jìn)一步說明一步降階法可以極大降低計(jì)算矩陣條件數(shù)，從而提高計(jì)算精度。

3.3 非齊次邊界條件下的梁方程

梁方程邊值問題由式（23）［10］表示為

式（23）的解析解由式（24）表示為

直接法和一步降階法計(jì)算的數(shù)據(jù)結(jié)果見表3。不同數(shù)量節(jié)點(diǎn)計(jì)算的相對誤差對數(shù)變化曲線如圖5所示。

表3 直接法和一步降階法計(jì)算矩陣條件數(shù)與計(jì)算的相對誤差Er

圖5 不同數(shù)量節(jié)點(diǎn)計(jì)算的相對誤差對數(shù)變化曲線圖

觀察圖5可以發(fā)現(xiàn)，隨著節(jié)點(diǎn)數(shù)的增加，直接法計(jì)算的誤差曲線呈“V”字形，曲線下降段的原因主要是問題的解析解為三角振蕩函數(shù)，采用較多的計(jì)算節(jié)點(diǎn)才能得到問題的高精度解；曲線上升段的原因是直接法計(jì)算矩陣的條件數(shù)迅速增加，導(dǎo)致直接法的計(jì)算精度持續(xù)下降。直接法和一步降階法計(jì)算矩陣條件數(shù)與節(jié)點(diǎn)乘方的關(guān)系曲線如圖6所示。

圖6 直接法和一步降階法計(jì)算矩陣條件數(shù)與節(jié)點(diǎn)乘方的關(guān)系曲線圖

圖6中的關(guān)系曲線再次得出了與圖2和4相同的結(jié)論，進(jìn)一步驗(yàn)證了一步降階法可以明顯降低計(jì)算矩陣條件數(shù)，從而提高計(jì)算精度。

4 結(jié)論

通過上述研究表明:

（1）數(shù)值求解四階梁方程，隨著計(jì)算節(jié)點(diǎn)數(shù)量的持續(xù)增加，直接法的計(jì)算精度逐步下降，直至失去數(shù)值意義，而一步降階法和兩步降階法的計(jì)算精度仍保持在10-10～10-12。相比直接法，降階法具有更高的計(jì)算精度和更好的數(shù)值穩(wěn)定性。

（2）由于邊界條件必須是未知函數(shù)和中間變量的約束，兩步降階法僅適用于求解兩端簡支邊界條件下的梁方程。求解兩端簡支的梁方程時，兩步降階法的計(jì)算精度高于一步降階法。

（3）求解梁方程時，計(jì)算矩陣的條件數(shù)與節(jié)點(diǎn)數(shù)存在著明顯的數(shù)量關(guān)系:直接法計(jì)算矩陣條件數(shù)與N7是同階的，而一步降階法計(jì)算矩陣條件數(shù)與N4是同階的。這說明一步降階法能夠通過降低求導(dǎo)的階數(shù)，有效降低計(jì)算矩陣的條件數(shù)，從而明顯改善計(jì)算精度。

（4）采用重心Lagrange插值微分矩陣離散四階梁控制微分方程，使得控制方程的離散過程非常簡便和直接。矩陣—向量形式的計(jì)算公式，非常適用于MATLAB計(jì)算程序的編寫，極大地提高了數(shù)值求解的效率。重心插值配點(diǎn)法的計(jì)算公式簡單，程序?qū)嵤┓奖?，是一種高精度的無網(wǎng)格數(shù)值分析方法。

［1］ Arnold D.N.，Awanou G.，Qiu W.，etal..Mixed finite elements for elasticity on quadrilateral meshes［J］.Advances in Computational Mathematics，2015，41（3）:553-572.

［2］ Atluri S.N.，Liu H.T.，Han Z.D.，et al..Meshless local petrov-Galerkin（MLPG）mixed collocation method for elasticity problems［J］.Cmes-computer Modeling in Engineering＆Sciences，2006，4（3）:141-152.

［3］ Kruse R.，Nguyenthanh N.，De Lorenzis L.，etal..Isogeometric collocation for large deformation elasticity and frictional contact problems［J］.Computer Methods in Applied Mechanics and Engineering，2015，296（1）:73-112.

［4］ Stevens D.，Power H.，Cliffe K.A.，et al..A meshless local RBF collocationmethod using integral operators for linear elasticity ［J］.International Journal ofMechanical Sciences，2014，88（1）: 246-258.

［5］ 李樹忱，王兆清.高精度無網(wǎng)格重心插值配點(diǎn)法——算法、程序及工程應(yīng)用［M］.北京:科學(xué)出版社，2012.

［6］ 王兆清，唐炳濤，李樹忱，等.求解間斷邊值問題的重心插值單元配點(diǎn)法［J］.山東建筑大學(xué)學(xué)報(bào)，2009，24（2）:115-118.

［7］ 李樹忱，王兆清，袁超.極坐標(biāo)系下彈性問題的重心插值配點(diǎn)法［J］.中南大學(xué)學(xué)報(bào)（自然科學(xué)版），2013，44（5）:2031-2040.

［8］ 王兆清，莊美玲，姜劍.非線性MEMS微梁的重心有理插值迭代配點(diǎn)法分析［J］.固體力學(xué)學(xué)報(bào)，2015，36（5）:453-459.

［9］ Yogesh G.，Manoj K..Numerical Method for Solving Boundary Value Problem Arising in Deflection of Beams［J］.Canadian Journal on Computing in Mathematics，Natural Sciences，Engineering and Medicine，2011，2（7）:166-169.

［10］Malek A.，Phillips T.N..Pseudospectral collocation methods for fourth-order differential equations［J］.IMA Journal of Numerical Analysis，1994，15（4）:523-553.

［11］Muite B.K..A numerical comparison of Chebyshev methods for solving fourth order semilinear initial boundary value problems ［J］.Journal of Computational and Applied Mathematics，2010，234（2）:317-342.

［12］Weideman J.A.，Reddy S.C..A MATLAB differentiation matrix suite［J］.ACM Transactions on Mathematical Software，2000，26 （4）:465-519.

［12］王兆清，唐炳濤，李樹忱，等.求解間斷邊值問題的重心插值單元配點(diǎn)法［J］.山東建筑大學(xué)學(xué)報(bào)，2009，24（2）:115-118，123.

［14］王兆清，李淑萍，唐炳濤.圓環(huán)變形及屈曲的重心插值配點(diǎn)法分析［J］.機(jī)械強(qiáng)度，2009，31（2）:245-249.

［15］王兆清，綦甲帥，唐炳濤.奇異源項(xiàng)問題的重心插值數(shù)值解［J］.計(jì)算物理，2011，28（6）:883-888.

［15］李淑萍，王兆清，唐炳濤.重心插值配點(diǎn)法求解初值問題［J］.山東建筑大學(xué)學(xué)報(bào)，2007，22（6）:481-485，506.

［17］王兆清，馬燕，唐炳濤.梁動力學(xué)問題重心有理插值配點(diǎn)法［J］.振動與沖擊，2013，32（22）:41-46.

（學(xué)科責(zé)編：趙成龍）

Barycentric interpolation collocation method based on depression of order for solving beam equations

Xu Zikang1，Wang Zhaoqing1，*，Sun Haosen2，et al.
（1.Institute of Mechanics，Shandong Jianzhu University，Jinan 250101，China；2.Editorial Department of Journal of Shandong Jianzhu University，Jinan 250101，China）

When the beam equations are solved by barycentric interpolation collocation method，computational accuracy will decline gradually as the number of nodes increases.The research on barycentric interpolation collocationmethod based on depression of order，can provide the new method that has a good numerical stability and high computational accuracy for beam equations.Based on barycentric Lagrange interpolation and its differentialmatrices，the formula of barycentric interpolation collocationmethod based on depression of order is derived.Numerical examples are given to verify the effectiveness of the proposed method.The result shows that the computational accuracy of depression of ordermethod remains the range of10-10～10-12as the number of nodes increases.When the beam equations with simply supported ends are solved，the computational accuracy of the two-step depression of order method is higher than the one-step depression of order method.The condition number of the directmethod is close to 7th power of the number of nodes，and the condition number of one-step depression of ordermethod is close to 4th power of the number of nodes.The depression of ordermethod can effectively reduce the condition number of computingmatrix such that improves the computational accuracy.By applying the computational formula with matrix-vector form，the program is easy to write and the computational efficiency of barycentric interpolation collocationmethod can beimproved remarkably.

beam equations；depression of order method；barycentric Lagrange interpolation；collocation method

O302

A

1673-7644（2017）03-0245-06

2017-04-25

國家自然科學(xué)基金面上項(xiàng)目（51379113）；國家自然科學(xué)基金項(xiàng)目（11471195）；山東省自然科學(xué)基金重點(diǎn)項(xiàng)目（ZR2016JL006）

徐子康（1993-），男，在讀碩士，主要從事工程數(shù)值分析方法等方面的研究.E-mail:xuzikang1993＠163.com

*：王兆清（1965-），男，副教授，博士，主要從事工程數(shù)值分析方法等方面的研究.E-mail:sdjzuwang＠126.com