Volterra積分-微分方程的高精度數(shù)值算法研究

陳 炎, 張新東

(新疆師范大學(xué) 數(shù)學(xué)科學(xué)學(xué)院,新疆 烏魯木齊 830017)

20世紀初,生物學(xué)家Ancona和數(shù)學(xué)家Volterra提出了經(jīng)典的“地中海-鯊魚模型”.1931年,Volterra在“地中海-鯊魚模型”基礎(chǔ)上提出了經(jīng)典的單種群增長積分-微分模型[1]

Volterra在研究種群動力學(xué)增長模型時將遺傳性影響因素引入, 得到著名的Volterra積分-微分方程[2]

其中,K(x,t)為積分核函數(shù),f(x)為自由項.Volterra型積分-微分方程在工程和應(yīng)用科學(xué)問題的研究中都有廣泛的應(yīng)用,如納米流體力學(xué)、擴散過程、中子擴散化學(xué)和電化學(xué)過程、熱質(zhì)傳遞、彈性力學(xué)與動態(tài)接觸、黏彈性和反應(yīng)器動力學(xué)、流行病建模等[3].正是由于Volterra型積分-微分方程在各個領(lǐng)域的廣泛應(yīng)用,所以關(guān)于Volterra積分-微分方程的數(shù)值研究也數(shù)不勝數(shù),比如譜方法[4-5]、有限差分法[6]、隱式Runge-Kutta方法[7-8]、歐拉法[9]、線性多步法[10]、多項式樣條配置法[11]和有限元方法[12]等,以上方法具有各自特點的同時也存在一些不足,例如有限元方法要求計算時劃分稠密的計算網(wǎng)格,有限差分法要求很小的差分步長,這將極大地降低分析問題的效率,而本文采用的重心插值配點法不用劃分計算網(wǎng)格,并且計算程序編寫非常簡便,是一種高精度、高效率的數(shù)值計算方法.

重心插值配點法已經(jīng)獲得了廣泛應(yīng)用,王兆清等[13]運用重心插值配點法分析了圓環(huán)變形及屈曲問題; 鄧楊芳等[14]運用重心插值配點法對Cahn-Hilliard方程進行了離散求解.過去一段時間,關(guān)于Volterra積分-微分方程的重心插值配點法研究層出不窮,但是目前很多數(shù)值研究都沒有改變積分項中的未知函數(shù)的形式,本文將未知函數(shù)變?yōu)橐浑A導(dǎo)數(shù)、二階導(dǎo)數(shù)和三階導(dǎo)數(shù),采用重心插值配點法進行數(shù)值模擬,并觀察其誤差變化.本文首先介紹了重心插值配點法,然后構(gòu)造了n階Volterra積分-微分方程的離散格式,最后使用Matlab對四個不同階數(shù)的數(shù)值算例進行了數(shù)值實現(xiàn),并驗證了方法的有效性.

1 重心插值配點法求解Volterra積分-微分方程

1.1 重心Lagrange插值

重心Lagrange插值是Lagrange插值的改進,2004年,Berrut和Trefethen提出了重心Lagrange插值,重心Lagrange 插值在一些特殊節(jié)點下具有很好的數(shù)值穩(wěn)定性[15].

設(shè)有n+1個不同的插值節(jié)點xj(j=0,1,…,n)和相對應(yīng)的一組實數(shù)fj, 則插值多項式就可以寫成Lagrange插值公式

(1)

其中,Lj(x)為Lagrange插值基函數(shù),

令l(x)=(x-x0)(x-x1)…(x-xn), 定義重心權(quán)

(2)

也就是wj=1/l′(xj), 則

(3)

將式(3)代入(1), 得到另一種Lagrange插值形式

(4)

利用插值常數(shù)1，可得下面恒等式

(5)

將式(5)兩邊分別除去式(4)的兩邊, 得到重心Lagrange插值公式

1.2 重心有理插值

