(四)高精度差分格式及其数值解的逼近程度分析

2
(1) (2)
−π
K=1
π
K=2
−π ≤ x ≤ π
模型方程及半离散方程(cont.)
半离散方程（ODE方程） Fj ∂u 1. + c =0 (3) ∂t ∆x F S ∂u 2. + c j = µ j2 (4) ∂t ∆x ∆x 其中 ∂u S j ∂ 2u 1 → , 2 → 2 ，例如F j = (u j +1 − u j −1 ), S j = u j +1 − 2u j + u j −1 ∆x ∂x ∆x ∂x 2 离散的初值： u ( x j ,0) = exp(ikx j ) Fj 离散方程的精确解及其导数(无∆x) : ) u ( x j , t ) = u (t ) exp(ikx j ) ) ) F j = ke ⋅ u (t ) exp(ikx j ) ) ) S j = k d ⋅ u (t ) exp(ikx j ) 代入（3），（4）
(
)
要求数值格式： ) kd →1 α2 ) ) ) kr → 0（0耗散要求；k r > 0表示格式有正耗散，k r < 0负耗散） α ) ki → 1 (对 色散误差，要求保持原波速c）
α
4
4.2 高精度差分格式
• 特别适用于光滑问题 • 减少网格点数 • 应用于湍流等多尺度问题
5
4.2.1 传统型差分格式
(四) 高精度差分格式及其数值解 四 的逼近程度分析
• 指大于二阶精度的格式 • 要求准确模拟小扰动量长时间、远距离 传播的速度和幅值 • 用于计算噪声、DNS等 • 有传统高精度差分和紧致差分两种
1
4.1模型方程及半离散方程
模型方程 ∂u ∂u 1. + c =0 ∂t ∂x ∂u ∂u ∂ 2u 2. + c =µ 2 ∂t ∂x ∂x c = const , µ = const 单个Fourier分量初值 : u ( x,0) = exp(ikx), 精确解： 1. u ( x, t ) = exp[ik ( x − ct )] 2. u ( x, t ) = exp(− µk 2t ) exp[ik ( x − ct )]
• 三阶迎风紧致格式对模型方程的逼近程 度：
ke = k r + iki (1 − cos α ) 2 sin α (8 + cos α ) kr = , ki = 5 + 4 cos α 5 + 4 cos α k r > 0，耗散型格式
22
4.3 格式的精度及分辨率
l ∂u ∂f l −1 ∂ f + = ∑ al ∆x ∂t ∂x ∂x l l
∂ 2 m −1 f 2 m −1 ∂ ∂f ∂u ∂ ∂u ∂f 2 m −1 ∂x = ⋅ ⋅ + + ∑ µ 2m ∆x ∑ µ 2 m+1∆x 2m ∂f ∂u ∂x ∂x m ∂t ∂x ∂x m ∂x = ∂ ~ µ2 ∂x ∂ ~ µ2 ∂x ~ ∂u µ ∂u ∂ ~ ∂u ∂ ~ ∂u ∂ 3 ∂x ∂ 2u + µ3 = µ2 + ⋅ ∂x ∂x ∂x ∂x ∂x ∂x ∂ 2u ∂x 2 ∂x 2 ∂u ∂ ∂ 2u + a ∂x ∂x ∂x 2
• 获得S需要求解矩阵方程组
13
二阶导数的紧致差分(cont.)
• 5点六阶紧致：
5 2 δ x S j = δ x2 (u j + δ x2u j ) 12 对角不占优 ? Sj +
14
一阶导数的对称型紧致差分
• 一般形式:
2 相容要求：
l ≥0 l ≥0 l l >1 l
∑α
F j + l + F j −l
• 六阶精度:
1 2 1 1 1 1 F j −1 + F j + F j +1 + δ x2 F j = δ x0u j − δ x0u j − (u j + 2 − u j − 2 ) 6 3 6 30 15 4
17
一阶导数的迎风紧致差分
• 有数值耗散，抑制高频振荡 • 一般形式： ∑α F
(2) (3)
J-3
J-2 J-1 J
J+1 J+2
J-2
J-1 J
J+1
J+2 J+3
7
传统型差分格式(cont)
迎风偏斜差分的应用： ∂u ∂u +c =0 ∂t ∂x Fj ∂u +c =0 ∂t ∆x cF j = c + F j+ + c − F j−
8
传统型差分对模型方程的逼近特性
• 所有中心型差分：无耗散 • 迎风偏斜差分：有耗散，但可能为负耗散 • 可从差分格式的精确解分析色散和耗散
• 若边界点S0和SN已知，可用解三对角矩阵方程 得到所有网格点上的差分
12
二阶导数的紧致差分(cont.)
• 一般形式
2 l ≥0 相容要求：
∑β
S j + l + S j −l
l
= ∑ bl
l >0
u j +l − 2u j + u j −l l2
∑ β = ∑b
l ≥0 l l >1
l
∂ 2k S ∂ 2 k + 2u Taylor展开左端的 2 k 和右端的∆x 2 2 k + 2 的系数相等 ∂x ∂x 得到代数方程组，求出系数β l 和bl .
(
)
(
)
(
)
10
4.2.2 紧致差分格式
• • • • • 同样的精度比传统差分的基架点少 截断误差的系数较小 二阶导数的(中心型)紧致差分 一阶导数的对称紧致差分 一阶导数的迎风紧致差分
11
二阶导数的紧致差分
• 传统型差分的截断误差项的再次离散 • 四阶紧致：
1 5 1 S j −1 + S j + S j +1 = δ x2u j 12 6 12
首项既是精度：l < l0 , al = 0, al0 ≠ 0, 精度为l0 − 1阶。 奇数阶导数和偶数阶导数分开：
2m 2 m +1 ∂u ∂f f f 2 m −1 ∂ 2m ∂ + = ∑ µ 2 m ∆x + ∑ µ 2 m +1∆x ∂t ∂x m ∂x 2 m ∂x 2 m +1 m 改写成：
逼近于一阶导数的四阶精度中心型差分 1 Fj = 8(u j +1−u j −1) ) − (u j + 2 −u j − 2 ) ) 12 五阶精度的迎风偏斜差分
[
]
(1)
F j+ =
1 − δ x (− 3u j + 2 + 27u j +1 + 47u j − 13u j −1 + 2u j − 2 ) 60 1 F j− = δ x+ (− 3u j − 2 + 27u j −1 + 47u j − 13u j +1 + 2u j + 2 ) 60
l
= ∑ al
l >0
u j + l − u j −l 2l
∑α = ∑ a (= 1)
∂k S ∂ k +1u 于j处Taylor展开左端的 k 和右端的∆x k +1 的系数相等 ∂x ∂x 得到代数方程组，求出系数α l 和al .
• 获得F需要求解矩阵方程组，左端点最好不超过3
15
一阶导数的对称型紧致差分(cont.)
18
一阶导数的迎风紧致差分(cont.)
• 五阶迎风紧致格式(5点）：
c ≥ 0: 3 + 2 + 1 − F j + F j −1 = δ x (− u j + 2 + 11u j +1 + 47u j + 3u j −1 ) 5 5 60 c < 0: 3 − 2 − 1 + F j + F j +1 = δ x (− u j − 2 + 11u j −1 + 47u j + 3u j +1 ) 5 5 60
微分方程精确解及其导数： ˆ u ( x, t ) = u (t ) exp(ikx) ∆xu = ik∆x ⋅ u (t ) exp(ikx) = k ⋅ u (t ) exp(ikx), k = i (k∆x) = iα ˆ ˆ x e e ∆x 2 u = −k 2 ∆x 2 ⋅ u (t ) exp(ikx) = −k ⋅ u (t ) exp(ikx), k = α 2 ˆ ˆ xx d d
• F+按j增加的方向求解, F-按j减少的方向求解
19
非等距网格上的格式
∂f ∂f • 计算精度比 = ⋅ξ x 要高。 ∂x ∂ξ • 三阶迎风紧致：
α j F j+−1 + F j+ = a j
f j +1 − f j hj + bj f j − f j −1 h j −1
θ j2 2 + 3θ j 1 αj = ,aj = ,bj = 1+θ j (1 + θ j ) 2 (1 + θ j ) 2 γ j F j−+1 + F j− = a j γj =
直接由网格点的函数值线性组合 N阶精度至少需要N+1个网格点 Taylor展开确定系数
∂u ∆x ∂x F j = ∑ al (u j +l +1 − u j +l ) Fj →