雍龍泉

(陜西理工大學(xué) 數(shù)學(xué)與計(jì)算機(jī)科學(xué)學(xué)院, 陜西 漢中 723001)

非線性兩點(diǎn)邊值問題如下：

(1)

其在物理學(xué)、流體力學(xué)、材料力學(xué)、波動(dòng)力學(xué)、控制理論等領(lǐng)域有著廣泛的應(yīng)用和重要的理論價(jià)值[1-3]。其研究主要分為兩個(gè)方面：一是研究解的存在性；二是計(jì)算其數(shù)值解。關(guān)于邊值問題解的存在性和唯一性的研究，在20世紀(jì)出現(xiàn)了大量文獻(xiàn)，至今仍不斷發(fā)表新的研究成果。但在解的存在性方面，目前研究較多的是正解的存在性[4-5]。

鑒于可以解析求解的問題甚少，所以邊值問題的求解是較為困難的。常見計(jì)算二階兩點(diǎn)邊值問題數(shù)值解的方法有打靶法、有限差分法等[1,6-8]。打靶法主要思路是恰當(dāng)選擇和調(diào)整初值的條件，對(duì)一系列初值問題進(jìn)行求解，使其接近給定的邊界條件。當(dāng)g(t,x,x′)=p(t)x′+q(t)x+r(t)時(shí)，問題(1)則為線性兩點(diǎn)邊值問題。針對(duì)線性兩點(diǎn)邊值問題，可以采用有限差分法將其轉(zhuǎn)化為線性方程組進(jìn)行求解。對(duì)于非線性兩點(diǎn)邊值問題，其離散化后所得到的方程組是非線性的，因此就需要研究求解非線性方程組的快速算法。

本文研究了無奇異點(diǎn)的非線性兩點(diǎn)邊值問題的數(shù)值解，通過將其離散化，得到非線性方程組：

(2)

其中:xi表示函數(shù)x(t)在離散點(diǎn)ti處的函數(shù)值，即xi=x(ti)；fi(x1,x2,…,xn):D?Rn→R是非線性函數(shù)，且充分可導(dǎo)；i=1,2,…,n(詳細(xì)轉(zhuǎn)化過程在本文第2部分)。

記向量x=(x1,x2,…,xn)，向量函數(shù)F(x)=(f1(x),f2(x),…,fn(x))T，則方程組(2)等價(jià)于如下非線性方程：

F(x)=0

(3)

下面采用高階牛頓法求解非線性方程組。

1 高階牛頓法

求解非線性方程常見的方法有牛頓迭代法、修正牛頓迭代法以及智能優(yōu)化算法等[9-15]。近年來，各類高階牛頓迭代法成為研究的熱點(diǎn)，相繼出現(xiàn)了五階牛頓迭代法[16-21]、七階牛頓迭代法[22-24]、八階牛頓迭代法[25-26]、九階牛頓迭代法[27-28]等。當(dāng)然，收斂階數(shù)越高，需要付出的計(jì)算代價(jià)也就越大。因此，在構(gòu)造高階牛頓迭代法求解非線性方程組時(shí)，既需要考慮收斂階，更需要考慮計(jì)算復(fù)雜性。文獻(xiàn)[21]中給出的五階牛頓迭代法公式簡單，且計(jì)算量少，因此已成功應(yīng)用于求解線性規(guī)劃等問題[29-30]。

五階牛頓迭代法步驟如下：

1) 輸入

非線性方程組F(x)=(f1(x),f2(x),…,fn(x))T=0；

非線性方程組對(duì)應(yīng)的雅克比矩陣F′(x)；

初始迭代點(diǎn)x(0)∈Rn，終止條件ε；

令k=0。

2) 開始計(jì)算

③x(k+1)=z(k)+([F′(x(k))]-1-2[F′(y(k))]-1)F(z(k))；

④ 如果(||x(k+1)-x(k)||+||F(x(k))||)<ε，停止，輸出計(jì)算結(jié)果，否則進(jìn)行下一步；

⑤k=k+1，轉(zhuǎn)步驟①。

3) 輸出計(jì)算結(jié)果

x*=x(k+1),F(x*)=F(x(k+1))

