葉林薇, 李書琰, 晏 云*

（1.閩南師范大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院,福建 漳州 363000；2.福建省粒計算及其應(yīng)用重點實驗室,福建 漳州 363000；3.蓼堤二中,河南 商丘 476900）

有限差分方法以其簡單、精確和通用性強等特點被廣泛運用于微分方程的求解過程中，有限差分方法的核心思想在于利用函數(shù)值的加權(quán)和來近似導(dǎo)數(shù)值，因此確定有限差分公式中的權(quán)重顯得尤其關(guān)鍵[1-7].求取有限差分公式權(quán)重的傳統(tǒng)方法是在某點進行Taylor展開，將Taylor級數(shù)展開式用矩陣的形式表示，再應(yīng)用Gauss 消元法以求取所需有限差分公式的權(quán)重，但該方法針對不同形式的差分公式，展開階數(shù)不同，且局限于等距網(wǎng)格，計算繁雜不夠高效[8].

為此，Bengt Fornberg 提出一系列基于插值思想的算法[9-12].該類型算法的主要好處是：1)插值法只要求節(jié)點互異，并沒有要求節(jié)點間等距，故而可用于構(gòu)造不等距網(wǎng)格的差分公式；2)插值的結(jié)構(gòu)性好，易于編程實現(xiàn)．1988年，他通過Lagrange型插值多項式，利用多項式系數(shù)求導(dǎo)的遞歸關(guān)系，推導(dǎo)出任意間距網(wǎng)格上任意階導(dǎo)數(shù)對應(yīng)的有限差分公式權(quán)重的計算方法，該算法主要是針對顯式有限差分格式[9].1998 年，他將其拓展到隱式有限差分格式[10].2006 年，針對散亂節(jié)點下的有限差分格式中的權(quán)重計算，他給出徑向基函數(shù)插值的方法[11].2020年，他基于帶有一階導(dǎo)數(shù)信息的Lagrange型Hermite插值，給出帶有一階導(dǎo)數(shù)的緊致差分格式的權(quán)重求解方法[12].以上這些方法相比Taylor展開都更加高效，且易編程實現(xiàn).

Bengt Fornberg 計算有限差分格式權(quán)重的主要思想是利用插值多項式的導(dǎo)數(shù)來逼近原函數(shù)的導(dǎo)數(shù).例如若已知n+1 個互異的插值節(jié)點x0,x1,…,xn的函數(shù)值f(xi)及直到mi階的導(dǎo)數(shù)值(不一定需要使用到)信息，則存在插值多項式具有一般形式為

對pm(x)求q(q=0,1,…)階導(dǎo)后，再代入逼近節(jié)點，可得到

顯然，顯式的有限差分格式權(quán)重，即ci，k可直接從上式右端獲取.另一方面，如對求導(dǎo)后的式子分別代入所需點再進行適當(dāng)組合可獲得緊致差分格式的一般形式為

在此啟發(fā)下，給出三種新的算法來求解在任意網(wǎng)格間距下的有限差分格式權(quán)重，它們分別基于Newton型插值多項式、帶有連續(xù)階數(shù)導(dǎo)數(shù)信息的Lagrange型Hermite插值多項式與帶有連續(xù)階數(shù)導(dǎo)數(shù)信息的Newton 型Hermite 插值多項式，因此分別稱為Newton 型算法，Lagrange 型Hermite 算法及Newton 型Hermite算法.特別的，后面兩種算法比文獻[12]中的算法適用的有限差分公式類型更廣.

1 Newton型算法

若已知f(x)在n+1個互異節(jié)點x0,x1,…,xn的函數(shù)值，則Newton型插值多項式[13]可表示為

其中：f[x0] =f(x0) 稱為f(x) 關(guān)于節(jié)點x0的零階差商.一般地，f[x0,x1,…,xk-1,xk] =稱為f(x)關(guān)于節(jié)點x0,x1,…,xk的k階差商.

觀察式(1)可發(fā)現(xiàn)，其每一項都可寫為差商f[x0,x1,…,xk]與乘積相乘的形式.不妨將乘積部分記為M0(x) =1,Mk(x) =(x-x0)(x-x1)…(x-xk-1),k=1,2,…,n, 顯然乘積部分為k次多項式且存在遞推關(guān)系Mk(x) =Mk-1(x)(x-xk-1).由此可將Newton型插值多項式重新表示為

對式(2)求q(q=0,1,…)階導(dǎo)，則有從而為得到所需權(quán)重，只需分別考慮差商部分各點函數(shù)值f(xi)前的系數(shù)以及乘積部分求導(dǎo)后的結(jié)果．

