2.1有限差分法基础(课堂PPT)

对于函数f(x)，通常意义下的导数（微商）定义为：
f (x ຫໍສະໝຸດ Baidux) f (x)
x lim
dx
x
lim
f (x) f (x dx) dx
x
lim
f (x dx) f (x dx) 2dx
当dx→0时，以上三种形式都是微商的正确定义。 如果dx是有限的，如何给出微商的近似定义？
8
2.3 差分格式
L(U)f(P), PD I B (U )g(P ), P D B
其中：以上差分方程是偏微分方程的有限差分近似，U是u的有限差分近似。 差分方程要求U在所有节点上是u的很好的近似，并且方程所给出的有限差 分近似解U是唯一的。
25
差分方程 例子：对流方程（双曲型）的初值问题
u au 0, xR, t 0
(20)
如果u(x,t)是满足方程（13）的光滑解，则
u t
a
u x
n
j
0
代入（20），可以看出，偏微分方程（13）在(xj, tn)处可以近似地 用下面的方程来代替：
un1 j
unj
a
un j1
unj
0
h
(21)
j 0,1,2,, n=0,1,2,
其中是u(xj,tn)的近似值。（21）式称为逼近方程（13）的有限差分方程，简 称差分方程。用到的节点如图所示，
3. 建立和求解差商方程组
14
差分格式的另一种推导
(i-1,j+1) (i,j+1) 2 j+1
(i+1,j+1)
j (i-1,j)
j-1 (i-1,j-1)
h2
3
0
h3
h1
(i,j)
h4 4
(i,j)
1 (i+1,j)
(i+1,j-1)
i-1
i
i+1
为了寻求更精确的差分格式，我们引入两个 待定常数，由泰勒展开，构造如下关系式
该方法是一种直接将微分问题变为代数问题的近似数值解法，数学 概念直观，表达简单，是发展较早且比较成熟的数值方法。
4
2.1 有限差分法基础
有限差分方法的基本原理
该方法将求解域划分为差分网格，用有限个网格节点代替连续的求 解域。有限差分方法以Taylor级数展开等方法，把控制方程中的导数 用网格节点上的函数值的差商代替进行离散，从而建立以网格节点 上的函数值为未知数的代数方程组。
j (i-1,j)
j-1 (i-1,j-1)
3 h2 0
(hi3,j)
h1
h44
(i,j)
(i+1,j+1)
1 (i+1,j) (i+1,j-1)
i-1
i i+1
22
（2） 第二类和第三类边界条件
过O点向边界G做垂线PQ交边界于Q， 交网线段VR于P，OP=ah，PR=bh，VP=ch
对于点O 因为P一般不是节点，其值应当以点和P、R点的 插值给出
L(u(P ))f(P ), P D I 边界条件表示为以下方程：
B (u(P ))g(P ), P D B
设D
I
和D
B
分别表示区域D的北部节点和边界节点，则下式表示了以
上偏微分方程的有限差分方程(finite-difference equations, or finite-
difference scheme)：
df(x)f(xx)f(x)
dx
x
（5）
（5）式称为向前差分格式（forward-difference formula）
由（2）式得到
d ( x)2d2
( x)3d3
( x)4d4
f(x)f(x- x) x f(x)d x
2 !
d x2f(x)3 !
d x3f(x)-
4 !
d x4f(x)
df(x)f(x)f(x x)+ O ( x)
2.从理论上研究有限差分模型的形态，以保证计算过程的可行性和计算结果的正确性 （1）解的相容性； （2）解的稳定性； （3）解的收敛性。
3. 如何数值求解差分方程组
6
2.2 网格剖分
• 网格剖分就是研究区域和边界的离散化 • 1.矩形分割 • 2.三角形分割 • 3.极网格分割
7
对地球物理问题的连续求解区域通过网格划分离散为空间上得一系 列网格点，接下来需要利用一定的差分格式对偏微分方程组中的导 数用差商进行近似，从而将偏微分方程组离散化为差分方程组。
(11)
d2
f(x x)2f(x)f(x- x)
dx2f(x)
( x)2
(12)
（12）式称为二阶导数的二阶精度中心差分形式。忽略Δx的四次方及更高阶项
12
总结： 对一阶导数 1、向前差分形式： 2、向后差分形式： 3、中心差分形式：
对于二阶导数
df(x)f(xx)f(x)
dx
x
d f(x)f(x)f(x-x)
汶川大地震的动力学成因
Hu, C., Y. Cai, and Z. Wang (2012), Effects of large historical earthquakes, viscous relaxation, and tectonic loading
on the 2008 Wenchuan earthquake, Journal of Geophysical Research, 117, B06410, doi:10.1029/2011JB009046.
f(x x )f(x ) xdf(x ) ( x )2d 2f(x ) ( x )3d 3f(x ) ( x )4d 4f(x )
d x
2 ! d x2
3 ! d x 3
4 ! d x4
（1） （2）
（1）式减去（2）式，得到：
df(x)=f(x x)f(x- x)+ O ( x2)
d x
2 x
第二章 有限差分法
主讲人：胡才博
中国科学院大学地球科学学院 中国科学院计算地球动力学重点实
验室
1
第二章 有限差分法
• 2.1 有限差分法基础 • 2.2 网格剖分 • 2.3 差分格式 • 2.4 差分方程 • 2.5 应用实例
2
解析方法的局限性
1. 地球内部介质，不仅存在纵向非均匀结构（一维地球模型）， 也存在横向非均匀结构（不同块体、断层系统）； 2. 几何模型也呈现出相当的复杂性； 3. 另外，边界条件和初始条件对于不同问题具有特殊性。
dx
x
df(x)f(xx)f(x-x)
dx
2x
单侧，一阶精度 单侧，一阶精度 对称，二阶精度
d2 f(x)f(x x)2f(x)f(x- x)
dx2
( x)2
二阶精度
13
• 定解问题的有限差分解法 1．离散
x = ih, y= jh, i= 0, ±1, ±2,…. ±n, h: 步长（正方形的边长） 2．根据泰勒级数建立差商格式：对于一 维情况:在x处的一阶导数可以用
• f(xi+h)-f(xi): 节点xi的一阶向前差分 • f(xi)-f(xi-h): 节点xi的一阶向后差分 • f(xi+h)-f(xi-h): 节点xi的一阶中心差分
• 前后是相对x轴正方向而言 f(x x )f(x ) xd d xf(x ) ( 2 x !)2d d x 2 2f(x ) ( 3 x !)3d d x 3 3f(x ) ( 4 x !)4d d x 4 4f(x )
（9）
（9）式中的O(Δx2)项表示忽略掉这些项带来的误差中的最大项和Δx2成正比。
由（9）式得到导数的二阶精度（second order accurate）近似为：
df(x)f(xx)f(x-x)
dx
2x
（10）
（10）式称为中心差分形式（central-difference formula）。
11
d
( x )2d 2
( x )3d 3
( x )4d 4
f(x x )f(x ) x f(x )
d x
2 !
d x2f(x )3 !
d x 3f(x )4 !
d x4f(x )
（1） （2）
（1）式和（2）式相加，得到：
d d x 2 2f(x)=f(x x) (2 fx ( )x 2)f(x- x)+ O ( x2)
(13)
t x
u(x,0) g(x), xR
(14)
假定以上问题的解u(x,t)是充分光滑的，由Taylor级数展开有：
26
利用（15）和（17）式，得到：
u (x j,tn 1 ) u (x j,tn ) a u (x j 1 ,tn ) h u (x j,tn ) u t a u x n j O ( h )
d x
x
由（7）式得到导数的另一个一阶精度近似：
d f(x)f(x)f(x-x)
dx
x
（8）
（8）式称为向后差分形式（backward-difference formula）。
（6） （7）
10
f(x x )f(x ) xd d xf(x ) ( 2 x !)2d d x 2 2f(x ) ( 3 x !)3d d x 3 3f(x ) ( 4 x !)4d d x 4 4f(x )
(SCI, IF: 3.303)
3
2.1 有限差分法基础
对于存在复杂介质和几何、特殊边界条件和初始条件的实际地质问题， 一般不存在解析解，需要近似的数值求解方法。 有限差分方法是地球物理方法中最常见的一种。
有限差分方法（Finite Difference Method, FDM）是计算机数值模拟最 早采用的方法，至今仍被广泛使用。 有限差分方法的基本特点
代入第二、三类边界条件
23
（2） 第二类和第三类边界条件 图中O与R重合
图中V与R点重合
24
2.4 差分方程
对于具体地球物理问题的偏微分方程组，利用上述差分格式，可以给出偏 导数的微商近似，进一步得到差分方程组。
设f(P)是内部区域DI上定义的一个函数，设L(u)是一个微分算子，则以下表 示了未知量u(P)的偏微分方程：
局部节点离散化方程 总体节点离散化方程 总体节点离散化方程
19
边界条件的离散化的处理
（1）第一类边界条件
狄利克莱问题
G g(s)
（a） 直接转移法
在图中网格是按正方形分割， 步长为h。0点为靠近边界G的一个网格节点， 1和2为边界节点。我们取最靠近0点的 边界节点1上的函数值作为0点的函数值。 即取φ0≈φ1。
15
一阶偏导数中心差分的推导 为了寻求更精确的差分格式，我们引入两个 待定常数，由泰勒展开，构造如下关系式
回代（1）中，舍去高阶项
（1）
16
二阶偏导数差分的推导 回代（1）中，舍去高阶项
（1）
二阶偏导数差分公式
17
一个例子： 等步长
18
一个例子：
f=0时，变为泊松方程 f=q=0时，变为拉普拉斯方程
df(x)f(x x)f(x) O ( x)
d x
x
（3） （4）
9
df(x)f(x x)f(x)O ( x)
dx
x
（4）
式中的O(x)项表示忽略掉的所有项中的最大项的量级 是Δx，也就是说，忽略掉这些项带来的误差中的最大 项和Δx成正比。
由（4）给出导数的一阶精度（first order accurate）近似为：
• 有限差分法以变量离散取值后对应的函 数值来近似微分方程中独立变量的连续 取值。
• 我们放弃了微分方程中独立变量可以取 连续值的特征，而关注独立变量离散取 5
有限差分法的主要内容
1. 建立地球物理问题的离散有限差分模型 （1）如何根据问题的特点将定解区域做网格划分； （2）如何在所有网格节点上用有限差分格式对导数求近似， 对函数、初始条件和边界条件求近似； （3）如何把原方程离散化为代数方程组，即有限差分方程组。
用Taylor级数展开可以给出微商的近似形式。
dx变为Δx
对于连续函数f(x)，它在相邻点上的值f(x+Δx)和f(x- Δx)可以用Taylor 级数展开为
f(x x )f(x ) xd d xf(x ) ( 2 x !)2d d x 2 2f(x ) ( 3 x !)3d d x 3 3f(x ) ( 4 x !)4d d x 4 4f(x )
d
( x )2d 2
( x )3d 3
( x )4d 4
f(x x )f(x ) x f(x )
f(x )
f(x )
f(x )
d x
2 ! d x2
3 ! d x 3
4 ! d x4
（1） （2）
如果Δx很小，f(x)可微，则以上级数收敛。 次数越高，收敛级数的项的绝对值越小。
由（1）得到，
f(x x)f(x) xd d xf(x) ( 2 x !)2d d x 2 2f(x) ( 3 x !)3d d x 3 3f(x) ( 4 x !)4d d x 4 4f(x)
这种方法称为直接转移法，又称为零次插值法。
20
（b） 线性插值法 先判断x方向的边界节点1和y方向的边界节点 2哪一个更靠近0点。
如果1更靠近0点，则可以用x方向的线性插 值给出0点的函数值
如果2更靠近0点，则可以用x方向的线性插 值给出0点的函数值
21
（c） 双向插值法 变步长二次偏导数
j+1 (i-1,j+1) (i,j+1) 2