计算土力学
主讲教师:张爱军
? 其他差分公式介绍
?以上讲的差分公式为中心差分公式,
其特点是需要采用前一步和后一步值
来进行计算。
?另外还由向前和向后差分方法,这些
方法的实质是采用本身值和前一步或
后一步值来计算一次微分
?以下列出这些差分的公式
方法 差分公式 适用范围
中心差分 精度较高
向前差分 精度差一些
向后差分 精度差
1,1,
2
1,,1,
22
2
2
i j i j
i j i j i j
fff
xh
f f ff
xh
??
??
??
?
?
???
?
?
1,,
2
2,1,,
22
2
i j i j
i j i j i j
fff
xh
f f ff
xh
?
??
??
?
?
???
?
?
,1,
2
,1,2,
22
2
i j i j
i j i j i j
fff
xh
f f ff
xh
?
??
??
?
?
???
?
?
方法 差分公式 适用范围
半步长中心
差分
精度较高
时间域上差分
用之
0,5,0,5,
2
0,5,,0,5,
22
4 ( 2 )
i j i j
i j i j i j
fff
xh
f f ff
xh
??
??
??
?
?
???
?
?
? 差分公式的截断误差
?差分公式实际上是一种近似解法,其
误差分析是非常重要的
?截断误差值分析是分析在采用差分公
式代替微分公式时,其误差有多大。
截断误差分析时进行差分方程收敛性
和稳定性分析的基础。
?向前差分公式截断误差
由 Taylor级数展开式得到:
将以上公式移项得到:
2 2 3 3
1,,232 ! 3 !
i ii
i j i j
xx x x x x
f h f h ff f h
x x x? ? ??
? ? ?? ? ? ? ?
? ? ?
2 2 3
1,,
23
2 ! 3 !
()
i ii
i
i j i j
xx x x x x
xx
ff f f h f
h x x x
h
f
x
h
O
?
? ??
?
? ? ? ?
? ? ? ?
? ? ?
?
??
?
截断误差
?向后差分公式截断误差
由 Taylor级数展开式得到:
将以上公式移项得到:
2 2 3 3
1,,232 ! 3 !
i ii
i j i j
xx x x x x
f h f h ff f h
x x x? ? ??
? ? ?? ? ? ? ?
? ? ?
2 2 3
,1,
23
2 ! 3 !
()
i ii
i
i j i j
xx x x x x
xx
ff f f h f
h x x x
f
O
x
h
h
?
? ??
?
? ? ? ?
? ? ? ?
? ? ?
?
??
?
因此:从上式可以看出,向前与向后差分是
对于某个方向上的截断误差是其步长 h的一
次函数。我们说向前和向后差分的精度是 一
阶 近似的。
? 中心差分公式截断误差
2 2 3 3
1,,232 ! 3 !
i ii
i j i j
xx x x x x
f h f h ff f h
x x x? ? ??
? ? ?? ? ? ? ?
? ? ?
2 2 3 3
1,,232 ! 3 !
i ii
i j i j
xx x x x x
f h f h ff f h
x x x? ? ??
? ? ?? ? ? ? ?
? ? ?
3 4 5
1,1,
3
2
52 3 ! 5 !
i ii
i j i j
xx x x x x
ff f f h f
h x x
h
x
??
? ??
? ? ? ?
? ? ?
? ? ?
从上式可以看出,中心差分公式较向前和向后
差分公式的精度提高一个等级。
可以推导出二阶微分的精度如下,
1,21,()
2
i
i j i j
xx
ff f hO
hx
??
?
? ???
?
2 2 3 3
1,,232 ! 3 !
i ii
i j i j
xx x x x x
f h f h ff f h
x x x? ? ??
? ? ?? ? ? ? ?
? ? ?
2 2 3 3
1,,232 ! 3 !
i ii
i j i j
xx x x x x
f h f h ff f h
x x x? ? ??
? ? ?? ? ? ? ?
? ? ?
2 4 4 6
1,,1,
2 2 4 6
22
4 ! 6 !
i i i
i j i j i j
x x x x x x
f f f f f h f
h x x
h
x
??
? ? ?
?? ???? ? ?
? 从上式可以看出,对于二阶偏导,向前、
向后和中心差分公式的精度是一致的。
? 从以上公式看出,对于以上差分公式,
当 h无限减少时,截断误差可以无限减少,
2
1,,1,
22
2
2,1,,
22
2
,1,2,
22
2
2
2
2
()
2
()
2
()
i
i
i
i j i j i j
xx
i j i j i j
xx
i j i j i j
xx
f f f f
O
hx
f f f f
O
hx
f f f f
O
hx
h
h
h
??
?
??
?
??
?
?? ?
??
?
?? ?
??
?
?? ?
??
?
同 理,
?但是注意的是,对于一个微分方程而
言,有两个方向、或三、四个方向的
差值公式,相互有影响,一个方向截
断误差对于其他方向还有影响,需要
进一步分析。 —— 结合差分格式分析。
§ 3.2 差分格式
? 将差分公式代入基本控制方程得到的方
程成为差分方程,也称为差分格式。对
于同一个微分方程组和定解条件可以建
立不同的差分方程(差分公式不同、使
用形式不同,方程不同),构造差分格
式也有不同的途径。
? 以扩散方程为例说明差分格式的建立
其中,a 为 常量 参数
u 为要求的未知量,可以是温度,
水压等等扩散的物理量
x,t分别为空间和时间自变量
这个方程是物理量随着时间和空间(一
维)的扩散计算的方程,我们讨论这个
方程的差分解法。--代表性强。
2
2 0 ( 0)
uuaa
tt
??? ? ?
??
设:在 x- t平面上 作平行于 x轴和 t轴的两组平
行线,形成 h× τ的差分网格,则:
其中:时间的步长为 τ
空间的步长为 h,
且称网格比为:
规定:差分表示为:
0,1,2
0,1,2
x j h j
t n n?
? ? ? ?
??
2h
???
nju
上标为时间差分的结点号
下标为空间差分的结点号
表示在时间为 n,空间为 j点的变量值
t
x
? 若对于时间 t方向上采用一次向前差分公
式,而对于 x方向上采用中心差分公式,
得到扩散方程的差分格式:
这称为求解扩散方程的显式格式。
? 若变化 x,t差分格式的 形式 就变成为:
这也称为求解扩散方程的显式格式。
1
11
2
2 0n n n n nj j j j ju u u u ua
h?
?
??? ? ???
11
11
2
2 0
2
n n n n n
j j j j ju u u u ua
h?
??
??? ? ???
? 若变化 x,t差分格式的 关系 就变成为:
这也称为求解扩散方程的 隐式 格式。
? 若再变化 x,t差分格式的 关系 就变成为:
这就称为求解扩散方程的加权 隐式 格式。其
中 θ为小于 1的加权系数。当 θ= 0.5时称为
Crank-Nicolson公式,还有 Douglas公式等多种
形式。
注:这些公式名称只对于扩散方程而言
11
1
11
1
2
2 0nnj j j j jn n nu u u u ua
h?
???
??
?? ? ?
??
1
1 1 1
1
12
111 ( ) [ ( 2 ) ( 1 ) ( 2 ) ] 0n n n n nnn
j j j j j j j
n
j
au u u u u u u u
h ???
??
? ? ? ?
??? ? ? ? ? ? ? ? ?
? 相对于在 t上的加权差分格式:
称为 Richtmyer隐式差分格式
? 这就引出差分解需要解决的两个问题:
?如何构造一个方程的差分格式?
?那种差分格式好,那种能用,那种不
能用?
1 1 1 1 1
11
2
231( ) 0
22
n n n n n n n
j j j j j j ju u u u u u ua
h??
? ? ? ? ?
??? ? ? ?? ? ?
? 构造差分格式的方法
? 数值微分法
这是我们以上讲的方法,就是用适当的差商
代替微商,从而得到任意逼近微分方程的差
分格式,常用、直接、简单。
? 积分插值法:从守恒定律出发构造差分格式。
? 待定系数法:选用形式确定而系数待定的差
分方程逼近微分方程,然后在截断误差可能
达到的范围内,按精度要求定出差分系数,
构造具体的差分格式。
—— 积分积分差值法和待定系数法自己看。
? 用“收敛性”和“稳定性”判断差分格式的优

