陳 炎, 張新東

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

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

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

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

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

1 重心插值配點(diǎn)法求解Volterra積分-微分方程

1.1 重心Lagrange插值

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

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

(1)

其中,Lj(x)為L(zhǎng)agrange插值基函數(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ù)進(jìn)行插值, 在2007年，F(xiàn)loater和Hormann提出了一種重心有理插值形式, 該形式無(wú)論是在等距節(jié)點(diǎn)下, 還是在第二類Chebyshev這種特殊分布的節(jié)點(diǎn)下都具有很高的計(jì)算精度[16].

將(2)重新定義為

(6)

式中的d是選定的一個(gè)正整數(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)造

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

(7)

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

(8)

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

(9)

將(9)寫成矩陣的形式

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

(10)

2 數(shù)值算例

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

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

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

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

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

表1 例1的計(jì)算誤差

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

圖1 等距節(jié)點(diǎn)下例1的計(jì)算誤差

圖2 Chebyshev節(jié)點(diǎn)下例1的計(jì)算誤差

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

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

表3 例2的計(jì)算誤差

圖3 等距節(jié)點(diǎn)下例2的計(jì)算誤差

圖4 Chebyshev節(jié)點(diǎn)下例2的計(jì)算誤差

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

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

表4 例3的計(jì)算誤差

圖5 等距節(jié)點(diǎn)下例3的計(jì)算誤差

圖6 Chebyshev節(jié)點(diǎn)下例3的計(jì)算誤差

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

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

表5 例4的計(jì)算誤差

圖7 等距節(jié)點(diǎn)下例4的計(jì)算誤差

圖8 Chebyshev節(jié)點(diǎn)下例4的計(jì)算誤差

3 結(jié)論

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