任全偉，莊清渠

(1．上海師范大學(xué)數(shù)理學(xué)院，上海 200234;2．華僑大學(xué)數(shù)學(xué)科學(xué)學(xué)院，福建泉州 362021)

0 引言

很多數(shù)學(xué)物理問題的研究最終都可歸結(jié)為微積分方程的研究．人們也一直在用微積分方程模型來解釋和預(yù)知各種自然現(xiàn)象．然而只有少數(shù)比較簡單的微積分方程能夠用初等解法進(jìn)行求解．大部分由實際問題建模得到的微積分方程往往具有比較復(fù)雜的形式，要找出解的解析表達(dá)式是比較困難甚至是不可能的．有些解雖然能解析表達(dá)，但表達(dá)式往往比較復(fù)雜而不切實用．況且實際問題往往也只要在一定的精度范圍內(nèi)求出微積分方程的解在某些點上的值，并不一定要求出它的解析表達(dá)式．因此研究微積分方程的數(shù)值解法對生產(chǎn)生活有著重要的現(xiàn)實意義，為微積分方程尋求快速有效的解法也成為當(dāng)今計算數(shù)學(xué)工作者的主要工作．

在眾多微積分方程的數(shù)值解法中，有限元方法是一種進(jìn)行工程計算的有效方法，自20世紀(jì)50年代起，在航空航天、水利、土木建筑、機(jī)械和流體力學(xué)等方面得到了廣泛的應(yīng)用［1－2］．在研究長為L的懸索橋中，當(dāng)豎向位移為y(x)，負(fù)重為p(x)時，可得如下方程［3］

其中:C1和C2都是正常數(shù);p(x)(≤0或≥0)是在區(qū)間［0，L］上的連續(xù)函數(shù)．為了問題的方便描述，假設(shè)在區(qū)間［0，L］上p(x)≤0．

對該四階微積分方程，文獻(xiàn)［4］給出了混合有限元方法的誤差分析;由于對式(1)進(jìn)行求解的主要難處在于非線性項的處理，文獻(xiàn)［5－6］介紹了一種簡單迭代格式來處理非線性項，并給出了有限元逼近的數(shù)值結(jié)果，但未給出誤差分析．而且在計算過程中，文獻(xiàn)［5－6］的迭代算法收斂速度不是很快，當(dāng)區(qū)間分割較小時，其計算過程將會耗用大量的時間．本文將引入Newton型迭代法［7－9］對非線性項進(jìn)行處理，可以較大提高迭代的收斂速度．另外，給出了有限元逼近的誤差分析，最后通過數(shù)值算例來說明算法的可行性和有效性．

1 有限元逼近

記I=(0，L)，對任意1≤p≤∞，定義Lp(I)={v;vLp＜∞}，其中

令φ=－y″，則式(1)可以寫成

式(2)的變分形式可以寫為:找(φ，y)∈H10(I)×H10(I)，使得

為了給出式(1)的有限元逼近形式，對區(qū)間I進(jìn)行均勻剖分，令h=L/n為區(qū)間分割的長度，即

設(shè)Pm為次數(shù)不超過m的多項式空間．構(gòu)造滿足方程式(1)的邊界條件的分片線性多項式空間X0h如下:

因X0h?H10(I)，即X0h是H10(I)的子空間．則(3)式的有限元逼近形式為:找(φh，yh)∈X0h×X0h，使得

下面證明(2)式弱解的唯一性．

定理1 若(φ1，y1)，(φ2，y2)都是(4)式的解，則它們是相等的．

證明 由假設(shè)與(4)式可得

2 迭代求解

在(4)式中由于積分項的存在，直接計算將會給求解帶來一些困難．因此，引入一種新的迭代格式來處理積分項．首先，記

則(4)式可等價地寫為:找(φh，yh)∈Xh0×Xh0，使得

原問題的求解可歸結(jié)為求解(9)式．在求解(9)式的過程中，可能會遇到下面兩種情況:

1)式(9)的解析解容易求得，則由式(8)可求出式(1)的解．

2)式(9)的解析解不容易求或者求不出，則需采用牛頓型迭代算法對其數(shù)值求解，進(jìn)而求解式(1)．

記f(ξ)=g(ξ)－ξ．則求解式(9)可歸結(jié)為求解f(ξ)=0，其Newton迭代格式為

因此，式(8)的迭代計算過程是:找(φkh，ykh)∈Xh0×Xh0，使得

由上式可知，在迭代過程中需計算f(ξk)，f'(ξk)的值．由于f(ξk)的顯式表達(dá)式通常難以求得，可采用數(shù)值積分來計算f(ξk)，即

對于f'(ξk)采用差商的形式進(jìn)行替代，兩種迭代格式分別稱之為Newton法1，Newton法2．