? 差分格式的“收敛性”
?,收敛性”就是:当步长无限减少时差分公
式是否能够无限接近原微分方程,其截断误
差是否趋于 0。
? 收敛性的数学定义为:
设 u为微分方程的准确解,unj为相应的差分
方程的解,如果步长,对于任意点
( j,n)有:
则:称差分是收敛的。
(,)nj j nu u x t?
0,0hl??
? 对于扩散方程的总体截断误差为:
设,u( xj,tn )为微分方程的解,ujn为差分
格式的解,令 E为差分格式的截断误差,记
结点( xj,tn)处的截断误差为 Ej,n,则:
1
11
2
1
11
2
0
( 1 2 ) ( )
n n n n n
j j j j j
n n n n
j j j j
u u u u u
a
h
u a u a u u
?
??
?
??
?
??
? ? ?
??
? ? ? ?
网格比
24 2
,24
2 42
24
2
(,) (,)
2 1 2
(,) (,)
2 1 2
()
j n j n
jn
j n
u x t u x th
Ea
tx
ux utah
tx
Oh
?
? ??
?
??
? ? ? ?
??
? ?
??
??
??
是 u,非 ujn
另外,根据截断误差的定义可知:
将上式变化成为:
令:
1 1 1
,2
2
2
1 1 1
2
(,) (,) (,) 2 (,) (,)
{ (,) (,) }
(,) (,) (,) 2 (,) (,)
j n j n j n j n j n
jn
j n j n
j n j n j n j n j n
u x t u x t u x t u x t u x t
Ea
h
u x t a u x t
tt
u x t u x t u x t u x t u x t
a
h
?
?
? ? ?
? ? ?
? ? ?
??
??
??
??
? ? ?
??
将真解代入差分
方程中
为 0
1 1 1,(,) ( 1 2 ) (,) [ (,) (,) ]j n j n j n j n j nu x t a u x t a u x t u x t E? ? ?? ? ?? ? ? ? ?
,(,)nj n j j ne u u x t?? 1,1 1(,)nj n j j ne u u x t?????
差分解 微分解
得到:
当 0< λ ≤1/2/a时:
,1,1,1,,( 1 2 ) ( )j n j n j n j n j ne a e a e e E? ? ?? ? ?? ? ? ? ?
1 11(1 2 ) ( )n n n nj j j ju a u a u u??? ??? ? ? ?
1 1 1,(,) ( 1 2 ) (,) [ (,) (,) ]j n j n j n j n j nu x t a u x t a u x t u x t E? ? ?? ? ?? ? ? ? ?
,1,1,1,,
,,
( 1 2 ) ( )j n j n j n j n j n
j n j n
e a e a e e E
eE
? ? ?
?
? ? ?? ? ? ? ?
??
均大于 0
得到:
得到:当 时,
因此:证明差分格式当 u光滑时是收敛的。
同样可以证明:不满足 0< λ ≤1/2/a时差分是不
收敛的。
对于其他差分格式,也可以得到其收敛条件,有
条件收敛的差分称为 条件收敛格式,无条件收
敛的差分格式为 无条件收敛格式 。
2,1,0,,()j n j j n j ne e n E n E n O h? ? ? ?? ? ? ? ? ? ?
,0 (,)nj n j j ne u u x t??即,
0,0hl??
? 差分格式的“稳定性”
差分计算是一个逐次计算的过程,前一
次计算的误差会向下一次计算传播,这
种误差是否越来越大,在次数达到无限
大时,误差不可控制称为“不稳定的”,
可以控制的称为“稳定的”。
数学定义:设初始层引入误差为
J=0,1,2,.,而第 n层的引入误差为 那么:
称为差分为稳定的,其存在常数 k,使得:
成立,其中,为误差的
范数 —— 度量某种向量、矩阵的尺度。
0j?
0njjk????? m axnnjj??? ?
?稳定性的证明较为复杂,不讲参考文
献:陆金莆,顾丽珍、陈景良,偏微
分方程的差分方法,高等教育出版社,
1988。
?其方法有:傅立叶方法、能量方法、
单调矩阵方法和离散格林( Green )
方法
? 几种差分方法的收敛性和稳定性
格式 差分公式 截断误差 稳定条件
显式格

