第二章 有限差分法
合集下载
有限差分法
两端都要给定边界条件（双程坐标） 。
9
(C) 双曲型方程：适当的边界条件和初始条件，与波动传 播的性质有关 如：一维对流方程
∂u ∂u +c =0 ∂t ∂x u (x ,0) = f (x )
解为 u (x , t ) = f (x − ct ) ，代表一个向右(c > 0 时)或向左 ( c < 0 时)传播的波形。必须在波形传来的一侧提供边界条 件（单程坐标） 。
10
不适定的例子：
utt + u xx = 0 u (x ,0) = u t (x ,0) = 0
拉普拉斯方程＋非闭域边界条件，解为 u (x , t ) ≡ 0 。 然而，若定解条件为 u (x ,0) = 0, ut (x ,0) =
u (x , t ) = 1 sin nx ，解为 n
1 sinh nt sin nx n
(
)
n n um+1 = um −
cτ n n um +1 − um −1 2h
(
)
设计算到第 n 步时的累积误差
n ~n εn = 计算值um − 差分法精确解um m
反之
n ~n um = εn + um m
15
则第 n+1 步的计算值
~n ~ n cτ u n − u n ~ ~ um+1 = um − m +1 m −1 2h cτ n cτ n n n = um − um +1 − um −1 + εn − εm +1 − εn −1 m m 2h 2h n = um+1 + εn +1 m
uin +1 − uin −1 uin+1 − uin +1 − uin −1 − uin−1 −α =0 Lh u = τ h2 ατ 2 ⎛ ∂ 2u ⎞ τ 2 ⎛ ∂ 3u ⎞ Ti = Lh u − Lu (x i , t n ) = 2 ⎜ 2 ⎟ + ⎜ 3 ⎟ − L 截断误差 6 ⎜ ∂t ⎟i h ⎜ ∂t ⎟i ⎝ ⎠ ⎝ ⎠
2 有限差分法
K e L (ba ) 1 L (ba ) Dh e | y (a ) y0 |, L 0 | yn 1 ( xn 1 ) yn 1 | L K | y ( a ) y | D ( b a ) h , L0 0
稳定性
在具体计算时步骤中出现的误差对结果的影响
前向差分
例题 求
2 2 [ 2 ( x, y ) 2 ( x, y )] f ( x, y ) x x
的中心差分逼近
解：令 ij ( xi , y j ) i 1, j ( xi h, y j ) i , j 1 ( xi , y j h)
四、工程应用
有限差分法求解问题的流程图
1 | iN iN ,j ,j |
| iN ,j |
1 iN ,j
四、工程应用
举例
3 2 1 4 5 6 9 8 7
4×4个网格， 9个内结点， 16个边界点
内点差分方程 点1： 点2： …… 点9：
6 2 12 26 41 h 0 5 3 13 1 42 h 2 0 20 18 4 8 49 h 2 0
1
2
0 3
4 (1 2 3 30 ) 2 3h 0
三、差分方程组求解
直接求解线性方程组
P.32
结点N<500
迭代法求解线性方程组 高斯-赛德尔迭代法求解线性方程组 超松弛迭代法求解线性方程组
return
四、工程应用
有限差分法求解问题的步骤
P.34
[ i 1, j 2ij i 1, j x
2
得：
第二章 有限差分法初步-1
2 2 1 1 T T T i , j i 1 , j 2 2 2 2 i , j 1 (y ) (x) (x) (x)
q T T （2.22） i , j 1 i , j 1 2 2 k (y) (y) 1 1
若x y ， q 0，则式（2.22）又被
O(x)
由式（2.8）可得
（2.9）
T ( x) T ( x x) x T ( x) T ( x) x 2!
O(x)
（2.10）
（2.9）+（2.10），得到
T ( x x) T ( x x) (x) T ( x) T ( x) 2x 3!
而用中心差商代替微商，其截断误差是与 同量级的小量；
(x)
2
中心差商的截断误差小于向右差商或向左差商。
（V）二阶差商
上述一阶差商一般仍是x的函数，对它们还可以 求差商。这种一阶差商的差商称为二阶差商， 它是二阶微商的近似，常用向右差商的向左差 商来近似二阶微商，即：
T ( x x) T ( x) T ( x) T ( x x) 2 d T x x 2 x dx
(x) 2 T ( x x) T ( x) xT ( x) T ( x) 2! 3 (x) 4 T ( x) O(x) （2.8） 3!
由式（2.7）可得
T ( x x) T ( x) x T ( x) T ( x) x 2!
2 2
（2.13）
对区域内各个点都成立的，当然对任意一个内节 点（i，j）也成立。
2T 在（i，j）处存在二阶偏微商 x 2 与 i , j
这些二阶偏微商所对应的差商可表示成：
椭圆型方程
(1.5)
注 此方程组尽管是高阶方程组，但每个方程未知数
最多有3个易于求解.
④ 对方程组 (1.4)～(1.5) 的解分析需要考虑以下几个问题：
(a) 解是否惟一？ (b) 当网格无限加密时，即 h 0 时，差分解 ui
是否收敛到真解 u (xi ) ? (c) 在何种度量下收敛？ (d) 收敛速度如何？ 为了解决如上问题，需要给出如下说明：
于是在 xi 将方程 (1.1) 写成
u (xi1) 2u (xi ) u (xi1) h2
q(xi )
u (xi )
f
(xi )
R
i(u),
(1.3)
其中
R
i(u)
h2 12
d
4u(x) dx4
i
O(h3 ).
舍去 R i(u) 得逼近方程 (1.1) 的差分方程为：
du dx
i
hi1 2
hi
d 2u dx2
i
O(h2
)
(2.3)
p(
x i
1
)
2
u(xi ) u(xi1) hi
p
du dx i1
2
hi2 24
p
d 3u
dx3
i1
2
O(h3)
p
du dx
取 x(1) x0 a, x(2) x1 , 得
2
(2.7) (2.8)
(2.9) (2.10)
W (a) W (x1 ) 2
x1
第二章 有限差分法的基本概念
x
间距h > 0称为空间步长，间距τ > 0称为时间步长．
2 用Taylor级数展开方法建立差分格式
设 f ( x) 在 x0 的某个邻域 U ( x0 , δ ) 内具有直 到n + 1阶的导数，则 ∀x ∈ U ( x0 , δ ) 有 f ( x) = f ( x0 ) + f ′( x0 )( x − x0 ) + + f ( n ) ( x0 ) ( x − x0 ) n + Rn ( x) n! Rn ( x)是余项，且Rn ( x) = o(( x − x0 ) n )( x → x0 )．
∫
xj +h 2
xj −h 2
[u (tn + τ , x) − u (tn , x)]dx
tn+1
∂u ∂u h = a ∫ [ (t , x j + 2 ) − (t , x j − h )]dt 2 tn ∂x ∂x 应用数值积分可得： [u (tn + τ , x j ) − u (tn , x j )]h ∂u ∂u h ≈ a[ (tn , x j + 2 ) − (tn , x j − h τ 2 )] ∂x ∂x
∂t = O(τ + h)
−(
∂u ( x j , tn )
我们也用“精度”一词说明截断误差． 一般，如果一个差分格式的截断误差Ｔ = O(τ q + h p )， 就说差分格式对时间t (τ )是q阶精度的， 对空间x(h)是p阶精度的． 特别，当p = q时，说差分格式是p阶精度的． 差分格式(1.13),(1.15),(1.17)都是对t (τ )一阶精度， 对x(h)二阶精度．而差分格式(1.11)是一阶精度格式．
第二章椭圆型方程的有限差分法
.
差分方程(1.6)当i 1,2, N 1,时成立，加上边值条件 就得到关于的线性代方数程组：
Lhui
ui1
2ui h2
ui1
qiui
fi ,i
1,2,
N 1,(1.8)
u0 , uN . (1.9)
它的解ui是u(x)于x xi的近似。称(1.8),(1.9)为逼近(1.1) (1.2)的差分方程或差分格。式
立 差 分 方 程 的 稳 定检性验。相 容 条 件 并 不。困我难们 曾
用Taylo展 r 式证明它都满足条相件容，并且估计了截
误 差 的 阶 。 因 此 我主们要的任 务 去 建 立 差式分的格稳
定 性 ， 即 建 立 形 (1.1如7)的 估 计 式 ， 称 之 为差关分于方
程解的先验估计。 .
的解u，由Taylo展 r 式可得
u(xi1)2u(xi )u(xi1) h2
d2u(x) [ dx2 ]i
1h22[h2dux(2x)]o(h3),(1.3)
其中[ ]i表示括号内函xi点 数取值。 于 是 在 可 (1.1)写 将成 方 程
u(xi1)2uh(2xi)u(xi1)q(xi)u(xi)f(xi)Ri(u)(, 其 中 Ri(u)1 h22 [h2du(2 xx)]o(h3), (1.5)
)
u(
xi1
)
q(
xi
)u(
xi
)
f (xi ) Ri (u) fi Ri (u)
与Lhui
ui1
2ui h2
ui1
qiui
fi
相减，得 Lh(u(xi ) ui ). Ri (u)
引进误差
ei u( xi ) ui , 则误差函数 eh( xi ) ei满足下列差分方程；
3第二章-有限差分方法基础
2.1.1 基本方程和定解问题
u t
2u x2
( 0)
求解域： (x, t) [0,1][0, ]
(2.1.1)
初始条件： u(x, 0) f (x)
边界条件： u(0, t) a(t), u(1, t) b(t)
(2.1.2)
方程(2.1.1)和初边条件(2.1.2)构成了一个适定的定解问题。
根据数学分析中的知识，我们知道
2u (x,t) lim u(x x,t) 2u(x,t) u(x x,t)
x2
x0
x2
所以，二阶导数可以近似为
2u
x
2
n
k
un k 1
2ukn x2
ukn
un k 1
2ukn
un k 1
称为二阶中心差分。
容易证明：
un k 1
2ukn
un k 1
t
)
ut
(
x,
t
)
lim
t 0
u(
x,
t
t
)u 2t
(
x,
t
t
)
其中，lim 后面的项称为差商(difference quotient)。 t 0
当t足够小时，可以用差商来近似导数。
即：
u(x,t t) u(x,t)
ut (x,t)
t
u(x,t) u(x,t t)
ut (x,t)
t
u(x,t t) u(x,t t)
The Elements of Computational Fluid Dynamics
第二章 有限差分方法基础
§2.1 有限差分方法概述 §2.2 导数的数值逼近方法 §2.3 差分格式的性质 §2.4 发展方程的稳定性分析
有限差分法
有限差分法有限差分法finite difference method微分方程和积分微分方程数值解的方法。
基本思想是把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。
然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。
有限差分法的主要内容包括：如何根据问题的特点将定解区域作网格剖分；如何把原微分方程离散化为差分方程组以及如何解此代数方程组。
此外为了保证计算过程的可行和计算结果的正确，还需从理论上分析差分方程组的性态，包括解的唯一性、存在性和差分格式的相容性、收敛性和稳定性。
对于一个微分方程建立的各种差分格式，为了有实用意义，一个基本要求是它们能够任意逼近微分方程，这就是相容性要求。
另外，一个差分格式是否有用，最终要看差分方程的精确解能否任意逼近微分方程的解，这就是收敛性的概念。
此外，还有一个重要的概念必须考虑，即差分格式的稳定性。
因为差分格式的计算过程是逐层推进的，在计算第n＋1层的近似值时要用到第n层的近似值，直到与初始值有关。
前面各层若有舍入误差，必然影响到后面各层的值，如果误差的影响越来越大，以致差分格式的精确解的面貌完全被掩盖，这种格式是不稳定的，相反如果误差的传播是可以控制的，就认为格式是稳定的。
只有在这种情形，差分格式在实际计算中的近似解才可能任意逼近差分方程的精确解。
关于差分格式的构造一般有以下3种方法。
最常用的方法是数值微分法，比如用差商代替微商等。
另一方法叫积分插值法，因为在实际问题中得出的微分方程常常反映物理上的某种守恒原理，一般可以通过积分形式来表示。
此外还可以用待定系数法构造一些精度较高的差分格式。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
(x
xi )3
1 4 f 4! x4
(x
xi )4
在式 2.2.1 中，分别取 x xi h, x xi h ,
得
2.2.1
f h2 2 f h3 3 f h4 4 f
fi1, j
fi, j
h( ) x
2
x2
6
x3 24 x4
fi1, j
fi,
j
h( f x
)
h2 2
2 f x2
h3 6
x
f y
fi1, j1
fi1, j1 fi1, j1 4hl
fi1, j1
2.2.14
4 f x 2 y 2
2 x2
2 f
y 2
1 h2l 2
( fi1, j1 2 fi1, j
fi1, j1 2 fi, j1 4 fi, j
2 fi, j1 fi1, j1 2 fi1, j fi1, j1 )
3 f y3
fi2, j 2 fi1, j 2 fi1, j 2l 3
fi2, j
2.2.12
4 f y 4
fi2, j 4 fi1, j 6 fi, j 4 fi1, j l4
fi2, j
另外，利用式样2.2.6 至 2.2.9， 可以导出混合导数的中心差分公 式
2.2.13
2 f xy
位于V内的结点称为内结点，位于V和S以外的结点称为外部结点，如图2-1所示。
除特殊情况外，我们一般只考虑位于V内和S上（简记为V+S）的结点，即内部结
点和边界结点。结点 xi , yi 简记为i, j
，函数f在此点的值简记为fi, j 。
y
外部虚拟结点
内部结点
边界结点
xi , yi
S
x0, y0
V
l
h
x
0
图2-1 差分网格的划分
导数的差分公式可从函数的Taylor级数展开式导出。以二元函数f x, y 为
例，在点 xi , yi 附近，函数 f x, y 沿x 方向可以展为Taylor级数如下：
f
x,
y
fi, j
f x
(x
xi )
1 2 f 2! x2
(x
xi )2
1 3 f 3! x3
3 f x3
2 f x2
f x
fi2, j 2 fi1, j 2 fi1, j 2h3
fi2, j
4 f x4
2 f x2
2 f
x2
fi2, j 4 fi1, j 6 fi, j 4 fi1, j h4
fi2, j
2.2.10 2.2.11
同样，利用式2.2.8 和2.2.9， 可得
3 f x3
h4 4 f 24 x4
2.2.2 2.2.3
假定h是充分小的，因而可以不计它的三次幂及更高次幂 的各项，则 式2.2.2及式2.2.3 简化为：
fi1, j
fi, j
h(f ) x
h2 2
2 f x2
fi1, j
fi, j
h(f ) h2 x 2
2 f x2
联立求解 f 及 2 f ，得差分公式 ： x x2
2.2.15
例如，我们可以在式 2.2.2 中把 h2 的项也略去不计，得出
f fi1, j fi, j
x
h
2.2.16
h 或者把式2.2.3 中的 2 项也略去不计，则得
f fi, j fi1, j
x
h
2.2.17
它们分别称为一阶导数的向前、向后差分公式。以这两个公式化为基础， 可以导出高阶导数向前向后的差分公式。这种差分公式虽然比较简单，但 除了对时间进行差分以外，很少采用。因为它们不具有对称性，应用时容
i, j为整数
其中x0 , y0 为x-y平面上任意一点，它可在区域V内，亦可在V外。h 为x方向
步长，l 为y方向步长。这两组平行线的交点称为结点。这两组平行线在V内组成的
网格称为差分网格，若l h称为矩形网格，l=h 则称为正方形网格。沿x和y方向
距离均不超过一个步长的两结点称为相邻结点。位于边界S上的结点称为边界结点，
2.2.4 2.2.5
f fi1, j fi1, j
x
2h
2 f x2
fi1, j 2 fi, j h2
fi1, j
2.2.6 2.2.7
在y方向，同理可得
f fi1, j fi1, j 2 fi, j fi1, j
y 2
l2
2.2.8 2.2.9
公式2.2.6 至2.2.9是基本的中心差分公式，可以从它们导出其他的中心 差分公式。例如，利用式2.2.6 和 2.2.7 可得
第二章 有限差分法
概述
差分法的基本思想是用差分网格离散求解域，用差分公式将科学问题的控制 方程（常微分方程或偏微分方程）转化为差分方程，然后结合初始及边界条件， 求解线性代数方程组。由于这种方法比较直观，容易编制程序，所以从20世纪40 年代以来，至今仍被广泛应用。
有限差分法最早是被用于工程科学。40年代后期，土工中的渗流及固结问题 开始采用差分法成功地解决了某些实际问题，如如土坝渗流及浸润线的求法，土 坝及地基的固结等；20世纪50年代及60年代初，弹性地基上的梁与板以及板桩也 用差分法来求解。60年代以后，由于有限元解法的灵活性以及边界元法的异军突 起，使差分法在土工中的应用暂时趋于停滞。然而在最近几年，差分法又有了新 进展。任意网格的差分，使这老方法又可以与有限元相匹敌。另外，在某些特定 的条件下，有限元法与差分法及边界元相结合来处理一个课题，比它们各自求解 更显出优越性。
总之，由于有限差分法的使用以及它与其它方法的相互配合，使求解 土工及其他工程学科课题又达到了一个新阶段，使数值方法解决问题的能 力提高到新的水平。
在应用有限差分法分析土工问题时需要重视下述三个问题：
1. 如何选取差分格式将控制微分方程离散成差分方程； 2. 如何保证差分方程的稳定性和收敛性； 3. 如何求解差分方程组。
易发生差错，而且h把2 l 2或 的项略去不计，精确度也较差。
又例如，我们还可以在式2.2.1中，分别取 x xi 2h, x xi 2h 得
fi2, j
fi,
j
2h
f x
2h2
2 f x2
前两个问题将在2.2和2.3节中介绍，第三个问题结合解题步骤在2.4 节介绍。最后通过一个实例分析介绍有限差分法的具体应用。
差分公式推导
建立差分公式前先要将求解域划分差分网格。以图2-1所示二维问题为例，在xy平面上作分别平行于x轴和y轴的两组平行线：
xi x0 ih y j y0 jl