祖福興， 徐績青， 李巖汀

(1. 重慶交通大學(xué) 國家內(nèi)河航道整治工程技術(shù)研究中心，重慶 400074； 2. 重慶交通大學(xué) 水利水運工程教育部重點實驗室，重慶 400074)

0 引 言

自20世紀40年代以來，求解非線性問題的計算方法一直是非線性科學(xué)中的重要課題。由于各種解法對于非線性問題的求解不具有封閉性，致使各方法往往只能對弱非線性問題有效。對于強非線性問題的求解，要么是針對個案問題采取特殊的計算技術(shù)來進行，要么是采用近似的追蹤方法來進行，因此在理論上沒有完備的通用方法。

梁的大撓度彎曲問題屬于強非線性問題，且常為四階微分方程并包含非冪次項，求解難度大，采用通常的有限元方法很難獲得令人滿意的結(jié)果。近年來，國內(nèi)外學(xué)者相繼提出了小波分析方法[1]、微分求積法[2]、保辛-保能的數(shù)值積分法[3,4]等，這些方法雖然各有優(yōu)點，但對于多元(多于4個變元，即多維空間的問題)強非線性微分方程，還需要引進某種近似，例如攝動法[5]等。由于這些近似仍會產(chǎn)生一些問題，尚需繼續(xù)實踐探討[6]。

徑向基函數(shù)以空間距離為自變量，由于其形式簡單、各向同性的優(yōu)點，近年來得到了很快的發(fā)展，并被廣泛應(yīng)用于散亂數(shù)據(jù)處理、微分方程求解等領(lǐng)域。同時，一種將徑向基函數(shù)與配點法結(jié)合的無網(wǎng)格法[7]也受到了廣泛的關(guān)注，但這種強格式無網(wǎng)格法通常在邊界處容易產(chǎn)生較大的數(shù)值震蕩。目前，通常采用在邊界處加密配點的不均勻配點方案、PDECB方法[8]、調(diào)整邊界點的權(quán)重[9]來減小計算誤差。徐績青等[10]提出了直接應(yīng)用徑向基函數(shù)配點法離散時間區(qū)域的動力問題求解方法，通過構(gòu)建適當?shù)牟逯岛瘮?shù)并添加初始條件實現(xiàn)了減小邊界處數(shù)值震蕩的目的。這種方法已成功用于Bratu型強非線性方程[11]和非線性動力系統(tǒng)的數(shù)值求解[12]，并取得了不錯的求解效果。筆者在此基礎(chǔ)上提出一種適用于梁大撓度彎曲問題求解徑向基函數(shù)逼近思想結(jié)合加權(quán)余量配點的高精度計算方法。

1 徑向基函數(shù)法

徑向基函數(shù)(radial basis function，RBF)是一類以徑向距離為自變量的基函數(shù)，記為：φ=φ(Rj)，其中，自變量Rj=‖x-xj‖2，為任意點x到節(jié)點xj的歐氏距離，且Rj∈Rd。而φ(Rj)在本質(zhì)上是一個定義在[0，+∞)上的一元函數(shù)，其簡單的形式為存儲和計算帶來了很大的優(yōu)勢，這個特點對于多元問題，如：板殼、空間結(jié)構(gòu)，甚至是高維問題，相比其他類型的基函數(shù)，例其具有明顯的優(yōu)勢。徑向基函數(shù)插值的收斂性非常好，不僅函數(shù)本身有較高的收斂階，其各階導(dǎo)數(shù)也有很好的收斂速度[13]。