重心有理插值就是利用有理函數(shù)進行插值, 在2007年，F(xiàn)loater和Hormann提出了一種重心有理插值形式, 該形式無論是在等距節(jié)點下, 還是在第二類Chebyshev這種特殊分布的節(jié)點下都具有很高的計算精度[16].

將(2)重新定義為

(6)

式中的d是選定的一個正整數(shù), 指標(biāo)集Jk={i∈O:k-d≤i≤k},O={0,1,…,n}, 則重心有理插值公式就可以寫為

重心Lagrange插值公式與重心有理插值具有相同的形式,但插值權(quán)重的選擇不同，重心Lagrange插值和重心有理插值的插值權(quán)分別由(2)和(6)確定.

1.3 Volterra積分-微分方程離散格式構(gòu)造

針對如下形式的α階Volterra積分-微分方程, 我們用重心插值配點法對其進行離散

(7)

其中,y(α)和y(β)分別為y的α階導(dǎo)數(shù)和β階導(dǎo)數(shù),α和β均為自然數(shù).將區(qū)間[a,b]離散為n+1個點x0,x1,…,xn,令y0,y1,…,yn為未知函數(shù)y(x)在節(jié)點x0,x1,…,xn處的函數(shù)值. 利用重心插值近似未知函數(shù)

(8)

將(8)代入(7), 并將離散節(jié)點xi(i=0,1,…,n)代入, 得到

(9)

將(9)寫成矩陣的形式

(D(α)-Ψ-ΦK)y=f，

(10)

2 數(shù)值算例

在數(shù)值算例模擬中, 表格中的絕對誤差和相對誤差與圖中的絕對誤差分別定義為

err=‖yn-ya‖2,rerr=‖yn-ya‖2/‖ya‖2,err′=yn-ya,

其中, ‖·‖為向量的歐式范數(shù);yn和ya分別為未知函數(shù)的數(shù)值解和解析解. 為了方便比較, 以下四個算例的解析解和核函數(shù)K(x,t)均保持不變, 改變積分項中的未知函數(shù)的形式和方程的階數(shù). 以下算例中的I為單位矩陣, 實現(xiàn)以下算例所使用的軟件是MATLAB R2021a, 電腦型號是聯(lián)想拯救者R7000P.

例1考慮一般形式的一階Volterra積分-微分方程

在21個等距節(jié)點, 21個Chebyshev節(jié)點, 8個高斯積分點, 有理插值參數(shù)d=10的條件下, 使用重心插值配點法的誤差見表1, 從表中數(shù)據(jù)我們可以看出,在等距節(jié)點下, 重心有理插值配點法的計算精度遠高于重心Lagrange 插值配點法，在Chebyshev節(jié)點下, 重心有理插值配點法的計算精度略高于重心Lagrange插值配點法; 但總體來說Chebyshev節(jié)點的計算精度要高于等距節(jié)點. 從圖1我們可以看到, 在等距節(jié)點下, 重心有理插值的計算誤差要比重心Lagrange插值小很多, 從圖2我們可以發(fā)現(xiàn), 在第二類Chebyshev節(jié)點下, 兩種配點法的誤差變化趨勢有著很大的不同, 重心Lagrange插值配點法的計算誤差較為穩(wěn)定, 而重心有理插值配點法的計算誤差波動較大. 本算例使用的8個高斯積分點的節(jié)點和權(quán)重如表2.

表1 例1的計算誤差

表2 8個高斯積分點和高斯積分權(quán)重

圖1 等距節(jié)點下例1的計算誤差

圖2 Chebyshev節(jié)點下例1的計算誤差

例2考慮如下形式的二階Volterra積分-微分方程

