求解有限差分公式权重的三种算法

叶林薇, 李书琰, 晏 云*

（1.闽南师范大学数学与统计学院,福建 漳州 363000；2.福建省粒计算及其应用重点实验室,福建 漳州 363000；3.蓼堤二中,河南 商丘 476900）

有限差分方法以其简单、精确和通用性强等特点被广泛运用于微分方程的求解过程中，有限差分方法的核心思想在于利用函数值的加权和来近似导数值，因此确定有限差分公式中的权重显得尤其关键[1-7].求取有限差分公式权重的传统方法是在某点进行Taylor展开，将Taylor级数展开式用矩阵的形式表示，再应用Gauss 消元法以求取所需有限差分公式的权重，但该方法针对不同形式的差分公式，展开阶数不同，且局限于等距网格，计算繁杂不够高效[8].

为此，Bengt Fornberg 提出一系列基于插值思想的算法[9-12].该类型算法的主要好处是：1)插值法只要求节点互异，并没有要求节点间等距，故而可用于构造不等距网格的差分公式；2)插值的结构性好，易于编程实现．1988年，他通过Lagrange型插值多项式，利用多项式系数求导的递归关系，推导出任意间距网格上任意阶导数对应的有限差分公式权重的计算方法，该算法主要是针对显式有限差分格式[9].1998 年，他将其拓展到隐式有限差分格式[10].2006 年，针对散乱节点下的有限差分格式中的权重计算，他给出径向基函数插值的方法[11].2020年，他基于带有一阶导数信息的Lagrange型Hermite插值，给出带有一阶导数的紧致差分格式的权重求解方法[12].以上这些方法相比Taylor展开都更加高效，且易编程实现.

Bengt Fornberg 计算有限差分格式权重的主要思想是利用插值多项式的导数来逼近原函数的导数.例如若已知n+1 个互异的插值节点x0,x1,…,xn的函数值f(xi)及直到mi阶的导数值(不一定需要使用到)信息，则存在插值多项式具有一般形式为

对pm(x)求q(q=0,1,…)阶导后，再代入逼近节点，可得到

显然，显式的有限差分格式权重，即ci，k可直接从上式右端获取.另一方面，如对求导后的式子分别代入所需点再进行适当组合可获得紧致差分格式的一般形式为

在此启发下，给出三种新的算法来求解在任意网格间距下的有限差分格式权重，它们分别基于Newton型插值多项式、带有连续阶数导数信息的Lagrange型Hermite插值多项式与带有连续阶数导数信息的Newton 型Hermite 插值多项式，因此分别称为Newton 型算法，Lagrange 型Hermite 算法及Newton 型Hermite算法.特别的，后面两种算法比文献[12]中的算法适用的有限差分公式类型更广.

1 Newton型算法

若已知f(x)在n+1个互异节点x0,x1,…,xn的函数值，则Newton型插值多项式[13]可表示为

其中：f[x0] =f(x0) 称为f(x) 关于节点x0的零阶差商.一般地，f[x0,x1,…,xk-1,xk] =称为f(x)关于节点x0,x1,…,xk的k阶差商.

观察式(1)可发现，其每一项都可写为差商f[x0,x1,…,xk]与乘积相乘的形式.不妨将乘积部分记为M0(x) =1,Mk(x) =(x-x0)(x-x1)…(x-xk-1),k=1,2,…,n, 显然乘积部分为k次多项式且存在递推关系Mk(x) =Mk-1(x)(x-xk-1).由此可将Newton型插值多项式重新表示为

对式(2)求q(q=0,1,…)阶导，则有从而为得到所需权重，只需分别考虑差商部分各点函数值f(xi)前的系数以及乘积部分求导后的结果．

首先考虑差商部分各点函数值f(xi)前的系数lk，i, 其中k代表k阶差商部分共涉及k+1 个互异的节点，从而k为整数且i=0,1,…,k.当k=0时，有f[x0] =f(x0), 即l0,0=1.当k＞0时，由差商可由函数值线性组合表出的性质[13]，可知

即在k阶差商f[x0,x1,…,xk]中，f(xi)前的系数为

接下来考虑对乘积部分Mk(x)求q阶导．若q=0, 即不求导的情况，Mk(0)(x) =Mk(x).若q＞0， 结合乘积部分递推关系式以及Leibniz高阶求导公式，则计算结果可分为以下三种情况，即

由Newton型插值多项式求导后的一般形式，即式(3)，以及上述差商部分和乘积部分的结论，即式(5)～(6)，可得q阶导数有限差分公式中f(xi)前的系数为

2 Lagrange型Hermite算法

若已知n+1个互异节点xi(i=0,1,…,n)上的函数值和直到mi阶的导数值，即则存在m(这里)次的Lagrange型Hermite插值多项式[13]，则有

根据Leibniz高阶求导公式，对式(9)求q(q=0,1,…)阶导有

下面分别对求和部分sk，i(x)与求积部分πi(x)的q阶导数情况进行考虑．

首先来看求和部分的导数(x).若q=0, 即在不求导的情况下，(x) =sk，i(x).若q＞0, 因为sk，i(x)是最低次数为k次且最高次数为mi次的多项式，所以求导结果可分为以下三种情况，即

