微分求积法(Differential Quadrature)学习笔记

1.DQ方法的基本原理1

微分求积法(DQ)是Bellman和他的同事在70年代初提出的。DQ方法是求导数近似的一种数值离散化方法,起源于传统积分求积的思想。

1.1 积分求积

一般地,积分∫abf(x)dx可以近似为:

(1.1)∫abf(x)dx=w1f1+w2f2+⋯+wnfn=∑k=1nwkfk

其中w1,w2,…,wn为权值系数,f1,f2,…,fn为函数在离散点a=x1,x2,…,xn=b处的值。方程(1.1)称为积分求积,它使用整个积分域中的所有函数值来近似一个有限积分上的积分。

1.2 微分求积

函数f(x)对 x在网格点xi处的一阶导数可以近似为整个域中所有函数值的线性和:

(1.2)fx(xi)=dfdx|xi=∑j=1Naij⋅f(xj),   i=1,2,…,N

其中aij 表示加权系数,N表示整个域中的网格点数。方程(1.2)称为微分求积(DQ).

2.基于多项式的微分求积法(PDQ)

在DQ近似中,如何确定加权系数很多工作都是基于多项式逼近的,因此相关的方法可以被称为基于多项式的微分求积法(PDQ)。

2.1 一阶导数加权系数的计算

Bellman 等人假定函数f(x)在区间[a,b]上足够光滑,因此它在任何网格点上的一阶导数 f(1)(x)可以用下面的公式来近似:

(2.1)fx(1)(xi)=dfdx|xi=∑j=1Naij⋅f(xj),   i=1,2,…,N

方程(2.1)中加权系数aij的确定是DQ近似的关键步骤。下面是一些方法

2.1.1 Bellman方法2

(2.2)gk(x)=xk,  k=0,1,…,N−1

显然有N个测试函数,那么需要N个网格点x1,x2,…,xN,可以得到关于aij的N×N个线性代数方程如下:

(2.3){∑j=1Naij=0∑j=1Naij⋅xj=1∑j=1Naij⋅xjk=k⋅xik−1,k=2,3,…,N−1    i=1,2,⋯,N.

由k=0,k=1,k=2,3,…,N−1分别代入得到。

方程组2.3有唯一解,因为它的矩阵是范德蒙形式:

V=[111⋯1x1x2x3⋯xNx12x22x32⋯xN2⋮⋮⋮⋱⋮x1N−1x2N−1x3N−1⋯xNN−1]

但是,当n很大的时候,矩阵是病态的,求解它的逆是困难的。在这种方法的实际应用中,通常选择n<13。

gk′(xi)=∑j=1Naijgk(xj)=aij,k=j,gk′(x)=LN(1)(x)(x−xk)−LN(x)(x−xk)2LN(1)(xk),gk′(xi)=LN(1)(xi)(xi−xj)LN(1)(xj)

 k=i时,根据L‘Hospital规则以及PN(x)满足的条件

(1−x2)LN″(x)−2xLN′(x)+N(N+1)LN(x)=0

得到:

aii=limx→xiLN(1)(x)(x−xi)−LN(x)(x−xi)2LN(1)(xi)=limx→xiLN(2)(x)(x−xi)+LN(1)(x)−LN(1)(x)2(x−xi)LN(1)(xi)=limx→xiLN(2)(x)2LN(1)(xi)=limx→xi2xLN′(x)−N(N+1)LN(x)(1−x2)LN(1)(xi)=2xi2(1−xi2)

这种方法不像第一种方法那样灵活,因为这种方法中网格点的坐标不能任意选择,所以在实际应用中通常采用第一种方法。

2.1.2 Quan and Chang’s 方法3

Quan and Chang用下面的拉格朗日插值多项式作为检验函数:

(2.6)gk(x)=M(x)(x−xk)⋅M(1)(xk),  k=1,2,…,N

其中:

(2.7)M(x)=(x−x1)(x−x2)…(x−xN)
(2.8)M(1)(xi)=∏k=1,k≠iN(xi−xk)

随后,通过在N个网格点上应用方程(2.6),得到了以下计算加权系数aij的代数公式:

(2.9a)aij=1xj−xi∏k=1,k≠i,jNxi−xkxj−xk, for j≠i
(2.9b)aii=∑k=1,k≠iN1xi−xk

显然可以得到gk(xj)=δjk

gk′(xi)=∑j=1Naijgk(xj)=aikgk′(x)=M′(x)(x−xk)−M(x)(x−xk)2⋅M(1)(xk)

当i≠j时

∴  gk′(xi)=M′(xi)(xi−xk)⋅M(1)(xk)=M(1)(xi)(xi−xk)⋅M(1)(xk)∴  aik=M(1)(xi)(xi−xk)⋅M(1)(xk)=1(xi−xk)⋅∏m=1,m≠iN(xi−xm)∏m=1,m≠kN(xk−xm)∴  aij=1xj−xi∏k=1,k≠i,jNxi−xkxj−xk

当i=j 时

aii=limx→xiM′(x)(x−xi)−M(x)(x−xi)2M(1)(xi)=limx→xiM″(x)2M(1)(xi)=limx→xi2∑m=1,m≠iN∏k=1,k≠i,mN(x−xk)2∏k=1,k≠iN(xi−xk)=∑k=1,k≠iN1xi−xk

另一种方法(也就是下面的(???)):令M(x)=N(x,xk)⋅(x−xk),k=1,2,…,N,这里N(xi,xj)=M(1)(xi)⋅δij.得到:

gk(x)=N(x,xk)M(1)(xk),k=1,2,…,N∴aij=N(1)(xi,xj)M(1)(xj)

又有

(2.10)M(m)(x)=N(m)(x,xk)⋅(x−xk)+m⋅N(m−1)(x,xk),m=1,2,…,N−1;k=1,2,…,N
∴N(1)(xi,xj)=M(1)(xi)xi−xj,i≠jN(1)(xi,xi)=M(2)(xi)2∴aij=M(1)(xi)(xi−xj)M(1)(xj),i≠jaii=M(2)(xi)2M(1)(xi)

2.1.3 Shu’s Approach4

偏微分方程的解可以用高次多项式精确逼近。现在。我们假设逼近多项式的次数是N−1,这个近似多项式通过向量加法和标量乘法运算构成了一个N维线性向量空间,可以用不同的形式表达。基多项式的四个典型集如下:

(2.11a)rk(x)=xk−1,k=1,2,…,N
(2.11b)rk(x)=LN(x)(x−xk)⋅LN(1)(xk),k=1,2,⋯,N
(2.11c)rk(x)=M(x)(x−xk)⋅M(1)(xk),k=1,2,⋯,N
(2.11d)r1(x)=1,rk(x)=(x−xk−1)⋅rk−1(x),k=2,3,…,N

其中LN(x)是N次的Legendre多项式,m(x)在方程(2.7)中定义。在四组基多项式中。方程(2.11b) 和(2.11c) 来自拉格朗日插值多项式,方程(2.11d) 来自牛顿插值多项式。方程(2.11b) 和(2.11c) 的区别在于网格点的分布。方程(2.11b)是方程(2.11c)的一个特例,因为它只在Legendre配点处有效。

此外,根据线性向量空间的性质,如果一组基多项式满足一个线性算子,那么另一组基多项式也满足。因此 aij满足以下方程,这个方程是由基多项式xk−1在k=1时得到的:

(2.12)∑j=1Naij=0或aii=−∑j=1,j≠iNaij

2.2 二阶导数加权系数的计算

对于二阶导数的离散化,引入了一个类似的近似形式:

(2.13)fx(2)(xi)=∑j=1Nbij⋅f(xj),   i=1,2,…,N

2.2.1 Quan and Chang’s 方法

(2.14a)bij=2xj−xi(∏k=1,k≠i,jNxi−xkxj−xk)(∑l=1,l≠i,jN1xi−xl),fori≠j
(2.14b)bii=2∑k=1,k≠iN−1[1xi−xk(∑l=k+1,l≠iN1xi−xl)]
易得 bij=gk″(xi)g″(x)=M″(x)(x−xk)2−2M′(x)(x−xk)+M(x)(x−xk)3M(1)(xk)∴gk″(xi)=M″(xi)(xi−xk)−2M(1)(xi)(xi−xk)2M(1)(xk)∴i≠j时,aij=gk″(xi)=2∑k=1,k≠iNM(1)(xi)/(xi−xk)(xi−xk)M(1)(xk)=2xj−xi(∏k=1,k≠i,jNxi−xkxj−xk)(∑l=1,l≠i,jN1xi−xl)

2.2.2 Shu‘s 方法

∵gk(x)=N(x,xk)M(1)(xk)∴bij=N(2)(xi,xj)M(1)(xj)又有N(2)(xi,xj)=M(2)(xi)−2N(1)(xi,xj)xi−xj,i≠jN(2)(xi,xi)=M(3)(xi)3
(2.15a)bij=M(2)(xi)−2N(1)(xi,xj)(xi−xj)M(1)(xj),i≠j,
(2.15b)bii=M(3)(xj)3M(1)(xj)

由aij,aii的定义可得:

(2.16a)bij=2aij(aii−1xi−xj),i≠j

同样,也可以根据基函数为xk−1时的定义得到:

(2.16b)∑j=1Nbij=0或bii=−∑j=1,j≠iNbij

2.3 Shu’s方法关于高阶导数的递推公式

(2.17)fx(m−1)(xi)=∑j=1Nwij(m−1)⋅f(xj),
(2.18)fx(m)(xi)=∑j=1Nwij(m)⋅f(xj)

其中i=1,2,…,N;m=2,3,…,N−1。同样,我们可以得到:

(2.19)wij(m−1)=N(m−1)(xi,xj)M(1)(xj)
(2.20)wij(m)=N(m)(xi,xj)M(1)(xj)

所以得到:

(2.21)N(m−1)(xi,xj)=wij(m−1)M(1)(xj),i≠j

通过(2.10),得到:

(2.22)N(m−1)(xi,xi)=M(m)(xi)m
(2.23)N(m)(xi,xj)=M(m)(xi)−mN(m−1)(xi,xj)xi−xj,i≠jN(m)(xi,xi)=M(m+1)(xi)m+1

所以:

(2.24)N(m)(xi,xj)=m[N(m)(xi,xj)−N(m−1)(xi,xj)]xi−xj,i≠j
(2.25)∴   N(m)(xi,xj)=m[wii(m−1)M(1)(xi)−wij(m−1)M(1)(xj)xi,xj,i≠j

从而有:

(2.26)wij(m)=m(aijwii(m−1)−wij(m−1)xi−xj),i,j=1,2,…,N;m=2,3,…,N−1

同样,i=j时的权值系数如下:

(2.27)∑j=1Nwij(m)=0或wii(m)=−∑j=1,j≠iNwij(m)
注意:gk″(xi)=∑j=1Nbijgk(xj)=∑m=1Naim∑j=1Namjgk(xj)∴bij=∑m=1Naimamj

 

3.广义积分求积法(GIQ)

广义积分求积(GIQ)是在与PDQ相同的概念下发展起来的。如果一个函数在整个区域内是光滑的,那么它可以用一个涉及整个区域内所有函数值的高次多项式来逼近。然后,通过对近似多项式进行积分,可以计算出函数在整个区域的一部分上的积分。因此,这种近似包含了整个区域中的所有函数值,而且精度很高,即使积分区域只包含两个点。显然,该方法的关键步骤是确定权重系数。下面将说明,GIQ的加权系数可以很容易地从PDQ中的一阶导数的加权系数中得到。

3.1 1维广义积分求积

假设f(x)在整个区域[a,b]上是连续的,且可以分解为N−1个区间,其网格点为a=x1,x2,…,xN=b.由于f(x)在整个区域上是连续的,所以它可以用一个(N−1)次多项式来逼近.特别地,当N个网格点的函数值已知时,f(x)可以用与所有网格点的函数值相关的拉格朗日插值多项式来逼近。因此,该逼近多项式在[xi,xj]上的积分可能涉及整环以外的函数值。作为一般情况,假设整个域的一部分上的f(x)的积分由整个域中的所有函数值与形式的线性组合来近似:

(3.1)∫xixjf(x)⋅dx=∑k=1Nckij⋅f(xk)

与PDQ中的分析类似,作为f(x)的近似的(N−1)次多项式构成了N维线性向量空间。因此,如果所有的基多项式都满足等式(3.1),那么空间中的任何多项式都满足。考虑拉格朗日插值多项式rk(x),k=1,2,…,N作为基函数,那么

(3.2)ckij=∫xixjrk(x)⋅dx

直接计算ckij的表达是困难的,所以采用如下方法:

(3.3)f(x)=du(x)dx

我们可以清楚地看到,如果f(x)是(N−1)次多项式,则u(x)应该是N次多项式。假设f(x)由以下形式的(N−1)次多项式逼近:

(3.4)f(x)=a0+a1x+a2x2+⋯+aN−1xN−1

其中a0,a1,aN−1为常数,对(3.4)从常数c到变量x作积分有:

(3.5)u(x)=∫cxf(t)dt+u(c)=F(x,c)+u(c)

其中

F(x,c)=x(a0+a12x+⋯+aN−1NxN−1)−c(a0+a12c+⋯+aN−1NcN−1)

显然,F(x,c)取决于N个常数。证明了F(x,c)构成N维线性多项式向量空间。它的一组基多项式可选为:

(3.6)pk(x)=(x−c)rk(x),k=1,2,…,N

相似于PDQ,可以得到:

(3.7)Fx(xi,c)=∑j=1Na―ijF(xj,c),i=1,2,…,N

将(3.6)代入(3.7)中得到:

(3.8a)a―ij=xi−cxj−cwij(1),j≠i
(3.8b)a―ii=wii(1)+1xi−c,
pk′(x)=rk(x)+(x−c)rk′(x)∴  (xi−c)rk′(xi)=a―ij(xj−c)       (i≠j)∴  a―ij=xi−cxj−crk′(xi)=xi−cxj−cwij(1)同样pk′(xi)=1+(xi−c)wii(1)            (k=i)∴a―ii=wii(1)+1xi−c

从公式(3.8a)−(3.8b)中可以看出,不能选择c作为网格点的坐标,从公式10.6可以看出,当c=x时,积分区域为零。因此,不必近似c=x的积分。另一方面,公式(3.8a)−(3.8b)可以写成向量形式:

(3.9){Fx}={A―}{F}

其中

{F}={F(x1,c),F(x2,c),…,F(xN,c)}T{Fx}={f}={Fx(x1,c),Fx(x2,c),…,Fx(xN,c)}T

将(3.4)和(3.6)插入(3.9)中,令

(3.10){fI}={F}={∫cx1f(x)dx,∫cx2f(x)dx,…,∫cxNf(x)dx}T

可以得到:

(3.11){f}=[A―]{fI}

令[WI]=[A―]−1,则有{fI}=WI{f},进一步:

(3.12)∫cxif(x)dx=∑k=1NwikIf(xk),   i=1,2,…,N

那么就得到

(3.13)ckij=wjkI−wikI

很明显GIQ中的权重系数是由PDQ中的一阶导数离散化得到的。

参考文献

[1] Shu, Chang. Differential Quadrature and Its Application in Engineering[J]

[2] Bellman R, Kashef BG, Casti J. Differential Quadrature and Its Application in Engineering[J]

[3] Quan JR, Chang CT. New insights in solving distributed system equations by the quadrature method—I. Analysis[J]

[4] Shu C, Chew YT, Richards BE. Generalized differential and integral quadrature and their application to solve boundary layer equations[J]