微分求积法

X ij = W j ( x) |x = xi = l j
[k ]
[k ]
[k ]
( xi )
[ − ∑ lm ] ( xi )
k m =1 m≠i
N
(i ≠ j )
式（1.19）对微分求积法适用，这样，当一阶导数的权系数确定后，求 高阶导数的权系数有两种方法可供选择：既可用递推公式（1.5），也 可以用公式（1.19）。
以采用相同的权函数。）
f n ( xi ) = fi n = ∑ X n（xi ) f j = ∑ X n ij f j j
j =1 j =1
（1.4）
我们称 X ij 为第i个节点的第j个n阶加权系数。这样N个节 点共有 N 2 个加权系数。下面推导各阶加权系数的递推关 系。（黑板上推）
n
X
n +1 ij
∑A
k =1
N
y jk
A y kn
A y ij
N N ∏ ( yi − y k ) ∏ ( y j − y k ) (i ≠ j ) k =1 k =1, j k≠ j k ≠i ' = l j ( yi ) = N 1 ∏ ( y − y ) (i = j ) k =1 j k k ≠i
−1
（1.22）
(k − 1) x
k −2 i
= ∑ Aij x k −1 j
n =1
N
0 0 G] = [ ⋮ 0
1 2 x1 1 2 x2 ⋮ ⋮ 1 2 xN
⋯ ( N − 1) x1N − 2 N ⋯ ( N − 1) x2 − 2 ⋱ ⋮ N ⋯ ( N − 1) xN − 2 N × N
dξ β − α = ，由式（1.14）可知 dx a − b
l j ( x ) = l，因此 j (ξ )
l 'j ( x ) = l 'j (ξ )
dξ β − α ' = ll (ξ ) dx b − a
2
β − α '' d ξ β − α '' l ''j ( x ) = = lj l j (ξ ) b − a dx b − a
（1.10）
M M ∂kg k [k ] [k] = h i p j = ∑ Ψ jm p m h i = ∑ Ψ [jm]g m ∂y k ij m =1 m =1
N M N M ∂ r+s g s] [ r ] [s ] [r] [s ] [r] = h i p j = ∑ Φ il h l ∑ Ψ jm p m = ∑ Φ il ∑ Ψ [jm g lm ∂x r ∂ys ij l =1 m =1 l =1 m =1
若采用均匀网格划分，横向M=100个节点，纵向N=50个节 点。节点间距h=0.02。则共有5000个节点。
T ( x, y ) = X ( x)Y ( y )
M ∂ 2T '' = X ( x)Y ( y ) = ∑ Bm ( x) X mY ( y ) 2 ∂x m =1 N ∂ 2T '' = X ( x)Y ( y ) = ∑ Bn ( y )Yn X ( x) ∂y 2 n =1
（1.26）
⋮
β − α [k ] l j ( x) = l j (ξ ) b−a
[k ]
k
当节点 ξ i 对应于正则化区域，即 ξi ∈ [ −1,1] ，时，式（1.25）简化 为
b−a b+a xi = ξi + 2 2
（1.27）
对于任意一维有限区域[a,b],只需要按上式划分N个互异节点 ，就可以通过公式（1.26）将正 则化区域的权系数转换为当前要求的权系数。 权系数可以转换的特性，给编程带来很大的方便。通常情况下，我们只 要预先计算一次正则化区域的权系数，并予以保存。在求解具体问题时， 只须将已有的权系数按式（1.26）进行变换，而无须重新计算。
N
(i = j )
根据插值函数 l j ( x )各阶导数的递推关系，可以得到二阶及二阶以上的导 数的权系数显示计算表达式
k −1 k −1 l [j ] ( x j ) k li[ ] ( xi ) l 'j ( xi ) xi − x j
(i = j )
（1.19）
对应于节点（ xi , y j ） ，有：
∂ 2Tij ∂x ∂y T
'' ij 2
= ∑ Bim X mY j = ∑ BimTmj
m =1 N m =1 N
M
M
∂ 2Tij
2
= ∑ B j nYn X i = ∑ B j nTin
n =1 n =1 M N
= ∑ Bi mTmj + ∑ B j nTin
(2i − 1)π xi = − cos 2N ( N − i )π xi = cos ( N − 1)
xi = ( N − 2)个Gauss-Legendre积分点
xi = − cos
( 2i − 3) π 2 ( N − 2)
4、算例
对于二维温场分布所满足的泊松方程
∇ 2T ( x, y ) = xy x ∈ [−1,1], y ∈ [−0.5, 0.5] 矩形边界区域T=0。
j =1
N
（1.16）
对比可得：
Aij = l 'j ( xi )
因此：
k =1 k ≠i , j
（1.17）
∏ ( x − x ) ∏( x
N N i k k =1 k≠ j
j
− xk )
(i ≠ j )
（1.18）
Aij = l ( xi ) =
' j
1 ∏ (x − x ) k =1 j k
k ≠i
, y) 对于二维函数 g ( x，通常可以用一维函数的积来表示，即
g ( x, y ) = h ( x ) p ( y )
（1.7）
两个一维函数在节点的k阶导数值可以由节点的函数值表 示：
hi = ∑ Φ hl
l =1
M
[k ]
N
[k ]
il
k
（1.8）
p j = ∑ψ [jm]p m
m =1
[k ]
（1.11）
（1.12） 1.12
各阶偏导数的加权系数都出自一维函数的加权系数， 因此我们对一维函数的加权系数进行详细地讨论。
2.权系数的确定
对于n个离散点，我们可以用拉格朗日插值函数来近似代替他们所满 足的函数关系。 定义：
l j ( xi ) = δ ij f ( x) = ∑ l j ( x) f ( x j )
（1.9）
式（1.8）和（1.9）中的 N，M分别为x和y方向划分的节点数。 由上述两式可以导出二维函数g（x，y）在节点 ( x i，y的各阶偏 j) 导数的求值公式：
N N ∂kg k [k ] [k ] = h i p j = ∑ Φ il h l p j = ∑ Φ[il ]g lj ∂x k ij l =1 l =1
若
，并记
f (x) =
'
，f = f ( x ) j j
∑W
j =1
N
(注意函数f ( x)的任意性，这就要求
j
(x) f j
权函数W j(x)只与节点值xi 有关，而 与节点处的函数值无关，也就是说 对应于相同节点的不同函数我们可
（1.2） （1.3）
fi =
'
∑
N
N
j =1
Aij f j
同理：
N
节点的选取
用微分求积法求解问题时，计算精度将由权函数的选取决定，必须选 取合理分布的节点。针对具体问题求解区域性质的不同，常采用非均匀 节点划分，对靠近边界处和物理性质变化激烈的区域可以选取较多的结 点。 下面是常用的几种非均匀节点划分方式。
xi = 2 i −1 −1 N −1
（1.28a） （1.28b） （1.28c） （1.28d） （1.28e）
加权系数取决于基函数的选取以及节点的划 分
如果采用如下形式的插值基函数：
则（1.14）式写作：
采用这种形式的插值基函数的方法称为调和微分求积法。 由此可以导出一阶加权系数的表达式为：
Aij = l ' j ( xi )
高阶加权系数可有一阶系数由递推关系得到。
下面我们从另外一个角度来计算加权系数： 由微分求积法的定义（1.2）式，如果我们需要对未知函数求N阶导数， 我们需要保证其N-1阶导数连续，只要满足这个条件，我们可以选取任意 函数簇f（x）来导出权系数。
j =1 N
（1.13）
式中 l j ( x ) 称为拉格朗日插值多项式，其形式为
x − xk l j ( x) = ∏ k =1 x j − xk
N k≠ j
（1.14）
由式（1.13）对f（x）求一阶导数，得到
f
从而有
'
( x ) = ∑ l 'j ( x ) f j
j =1
N
（1.15）
fi ' = f ' ( xi ) = ∑ l 'j ( xi ) f j
A x ij
M M ∏ ( xi − x k ) ∏ ( x j − x k ) ( i ≠ j ) k =1 k =1, j k≠ j k ≠i ' = l j ( xi ) = M 1 ∏ ( x − x ) (i = j ) k =1 j k k ≠i
B jn =
（1.23）
1 x1 1 x2 [V ] = ⋮ ⋮ 1 xN
ห้องสมุดไป่ตู้
⋯ x1N −1 N ⋯ x2 −1 ⋱ ⋮ N ⋯ xN −1 N × N
（1.24）
3、不同求解域权系数的转换关系
实际应用中为了避免对不同的求解区域多次编程计算，我们可以采 用变量代换的方法将新变量转换为已有变量。 对于变量 x ∈ [ a, b ] ，与变量 ξ ∈ [α , β ] 常采用线性变换关系
f ( x ) = x k −1 ( k = 1, 2, ⋯, N )
(k − 1) x
k −2 i
(i, k = 1, 2,⋯ , N )
（1.20）
= ∑ Aij x k −1 j
n =1
N
将方程（1.20）写成矩阵形式，得到
[G ] = [ A][V ]
所以
（1.21）
[ A] = [G ][V ]
m =1 n =1
这样原方程化为如下代数方程组：
∑B
m =1
M
i m mj
T + ∑ B j nTin = xi y j
n =1
N
每个节点对应一条方程，共MN条方程。 对于边界上：
T1 j = Ti1 = TMj = TiN = 0
下面计算相应的二阶加权系数：
Bi m =
∑A
k =1
M
x ik
A x km
Aij = W j ( xi )
Bij = ∑ Aik Akj
k =1
N
N
（1.6）
N k =1
Cij = ∑ Aik Bkj = ∑ Bik AKj
k =1
N
Dij = ∑ Aik Ckj = ∑ Ckj Aik
k =1 k =1
N
可以看出，关键是确定一阶导数的加权系数，也就是插值 基函数 W j ( x) 。
微分求积法基本原理
• 目录
• • • • 微分求积法的定义 权系数的确定 不同求解域权系数的转换关系 一个算例
1.定义
不失一般性，考虑一维函数f(x),它在区间[a,b]上连续可微，则有 (1.1)
式中L是线性微分算子； 是插值基函数， 为N个互异节点 中第j节点的坐标值。
对于定义在D上的连续函数，将全域上各点的函数值加权 求和来表示各节点处的导数值的方法称为微分求积法。若 将函数值看做基函数，可以理解为对插值基函数加权求和。
x= b−a α +β ξ− 2 β −α a+b + 2
（1.25）
则区间[a,b]上N个互异节点 a = x1 < x2 < ⋯ < xN = b ，可以一一确 定出区间[α,β] 上对应的N个节点 α = ξ1 < ξ 2 < ⋯ < ξ N = β 。 由式（1.25）可得
= ∑ Aik X n kj
k =1
N
（1.5）
实际应用中，我们一般最高求导为四阶， 我们写出（1.4）式和（1.5）式的前面 四项：
N ' f i = ∑ Aij f j j =1 [ 2] N f i = ∑ Bij f j j =1 N [3] f = C f ∑ ij j i j =1 N f [ 4] = D f ∑ ij j i j =1