接下来考虑求积部分的导数πi(q)(x).记hi(x) =(x-xi)mi+1, 根据多个函数相乘的高阶求导公式[14]有

由Lagrange 型Hermite 插值多项式求导后的一般形式，即式(10)，以及上述求和部分与求积部分的结论，即式(11)～(12)，可得q阶导数有限差分公式中f(k)(xi)前的系数为

3 Newton型Hermite算法

若已知n+1 个互异节点xi(i=0,1,…,n) 上的函数值和直到mi阶的导数值，即f(xi),f′(xi),…,f(mi)(xi),可得到次的Newton型Hermite插值多项式[13].则有

对式(15)求q(q=0,1,…)阶导，有

观察式(16)结构发现获得有限差分格式权重的关键在于分析每项中差商部分f(k)(xi)(k=1,2,…,mi)前的系数与乘积部分求导后的结果，下面分别对这两部分进行分析．

先来考虑差商部分f(k)(xi)前的系数．将差商中各节点重数之和记为μ, 则

当1 ≤μ≤m0+1时，即仅考虑节点x0的差商，此时μ=k+1.由k阶导数与k阶差商之间的关系，可得

显然，在f([x0,k+1])中，f(k)(x0)前的系数为

此时存在互异节点x0,x1,…,xj(j=1,2,…,n)的差商公式[15]可表示为

接下来考虑乘积部分求导后的结果.不妨令乘积部分为(x),μ(0 ≤μ≤m)表示某一乘积部分中多项式的次数，则乘积部分的通项可表示为

观察此可发现，乘积部分存在如下的递推关系

接下来对(x)求q(q=0,1,…)阶导数.若q=0,即不求导的情况，若q＞0, 结合乘积部分递推关系式(22)和Leibniz高阶求导公式，求导结果共分为以下两种情况，即有

由Newton 型Hermite 插值多项式求导后的一般形式，即式(16)，以及上述差商部分和乘积部分的结论，即式(18)、式(20)和式(23)，可得q阶导数有限差分公式中f(k)(xi)前的系数为

4 数值实验

下面选取三个具有代表性的算例来说明三种算法如何确定显式有限差分格式的权重、隐式紧致有限差分格式的权重以及带有连续阶导数信息的有限差分格式的权重．

首先，对于显式有限差分格式的权重的求解，以二阶导数中心差分格式为例，不妨选取节点xi=-1,0,1, 则二阶导数中心差分格式表现形式为

于是可在节点处利用Newton型算法以获取权重c(02),c(12),c(22)来近似导数f0′′, 即有

由式(7)知，所需信息有插值节点、逼近节点以及求导阶数.针对二阶导数中心差分格式，插值节点为xi=-1,0,1, 那么逼近节点为=0, 求导阶数q=2, 将上述相关数据应用到Newton型算法中，可计算出权重为

显然代入权重后的式(26)即为二阶导数的中心差分格式.利用文献[9]中算法可得到相同的权重结果．

特别地，若使用Lagrange 型Hermite 算法或Newton 型Hermite 算法求解二阶导数中心差分格式的权重，根据式(13)和式(24)可知所需信息与Newton型算法相比，除插值节点、逼近节点以及求导阶数以外，还需已知各插值节点所对应最高阶导数的阶数.通过式(26)知在计算二阶导数中心差分格式的权重中，只涉及到函数值，没有涉及到导数值，因此各插值节点所对应最高阶导数的阶数均为零，即mi=0,0,0.这表明，Lagrange型Hermite算法和Newton型Hermite算法是Newton型算法的推广.

其次，考虑隐式紧致有限差分格式权重的求解，以一阶导数的三点四阶紧致差分格式为例．这里同样选取节点xi=-1,0,1, 则一阶导数的三点四阶紧致差分格式表现形式为

为得到此格式，可先考虑如下式子

为能够使用Lagrange型Hermite算法或Newton型Hermite算法来确定式(29)的权重，将上式整理成

这样问题就转化为用点xi=-1,0,1 处的相应函数值及一阶导数值的加权和来近似五阶导数f0(5).此时，插值节点xi=-1,0,1, 逼近节点=0, 求导阶数q=5,各插值节点所对应最高阶导数的阶数mi=1,1,1, 将上述相关数据代入Lagrange型Hermite算法或Newton型Hermite算法中可计算出

利用文献[7]中的算法也可得出一致的权重结果.接着将权重代回式(30)后，等式两端同乘h4/180并舍去截断误差h4f(5)0, 即可得一阶导数的三点四阶紧致差分格式.

最后，考虑文献[16]中的公式

式(32)涉及到函数值、一阶导数值、二阶导数值和三阶导数值，是没有办法用文献[12]中的方法得到的.而根据Lagrange型Hermite算法和Newton型Hermite算法，可以得到

5 结论

Newton型算法、Lagrange型Hermite算法及Newton型Hermite算法均可通过计算机编程实现，免去大量的手动计算， 从而可快速得到所需权重．其中Newton 型算法与文献[9]中的方法相比，每新增一个节点，Newton 型算法计算工作不必重头开始，计算量少些(但鉴于如今的计算机速度，在得到权重的时间上已经看不到两者区别)． Lagrange 型Hermite 算法和Newton 型Hermite 算法与文献[12]中的方法相比，它们除可以用于只含有一阶导数的情况，还适用于出现多个不同阶导数的情况．