首先考慮差商部分各點函數(shù)值f(xi)前的系數(shù)lk，i, 其中k代表k階差商部分共涉及k+1 個互異的節(jié)點，從而k為整數(shù)且i=0,1,…,k.當(dāng)k=0時，有f[x0] =f(x0), 即l0,0=1.當(dāng)k＞0時，由差商可由函數(shù)值線性組合表出的性質(zhì)[13]，可知

即在k階差商f[x0,x1,…,xk]中，f(xi)前的系數(shù)為

接下來考慮對乘積部分Mk(x)求q階導(dǎo)．若q=0, 即不求導(dǎo)的情況，Mk(0)(x) =Mk(x).若q＞0， 結(jié)合乘積部分遞推關(guān)系式以及Leibniz高階求導(dǎo)公式，則計算結(jié)果可分為以下三種情況，即

由Newton型插值多項式求導(dǎo)后的一般形式，即式(3)，以及上述差商部分和乘積部分的結(jié)論，即式(5)～(6)，可得q階導(dǎo)數(shù)有限差分公式中f(xi)前的系數(shù)為

2 Lagrange型Hermite算法

若已知n+1個互異節(jié)點xi(i=0,1,…,n)上的函數(shù)值和直到mi階的導(dǎo)數(shù)值，即則存在m(這里)次的Lagrange型Hermite插值多項式[13]，則有

根據(jù)Leibniz高階求導(dǎo)公式，對式(9)求q(q=0,1,…)階導(dǎo)有

下面分別對求和部分sk，i(x)與求積部分πi(x)的q階導(dǎo)數(shù)情況進行考慮．

首先來看求和部分的導(dǎo)數(shù)(x).若q=0, 即在不求導(dǎo)的情況下，(x) =sk，i(x).若q＞0, 因為sk，i(x)是最低次數(shù)為k次且最高次數(shù)為mi次的多項式，所以求導(dǎo)結(jié)果可分為以下三種情況，即

接下來考慮求積部分的導(dǎo)數(shù)πi(q)(x).記hi(x) =(x-xi)mi+1, 根據(jù)多個函數(shù)相乘的高階求導(dǎo)公式[14]有

由Lagrange 型Hermite 插值多項式求導(dǎo)后的一般形式，即式(10)，以及上述求和部分與求積部分的結(jié)論，即式(11)～(12)，可得q階導(dǎo)數(shù)有限差分公式中f(k)(xi)前的系數(shù)為

3 Newton型Hermite算法

若已知n+1 個互異節(jié)點xi(i=0,1,…,n) 上的函數(shù)值和直到mi階的導(dǎo)數(shù)值，即f(xi),f′(xi),…,f(mi)(xi),可得到次的Newton型Hermite插值多項式[13].則有

對式(15)求q(q=0,1,…)階導(dǎo)，有

觀察式(16)結(jié)構(gòu)發(fā)現(xiàn)獲得有限差分格式權(quán)重的關(guān)鍵在于分析每項中差商部分f(k)(xi)(k=1,2,…,mi)前的系數(shù)與乘積部分求導(dǎo)后的結(jié)果，下面分別對這兩部分進行分析．

先來考慮差商部分f(k)(xi)前的系數(shù)．將差商中各節(jié)點重數(shù)之和記為μ, 則

當(dāng)1 ≤μ≤m0+1時，即僅考慮節(jié)點x0的差商，此時μ=k+1.由k階導(dǎo)數(shù)與k階差商之間的關(guān)系，可得

顯然，在f([x0,k+1])中，f(k)(x0)前的系數(shù)為

此時存在互異節(jié)點x0,x1,…,xj(j=1,2,…,n)的差商公式[15]可表示為

接下來考慮乘積部分求導(dǎo)后的結(jié)果.不妨令乘積部分為(x),μ(0 ≤μ≤m)表示某一乘積部分中多項式的次數(shù)，則乘積部分的通項可表示為

觀察此可發(fā)現(xiàn)，乘積部分存在如下的遞推關(guān)系

接下來對(x)求q(q=0,1,…)階導(dǎo)數(shù).若q=0,即不求導(dǎo)的情況，若q＞0, 結(jié)合乘積部分遞推關(guān)系式(22)和Leibniz高階求導(dǎo)公式，求導(dǎo)結(jié)果共分為以下兩種情況，即有

由Newton 型Hermite 插值多項式求導(dǎo)后的一般形式，即式(16)，以及上述差商部分和乘積部分的結(jié)論，即式(18)、式(20)和式(23)，可得q階導(dǎo)數(shù)有限差分公式中f(k)(xi)前的系數(shù)為

4 數(shù)值實驗