根據(jù)作用范圍，徑向基函數(shù)可分為2大類：①全局徑向基函數(shù)(globally supported RBF)，這類RBF大多是在一定的數(shù)學(xué)和物理背景下提出的，因此很受歡迎，但它們在計算過程中經(jīng)常會產(chǎn)生病態(tài)稠密矩陣，且應(yīng)用全局徑向基函數(shù)的求解精度嚴重依賴于系數(shù)的選?。虎诰o支撐徑向基函數(shù) (compactly supported RBF)[14,15]，它們在計算中通常形成帶狀稀疏陣。與其他類型的緊支徑向基函數(shù)相比，由吳宗敏[13,15]教授所構(gòu)造的局部緊支撐徑向基函數(shù)，國際上稱為Wu函數(shù)，是嚴格正定的。筆者應(yīng)用Wu函數(shù)作為插值基函數(shù)，式(1)～式(3)這類徑向基函數(shù)按照在給定維數(shù)空間正定和連續(xù)性條件進行選用，具體原則及數(shù)學(xué)性質(zhì)見文獻[13，15-16]：

(1)

(2)

(3)

式中：r=Rj/Rmax, j(Rmax, j為定義在節(jié)點xj處的支撐域半徑)；(1-r)+定義為：

1.1 RBF方法基本過程

徑向基函數(shù)求解的一階非線性微分方程如式(4)[10,11]：

(4)

將求解區(qū)域Ω用n個節(jié)點xj(j=1，2，…，n)離散，函數(shù)y(x)在區(qū)域Ω里的近似解可用一個基本的插值函數(shù)表示，如式(5)：

(5)

式中：αj為各離散點的權(quán)重系數(shù)；α為權(quán)重系數(shù)向量，α=[α1，α2，…，αn]T；Φ(x)為插值基函數(shù)向量，Φ(x)=[Φ1(x)，Φ2(x)，…，Φn(x)]T。

令近似函數(shù)yh(x)在各個節(jié)點xj處的取值與函數(shù)y(x)在相同節(jié)點處的真實值yj相等，即yh(xj)=yj，則可以得到n維線性方程組：

Aα=y

(6)

y=[y1,y2,…,yn]T

由式(6)得到權(quán)重系數(shù)α=A-1y，代入式(5)可得：

yh(x)=ΦT(x)A-1y=N(x)y

(7)

令N(x)=ΦT(x)A-1，由于yh(x)是解析表達式，故可求得一階導(dǎo)數(shù)為

y′(x)=N′(x)y

(8)

式中：N′(x)為時間特征函數(shù)(類似于形函數(shù))。

可見，N(x)起到了形函數(shù)的作用。

將式(7)、式(8)代入方程(4)，則在求解區(qū)域Ω內(nèi)的原平衡方程，被離散成n個非線性代數(shù)方程的方程組(9)：

f(y,N′(x)y)=0

(9)

通常采用迭代算法來求解該非線性方程組，因此，可直接替換或添加邊值條件所對應(yīng)的方程，如將方程組(9)中x=a和x=b對應(yīng)的第1個和第n個方程替換為y1=a，yn=b，解該方程組即可求得y=[y1，y2，…，yn]T。

但是數(shù)值實驗顯示，這種徑向基函數(shù)方法的計算結(jié)果越靠近邊界誤差越大，最大誤差通常出現(xiàn)在邊界附近區(qū)域[17]。

1.2 誤差改進

針對微分方程(4)構(gòu)造插值函數(shù)如式(10)：

(10)

式(10)相比式(5)增加了邊界點上的一階導(dǎo)數(shù)項。同樣，將各配點坐標代入式(10)，可得到n維線性方程組，此時系數(shù)矩陣A不再是方陣，A-1是對應(yīng)的廣義逆矩陣。按照與式(7)～式(9)同樣的步驟，將方程(4)離散為一組非線性代數(shù)方程組，并代入邊界條件。為了使邊界點上不僅滿足邊界條件同時也滿足微分方程，考慮增加一階導(dǎo)數(shù)的邊界條件：

該邊界條件可由控制方程(5)求解。在離散過程中表達為

求解該n+ 2維的代數(shù)方程組可得微分方程的解。

為了克服經(jīng)典的牛頓法對初始點敏感和收斂性不佳等問題，筆者建議可聯(lián)合運用高斯-牛頓迭代法[18]和生物遺傳算法[19]來求解非線性代數(shù)方程組。

2 算法實施

以3.1的算例為例，本文方法在強非線性梁的微分方程邊值問題中的算法步驟如下：