算法說明：實(shí)際計(jì)算過程中，若雅克比矩陣奇異或接近奇異，可以采用阻尼牛頓法進(jìn)行處理。文獻(xiàn)[21]證明了該方法在求解非線性方程組時(shí)具有五階收斂速度。以下應(yīng)用該五階牛頓法求解非線性兩點(diǎn)邊值問題。

2 應(yīng)用于非線性兩點(diǎn)邊值問題

下面給出兩個(gè)非線性兩點(diǎn)邊值問題，通過將其轉(zhuǎn)化為非線性方程組后，采用上述五階牛頓迭代法進(jìn)行求解，程序采用Matlab R2009a編寫。設(shè)置x(0)=(1,1,…,1)∈Rn，終止條件ε=1×10-15。

算例1考慮如下非線性兩點(diǎn)邊值問題：

對(duì)區(qū)間t∈[0,1]做離散化處理，插入n個(gè)點(diǎn)：

這些分點(diǎn)滿足

0=t0

于是區(qū)間被分成n+1等分，記

x(t0)=x0=0,x(t1)=x1,x(t2)=

x2,…,x(tn)=xn,x(tn+1)=xn+1=1

下面近似計(jì)算函數(shù)x(t)在這些插入點(diǎn)處的函數(shù)值xi,i=1,2,…,n。

利用有限差分法可知

得到如下有n個(gè)變量(xi,i=1,2,…,n)的非線性方程組：

fi(x1,x2,…,xn)=xi-1-2xi+xi+1+h2xi=0,

i=1,2,…,n

即：

若插入n=15個(gè)點(diǎn)，采用上述五階牛頓迭代法進(jìn)行計(jì)算，表1給出了這15個(gè)插入點(diǎn)處的數(shù)值解和精確解。

表1 算例1對(duì)應(yīng)的數(shù)值解與精確解(n=15)

圖1、2分別給出了n=15和n=30對(duì)應(yīng)的數(shù)值解和精確解的示意圖。

圖1 n=15時(shí)對(duì)應(yīng)的數(shù)值解和精確解的示意圖

圖2 n=30時(shí)對(duì)應(yīng)的數(shù)值解和精確解的示意圖

算例2考慮如下非線性兩點(diǎn)邊值：

該問題的解析解不易求出，下面直接計(jì)算數(shù)值解。

對(duì)區(qū)間t∈[0,1]做離散化處理，插入n個(gè)點(diǎn)：

這些分點(diǎn)滿足

0=t0

于是區(qū)間被分成n+1等分，記

x(t0)=x0=0,x(t1)=x1,x(t2)=

x2,…,x(tn)=xn,x(tn+1)=xn+1=1

下面近似計(jì)算函數(shù)x(t)在這些插入點(diǎn)處的函數(shù)值xi,i=1,2,…,n。

利用

i=1,2,…,n

得到如下n個(gè)變量(xi,i=1,2,…,n)的非線性方程組：

即：

若插入n=15個(gè)點(diǎn)，采用上述五階牛頓迭代法進(jìn)行計(jì)算，表2給出了這15個(gè)插入點(diǎn)處的數(shù)值解和和對(duì)應(yīng)的fi(x)值,i=1,2,…,15。

圖3、4分別給出了n=15和n=30的數(shù)值解示意圖。

表2 算例2對(duì)應(yīng)的數(shù)值解(n=15)

圖3 n=15時(shí)對(duì)應(yīng)的數(shù)值解的示意圖

3 結(jié)束語

本文通過把非線性兩點(diǎn)邊值問題的數(shù)值解轉(zhuǎn)化為非線性方程組，進(jìn)而采用高階牛頓迭代法進(jìn)行求解，且對(duì)文獻(xiàn)[21]中的數(shù)值結(jié)果進(jìn)行了更正。數(shù)值計(jì)算結(jié)果表明：該方法計(jì)算速度快(即使插入200個(gè)分點(diǎn)，計(jì)算耗時(shí)也不到2 s)，且計(jì)算精度高，對(duì)此類問題較為有效。進(jìn)一步可以考慮把該方法應(yīng)用到實(shí)際工程問題中(如需源代碼，請(qǐng)與筆者聯(lián)系)。