下面選取三個具有代表性的算例來說明三種算法如何確定顯式有限差分格式的權(quán)重、隱式緊致有限差分格式的權(quán)重以及帶有連續(xù)階導(dǎo)數(shù)信息的有限差分格式的權(quán)重．

首先，對于顯式有限差分格式的權(quán)重的求解，以二階導(dǎo)數(shù)中心差分格式為例，不妨選取節(jié)點xi=-1,0,1, 則二階導(dǎo)數(shù)中心差分格式表現(xiàn)形式為

于是可在節(jié)點處利用Newton型算法以獲取權(quán)重c(02),c(12),c(22)來近似導(dǎo)數(shù)f0′′, 即有

由式(7)知，所需信息有插值節(jié)點、逼近節(jié)點以及求導(dǎo)階數(shù).針對二階導(dǎo)數(shù)中心差分格式，插值節(jié)點為xi=-1,0,1, 那么逼近節(jié)點為=0, 求導(dǎo)階數(shù)q=2, 將上述相關(guān)數(shù)據(jù)應(yīng)用到Newton型算法中，可計算出權(quán)重為

顯然代入權(quán)重后的式(26)即為二階導(dǎo)數(shù)的中心差分格式.利用文獻[9]中算法可得到相同的權(quán)重結(jié)果．

特別地，若使用Lagrange 型Hermite 算法或Newton 型Hermite 算法求解二階導(dǎo)數(shù)中心差分格式的權(quán)重，根據(jù)式(13)和式(24)可知所需信息與Newton型算法相比，除插值節(jié)點、逼近節(jié)點以及求導(dǎo)階數(shù)以外，還需已知各插值節(jié)點所對應(yīng)最高階導(dǎo)數(shù)的階數(shù).通過式(26)知在計算二階導(dǎo)數(shù)中心差分格式的權(quán)重中，只涉及到函數(shù)值，沒有涉及到導(dǎo)數(shù)值，因此各插值節(jié)點所對應(yīng)最高階導(dǎo)數(shù)的階數(shù)均為零，即mi=0,0,0.這表明，Lagrange型Hermite算法和Newton型Hermite算法是Newton型算法的推廣.

其次，考慮隱式緊致有限差分格式權(quán)重的求解，以一階導(dǎo)數(shù)的三點四階緊致差分格式為例．這里同樣選取節(jié)點xi=-1,0,1, 則一階導(dǎo)數(shù)的三點四階緊致差分格式表現(xiàn)形式為

為得到此格式，可先考慮如下式子

為能夠使用Lagrange型Hermite算法或Newton型Hermite算法來確定式(29)的權(quán)重，將上式整理成

這樣問題就轉(zhuǎn)化為用點xi=-1,0,1 處的相應(yīng)函數(shù)值及一階導(dǎo)數(shù)值的加權(quán)和來近似五階導(dǎo)數(shù)f0(5).此時，插值節(jié)點xi=-1,0,1, 逼近節(jié)點=0, 求導(dǎo)階數(shù)q=5,各插值節(jié)點所對應(yīng)最高階導(dǎo)數(shù)的階數(shù)mi=1,1,1, 將上述相關(guān)數(shù)據(jù)代入Lagrange型Hermite算法或Newton型Hermite算法中可計算出

利用文獻[7]中的算法也可得出一致的權(quán)重結(jié)果.接著將權(quán)重代回式(30)后，等式兩端同乘h4/180并舍去截斷誤差h4f(5)0, 即可得一階導(dǎo)數(shù)的三點四階緊致差分格式.

最后，考慮文獻[16]中的公式

式(32)涉及到函數(shù)值、一階導(dǎo)數(shù)值、二階導(dǎo)數(shù)值和三階導(dǎo)數(shù)值，是沒有辦法用文獻[12]中的方法得到的.而根據(jù)Lagrange型Hermite算法和Newton型Hermite算法，可以得到

5 結(jié)論

Newton型算法、Lagrange型Hermite算法及Newton型Hermite算法均可通過計算機編程實現(xiàn)，免去大量的手動計算， 從而可快速得到所需權(quán)重．其中Newton 型算法與文獻[9]中的方法相比，每新增一個節(jié)點，Newton 型算法計算工作不必重頭開始，計算量少些(但鑒于如今的計算機速度，在得到權(quán)重的時間上已經(jīng)看不到兩者區(qū)別)． Lagrange 型Hermite 算法和Newton 型Hermite 算法與文獻[12]中的方法相比，它們除可以用于只含有一階導(dǎo)數(shù)的情況，還適用于出現(xiàn)多個不同階導(dǎo)數(shù)的情況．