隐式格
式 无
加权隐
式格式 无
Crank-
Nicolson 上式中若,无
Richtmy
er隐式
格式

1 1121 ( ) ( 2 ) 0n n n n nj j j j jau u u u uh? ? ??? ? ? ? ?2()Oh? ? 1
2a??
1 1 1 11121 ( ) ( 2 ) 0n n n n nj j j j jau u u u uh? ? ? ? ???? ? ? ? ?2()Oh? ?
1
112
11
1 ( ) [ ( 2 )
(1 ) ( 2 ) ] 0
n n n n n
j j j j j
n n n
j j j
au u u u u
h
u u u
??
?
?
??
??
? ? ? ?
? ? ? ? ?
12?? 2()Oh? ?
1 1 1 1 111
2
231( ) 0
22
n n n n n n nj j j j j j ju u u u u u ua
h??
? ? ? ? ???? ? ? ?? ? ?
§ 3.3 差分方法的解题步骤
? 解题步骤
?对求解域进行网格划分
?选择差分格式
?边界结点和内结点采用不同的差分格

?建立差分方程
?求解差分方程
? 举例:软粘土地基非线性一维固结分析
qu
O
H
透水
kv,mv,cv=kv/γw/mv
Z
'
0
0 '
0
0
lo g ( )
lo g ( )
v
v
k
v
e e c
k
e e c
k
?
?
??
??
非线性特征
固结系数
e
σ’( 对数坐标 )
kv
e( 对数坐标 )
得到:
其中,cc 为压缩指数,ck为渗透指数( e~
logkv的斜率),cv为固结系数,下标 0表
示初始状态。
在瞬时均布荷载作用下,不
不随着时间的变化而变化,总应力不变,
则:
'
0
0 '
'
0
0 '
()
()
c
k
c
c
vv
vv
kk
mm
?
?
?
?
?
?
'' 0 uu q c on t? ? ?? ? ? ? ?
' u
tt
?????
??
又:一维饱和土体的连续方程为:
得到:
其中,V 表示与固结度有关的函数
Tv表示时间因素,c为常数
现在看如何用差分法进行计算:
0
11()
1v w
uek
z z e t?
? ? ??
? ? ? ?
2
21
2 ()
c
v
V c V VV
z V z T
?? ? ???
? ? ?
运用 Crank-Nicolson差分格式离散 Z方向上
的微分,用向前差分离散时间方向上的
微分,得到:
化解得到:
1 1 1 1 1 1
1 1 1 1 1 1
21
1
11
22
()
2 ( ) 2
()
j j j j j j j j
i i i i i i i i
j
i
jj
jc ii
i
V V V V V V V Vc
z V z
VV
V
T
? ? ? ? ? ?
? ? ? ? ? ?
?
?
??
? ? ? ? ? ?
?
??
?
?
?
1
11
1 1 1 1 1 1 1 2
1 1 1 11
1
2 [ 1 ( ) ]
1
2 [ 1 ( ) ] ( )
2
j j c j j
i i i i
j j c j j j j
i i i i i ij
i
V V V V
c
V V V V V V
V
?
?
?
??
? ? ? ? ? ? ?
? ? ? ??
? ? ?
? ? ? ? ? ? ?
式中,ΔZ,ΔT为步长,λ= ΔT/ (ΔZ)2
i为空间结点序号,j为时间结点序号
上式实际上是一个递推公式,又由初始条
件,结合边界条件
得到:
0 1iV ? 0 1 1,( )j j jnnV b V V???? 单 面 排 水
1 1 1
2 2 2
3 3 3
1 1 1
1
11
11
.,,,,
11
0
n n n
n n n
A V B
A V B
A V B
A V B
A V B
? ? ?
? ? ? ? ? ?
? ? ? ? ? ?
? ? ? ? ? ?
? ? ? ? ? ???
?? ? ? ? ??
? ? ? ? ??
? ? ? ? ??
? ? ? ? ??
? ? ? ? ??? ? ? ? ? ?
由上线性方程就可以求出各点的 V值,也
就可以求出各点的固结度值。
经过验证,该法与解析解符合的很好。
Any questions?