1)建立微分方程，如式(11)。

2)確定求解區(qū)域Ω和配點距離Δx，以及配點數(shù)目n。

3)選擇合適的緊支撐正定徑向基函數(shù)φ(r)：

式中：xj為各配點處的位置；x為任意位置；Rmax, j為各配點的支撐域半徑，建議以配點位置為中心，覆蓋整個計算區(qū)間。

4)根據(jù)各配點位置計算n階特征矩陣A和其逆矩陣B=A-1。

5)任意位置x處的撓度用徑向基函數(shù)表示：

ωh(x)=ΦT(x)Bω=N(x)ω

求ωh(x)的各階導(dǎo)數(shù)。若為泛函表達式，如3.2的算例，則根據(jù)任意位置x處撓度的徑向基函數(shù)表達式：ωh(x)=N(x)ω，采用高斯-勒讓德積分離散原表達式，從而將泛函表達式轉(zhuǎn)化為代數(shù)表達式。

6)將配點位置xj處撓度的各階導(dǎo)數(shù)的徑向基函數(shù)表達式、轉(zhuǎn)化后的泛函表達式和外荷載的值代入微分方程，建立非線性方程組。在此過程中，為了減小數(shù)值振蕩，取得滿意的計算結(jié)果，邊界點通常還需要增加附加方程，以滿足相應(yīng)的導(dǎo)數(shù)邊界條件和邊界點處的微分方程，帶入邊界條件和附加的邊界條件。

7)求解方程組，得到各配點位置的撓度值。

3 算例分析

3.1 非線性梁的彎曲

求在彎矩

和梁一端作用軸向集中力P=100 N的聯(lián)合作用下，簡支梁的撓曲變形 。簡支梁的長度l=1 m，抗彎剛度EI=50 N·m2。若不考慮軸向變形，撓曲微分方程可以寫為：

(11)

該非線性邊值問題(11)有解析解(單位：m)：

可以看出，由于方程(11)中含有非冪次非線性項，難以處理。運用徑向基函數(shù)配點法進行求解，若采用式(5)作為插值函數(shù)，將產(chǎn)生較大的數(shù)值振蕩。針對問題(11)求解目標是撓度的二階導(dǎo)數(shù)，筆者借鑒彈塑性靜力學(xué)的處理方法，提出撓度及撓度的二階變率聯(lián)合插值的徑向基函數(shù)表達式：

(12)

為避免式(12)中徑向基函數(shù)的導(dǎo)數(shù)使得系數(shù)矩陣條件數(shù)增大，應(yīng)用輔助函數(shù)ψ(x)來分別取代式(12)中的二階導(dǎo)數(shù)項，則有

(13)

取Δx=0.1，配點數(shù)n=11，φ(x)和ψ(x)分別按式(1)、式(2)計算，其中ψ(x)為輔助徑向基函數(shù)。代入各配點坐標，按照式(6)、式(7)類似的步驟，可求得形函數(shù)

N=

撓度的一階和二階變率可分別表示為：

配點法與此相對應(yīng)，不但要滿足微分方程，而且還必須滿足簡2端邊界上撓度、撓度二階變率的條件，即在邊界上也滿足微分方程。

算例中，直接調(diào)用MATLAB中的fsolve函數(shù)(最小二乘法)來求解非線性方程組。計算結(jié)果如圖1，可以看出數(shù)值解與解析解吻合得很好。

圖1 簡支彈性幾何非線性梁的縱橫彎曲Fig. 1 Combined vertical and horizontal bending of simply supported elastic beam with geometrical nonlinearity

表1為筆者提出的方法在配點處的計算值和相對誤差?？梢?，文中方法具有計算過程簡單、精度高的特點，其計算精度與修正小波伽遼金法不相上下[20]。

表1 各配點計算結(jié)果與相對誤差Table 1 Computed results and relative errors of each collocation point

3.2 非線性彈性地基梁的大撓度彎曲