其中:h為區(qū)間分割的長度，通過求解式(11)可求得式(1)的數(shù)值解．

3 誤差估計

首先，令P01，hy，φ∈X0h分別表示y，φ∈H10(I)到空間X0h上的投影算子，且定義如下:

由文獻(xiàn)［11］可知如下結(jié)論成立

其中C是正常數(shù)．此外，由式(13)、(14)可得類似文獻(xiàn)［12］中引理3．1的證明過程，易得如下結(jié)論．

引理2 若(φ，y)，(φh，yh)分別是式(3)、(4)的解，P01，h，分別是式(15)、(16)中定義的投影算子，則可得

證明:對于式(20)，在式(19)中令η=P01，hy－yh，由式(15)和引理1可得

在上式中令φ =P~01，hφ － φh，由引理1可得

證明完畢．

這里，O(h)表示h的高階無窮小．在證明過程中，不同式子中M可能代表不同的常數(shù)．證明 對于式(22)，由引理1，引理2和三角不等式可得:

即

式可知:

由式(17)可得 (y－yh)'=O(h)．

對于式(23)，由引理1及引理2可知

由式(17)和(22)，可得 φ－φh=O(h)，證明完畢．

4 數(shù)值結(jié)果

例1 p(x)=－π4sin(πx)－π2sin(πx)－2/π．此時方程式(1)的精確解為 y(x)=－sin(πx)，φ =－π2sin(πx)．

首先驗證算法的收斂性．表1給出的是TOL=10－14時誤差隨h的變化關(guān)系．由表易見，數(shù)值解關(guān)于h平方收斂．

再取h=0．01與TOL=10－8和TOL=10－14進(jìn)行計算，相應(yīng)的數(shù)值結(jié)果在表2～3中給出．由表可見，在節(jié)點個數(shù)和迭代精度相同的情況下，Newton型迭代法的收斂速度比Semper迭代法快得多．但是Newton型迭代法需要求解方程(8)兩次(ξk和ξk+h或ξk+h2)．

表1 TOL=10－14，Newton法1所得(例1)的數(shù)值結(jié)果Tab．1 Numerical results using Newton method 1(Example 1)with TOL=10 －14

表2 TOL=10－8，h=0．01 的數(shù)值結(jié)果Tab．2 Numerical results with TOL=10 －8 and h=0．01

例2 p(x)=11/4 －(6e+6e－1－12)/(e1－e－1)－6x，此時(1)的精確解y(x)=5x－(6ex－6e－x)/(e － e－1)+x3，φ =(6ex－ 6e－x)/(e － e－1)－ 6x．

表4～5分別給出了當(dāng)TOL=10－6，h=0．01、0．001時的數(shù)值結(jié)果．由表可見，在節(jié)點個數(shù)和迭代精度相同的情況下，Newton型迭代法的收斂速度比Semper迭代法快得多．

表3 TOL=10－14，h=0．01 的數(shù)值結(jié)果Tab．3 Numerical results with TOL=10 －14 and h=0．01

表4 TOL=10－6，h=0．01 的數(shù)值結(jié)果Tab．4 Numerical results with TOL=10 －6 and h=0．01

表5 TOL=10－6，h=0．001 的數(shù)值結(jié)果Tab．5 Numerical results with TOL=10 －6 and h=0．001

在實際應(yīng)用中，有時難以求得式(1)的精確解．下面從式(1)精確解未知的情況下來說明上述迭代算法的適用性．

例3 p(x)=－x2，此時式(7)的精確解難以求得．

圖1，圖2分別給出了當(dāng)h=0．1、0．01，TOL=10－8時由Newton法1所得的數(shù)值解φh和yh

圖1 h=0．1時，(φ，y)的數(shù)值解:φh，yhFig．1 Numerical solution of(φ，y)with h=0．1:φh，yh

圖2 h=0．01時，(φ，y)的數(shù)值解:φh，yhFig．2 Numerical solution of(φ，y)with h=0．01:φh，yh

5 結(jié)論

用逐次線性有限元和牛頓迭代法相結(jié)合數(shù)值求解四階微積分吊橋模型．給出了詳細(xì)的迭代計算過程以及誤差估計．由數(shù)值結(jié)果可知:

1)采用Newton型迭代法所得的數(shù)值結(jié)果和理論分析是一致的;

2)在相同的條件下，Newton型迭代法收斂速度比Semper所采用的迭代法收斂速度快的多．需要指出的是，由于Newton型迭代法需計算式(13)或(14)值，因此在計算過程中需要求解兩次方程組(ξk和ξk+h或ξk+h2);

3)迭代法可用于求解現(xiàn)實生活中的一些具體問題．