在21個等距節(jié)點,21個Chebyshev節(jié)點,8個高斯積分點,有理插值參數(shù)d=10的條件下,使用重心插值配點法的誤差見表3,從表中數(shù)據(jù)我們可以看出,在等距節(jié)點下， 重心有理插值配點法的計算精度遠高于重心Lagrange插值配點法，在Chebyshev節(jié)點下,重心有理插值配點法的計算精度略高于重心Lagrange插值配點法; 但總體來說Chebyshev節(jié)點的計算精度要高于等距節(jié)點.從圖3我們可以看到,重心有理插值的計算誤差要比重心Lagrange插值小很多;從圖4我們可以發(fā)現(xiàn),在第二類Chebyshev節(jié)點下,重心有理插值配點法的計算誤差較為穩(wěn)定,而重心Lagrange插值配點法的計算誤差變化趨勢較大.本算例使用的8個高斯積分點的節(jié)點和權(quán)重同表2.

表3 例2的計算誤差

圖3 等距節(jié)點下例2的計算誤差

圖4 Chebyshev節(jié)點下例2的計算誤差

例3考慮如下形式的三階Volterra積分-微分方程

在21個等距節(jié)點,21個Chebyshev節(jié)點,8個高斯積分點,有理插值參數(shù)d=10的條件下,使用重心插值配點法的誤差見表4,從表中數(shù)據(jù)我們可以看出,在等距節(jié)點下,重心有理插值配點法的計算精度遠高于重心Lagrange插值配點法,在Chebyshev節(jié)點下,重心有理插值配點法的計算精度略高于重心Lagrange插值配點法; 但總體來說Chebyshev節(jié)點的計算精度要高于等距節(jié)點.從圖5我們可以看到,在等距節(jié)點下,重心有理插值的計算誤差要比重心Lagrange插值小很多,從圖6我們可以發(fā)現(xiàn),在第二類Chebyshev節(jié)點下,重心有理插值配點法的誤差較為穩(wěn)定.本算例使用的8個高斯積分點的節(jié)點和權(quán)重同表2.

表4 例3的計算誤差

圖5 等距節(jié)點下例3的計算誤差

圖6 Chebyshev節(jié)點下例3的計算誤差

例4考慮如下形式的四階Volterra積分-微分方程

在21個等距節(jié)點, 21個Chebyshev節(jié)點, 8個高斯積分點, 有理插值參數(shù)d=10的條件下, 使用重心插值配點法的誤差見表5, 從表中數(shù)據(jù)我們可以看出, 在等距節(jié)點下, 重心有理插值配點法的計算精度遠高于重心Lagrange插值配點法, 在Chebyshev節(jié)點下, 重心有理插值配點法的計算精度略高于重心Lagrange插值配點法，但總體來說Chebyshev節(jié)點的計算精度要高于等距節(jié)點. 從圖7我們可以看到, 在等距節(jié)點下, 重心有理插值的計算誤差要比重心Lagrange插值小很多, 從圖8我們可以發(fā)現(xiàn), 在第二類Chebyshev節(jié)點下, 重心有理插值配點法的計算誤差較為穩(wěn)定, 而重心Lagrange插值配點法的計算誤差波動較大. 本算例使用8個高斯積分點的節(jié)點和權(quán)重同表2.

表5 例4的計算誤差

圖7 等距節(jié)點下例4的計算誤差

圖8 Chebyshev節(jié)點下例4的計算誤差

3 結(jié)論

本文運用重心插值配點法對Volterra積分-微分方程進行離散求解,并給出四個不同階數(shù)的數(shù)值算例的實驗結(jié)果.通過以上四個算例的誤差分析可知,當(dāng)Volterra積分-微分方程的積分項中含有未知函數(shù)的導(dǎo)數(shù)時,重心插值配點法依然能保持計算穩(wěn)定性和較高的計算精度.通過比較以上四個不同階數(shù)算例的誤差結(jié)果可知,隨著方程階數(shù)和積分項中未知函數(shù)的導(dǎo)數(shù)階數(shù)越來越大,誤差也越來越大,其原因是階數(shù)越高,計算程序中的循環(huán)次數(shù)就越多,計算機完成數(shù)值模擬所需要進行的運算就越多,從而導(dǎo)致誤差積累變大.