選擇一個置于非線性彈性地基上的軸向不可自由伸縮簡支彈性梁[20]，梁的長度l=1 m，抗彎剛度EI=1 N·m2，抗拉(壓)剛度EA=100 N，地基剛度系數(shù)k=100 N/m4，則作用于梁上的橫向分布力p(x)(單位：N/m)為：

梁的撓曲線微分方程和邊界條件：

做無量綱化處理，梁的撓曲線微分方程及邊界條件如方程(14)，此方程含有泛函表達式，增加了求解的難度：

(14)

式中：H、F均為無量綱參數(shù)，H=Al2/I，F(xiàn)=kl3/(EI)；P(x)=p(x)l3/(EI)。

將方程(14)中的幾何非線性項做近似處理(弱非線性時可行)：

則，此時近似系統(tǒng)的理論解為：

運用徑向基函數(shù)配點法，取Δx=0.05，均勻配點數(shù)n=21。內(nèi)插2個高斯積分點，采用三點高斯-勒讓德積分離散泛函表達式，將泛函表達式化為代數(shù)表達式引入非線性方程組。為便于程序?qū)崿F(xiàn)，高斯積分點既用于數(shù)值積分，又用于徑向基函數(shù)配點計算，總體方案為23個非均勻配點。

為了避免邊界出現(xiàn)數(shù)值振蕩，提出撓度、撓度的二階變率、撓度的四階變率聯(lián)合插值的徑向基函數(shù)表達式：

(15)

為避免式(15)中徑向基函數(shù)的高階導(dǎo)數(shù)使得系數(shù)矩陣條件數(shù)增大，應(yīng)用輔助函數(shù)θ(x)和ψ(x)來分別取代式中的二階和四階導(dǎo)數(shù)項，則有插值函數(shù)(16)：

(16)

式中：θ(x)、ψ(x)、φ(x)分別按式(1)、式(2)、式(3)計算。

由式(16)計算出ω=[ω1，ω2，…，ωn，v1，v2，v3，v4]T(v1、v2、v3、v4為附加未知量，分別為梁的兩個邊界上撓度的二階變率和四階變率)。

配點法與此相對應(yīng)，不但要滿足微分方程，而且還必須滿足兩端邊界上撓度、撓度二階變率、撓度四階變率的條件，即在邊界上也滿足微分方程，附加四階變率邊界條件的約束：ω(4)(0)=ω(4)(1)=0，這樣才能很好地消除數(shù)值振蕩的影響，獲得滿意的結(jié)果。

從圖2可以看出，數(shù)值解與解析解吻合得很好。

圖2 簡支非線性彈性地基梁的橫向大撓度彎曲Fig. 2 Transverse large deflection bending of simply supported nonlinear elastic foundation beam

圖3給出了各配點的相對誤差，最大誤差為0.67%；本算例用ANSYS軟件求解，相對誤差高達10%～12%，應(yīng)用修正小波伽遼金法求解，相對誤差也有3%～6%[20]。筆者提出的方法具有計算過程簡單、精度高的優(yōu)越性。

圖3 各配點計算結(jié)果的相對誤差Fig. 3 Relative error of each collocation point

4 結(jié) 論

筆者將徑向基函數(shù)逼近與配點法相結(jié)合，提出了對應(yīng)的插值表達式，構(gòu)造了一種適用于強非線性梁求解的數(shù)值算法。算例分析表明，與傳統(tǒng)的數(shù)值計算方法相比，筆者提出的方法具有以下優(yōu)越性：

1)無網(wǎng)格，計算過程簡，計算精度高。

2)針對具體的工程問題，配點位置可以根據(jù)需要靈活選取均勻或非均勻配點，且在較少配點的情況下也能夠獲得較高的計算精度，針對非線性問題文中方法具有一定優(yōu)勢。

3)關(guān)于徑向基函數(shù)的選擇，緊支柱正定徑向基函數(shù)相比國外的全局徑向基函數(shù)，前者的計算結(jié)果不依賴于參數(shù)的選擇，可以根據(jù)計算精度和效率的要求設(shè)置支撐域半徑，更適合在工程實踐中推廣。