显式差分和隐式差分

一种隐式差分格式的程序实现
1
2
3
4
5
6
求解区域：
边界条件：
初始条件：
内部节点：
A = sparse(nx,nx); for i=2:nx-1 A(i,i-1) = -s; A(i,i ) = (1+2*s); A(i,i+1) = -s; end
边界节点：
A(1 ,1 ) = 1; A(nx,nx) = 1;
相容性是比较容易满足的。在此基础上，如果满足了稳定性条件，差分格式的 收敛性就自动满足。
2.5 有限差分法实例
U=100
U=0
2U 0
j+1 (i-1,j+1) (i,j+1) 2
j (i-1,j)
j-1 (i-1,j-1)
3 h2 0 (hi3,j) h1 h44 (i,j)
(i+1,j+1)
T in 1 tT in 2 ( x )2(( T i n 1 1 2 T in 1 T i n 1 1 ) ( T i n 1 2 T in T i n 1 ))
s T i n 1 1 ( 2 2 s ) T i n 1 s T i n 1 1 s T i n 1 ( 2 2 s ) T i n s T i n 1
内部节点：
for j=2:n-1 for i=2:m-1;
a((j-1)*m+i,(j-1)*m+i+1)=1; a((j-1)*m+i,(j-1)*m+i-1)=1; a((j-1)*m+i,j*m+i)=1; a((j-1)*m+i,(j-2)*m+i)=1; a((j-1)*m+i,(j-1)*m+i)=-4; end end
Tj1Tj tf(Tj,tj)
Tj1Tj- tTj=T（ j 1- t）
Tj
T（0 1-
t
）j
现在改用向后差分形式进行近似，得到隐式差分格式：
显式差分格式
dT dt ttj
Tj Tj1O(t) t
Tj Tj1tf(Tj,tj)
f (T ,t) T
Tj Tj1(1t)1
隐式差分格式
Tj
T0(1
t
边界节点：
A矩阵非零系数减少， 同时引入第一类边界，
方程右端项B向量出现 非零元素。
AX B
组建A和B矩阵，求解线性方程组得到X
AA(135,135)
XX(135,1)
BB(135,1)
%Matlab 2D clear; clc; figure('color','w');
a=zeros(135,135); for i=1:135
显式与隐式差分格式
主讲人：胡才博 中国科学院大学地球科学学院 中国科学院计算地球动力学重点实验
室
显式和隐式：求解问题与时间相关
• 显式差分格式（explicit difference scheme） 差分方法中可逐层逐点分别求解的格式。
• 特点
Y(tt)F(Y(t))
• 1. 不联立解方程；
• 2.时间步长和空间步长的选择受限制。通常要求时间步长 足够小。
)j
可以验证，当时间步长趋近于零时，以上近似解趋于解析解。因此，
该格式收敛。 稳定性条件：T单调减小的条件
0
1 1 t
1
Tj Tj1(1t)1
0
1 1 t
1
重点考察差分格式的收敛性
隐式差分格式 无条件稳定
当dt=1.25，tau=0.7时，显式差分格式不稳定，结果振荡； 隐式差分格式稳定，结果不精确。
Ac(n1) Bc(n)
A = sparse(nx,nx);
1 0
0
0
0 0
s
2 2s
s
0
0
0
0 s 2 2s s 0 0
A
0
0
s 2 2s s
0
0 0 0 0
0 s 2 2s s
0
0
0 1
内部节点：
for i=2:nx-1 A(i,i-1) = -s; A(i,i ) = (2+2*s); A(i,i+1) = -s; end
1 1
1 0
0
0
0 0
s
2 2s
s
0
0
0
0 s 2 2s s
0 0
B
0
0
s
2 2s
s
0
0 0
0
s 2 2s s
0 0
0
0
0 1
T T
n 1
n 2
c (n )
T T
n 3
n 4
T T
n 5
n 6
s T i n 1 1 ( 2 2 s ) T i n 1 s T i n 1 1 s T i n 1 ( 2 2 s ) T i n s T i n 1
显示差分格式示意图
2. 隐式差分格式：
时间一阶精度 空间二阶精度
隐式有限差分格式
Crank-Nicolson 隐式差分格式
T in 1 tT in 2 ( x )2(( T i n 1 1 2 T in 1 T i n 1 1 ) ( T i n 1 2 T in T i n 1 ))
隐式差分格式（implicit difference scheme）
特点
G (Y(t),Y(t t))0
1. 时间步长和空间步长的选择不受限制；
2. 需要联立解方程组
例子： 1. 显式差分格式：
左端：n+1时刻的值； 右端：n时刻的值。
特点：结构简洁，直接求解，求解速度快。
但是，时间步长需满足：
显式差分格式才能得到稳定的数值解，否则，数值解将会不稳定而振荡。
1 (i+1,j) (i+1,j-1)
i-1
i i+1
U=0
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 U=0
U ( i 1 , j ) U ( i 1 , j ) U ( i , j 1 ) U ( i , j 1 ) 4 U ( i , j ) 0
j+1 (i-1,j+1) (i,j+1) 2
Tair
Tcap
一阶常微分方程的数值解 首先对时间和温度进行离散：
dT f (T , t) dt
利用向前差分形式：
tj t0jt,Tj T(tj)
dT dt ttj
Tj1Tj O(t) t
得到以下的显式差分格式：
Tj1Tj tf(Tj,tj)
利用向前差分格式：
dT dt ttj
Tj1Tj O(t) t
Cui, 2013
Cui, 2013
Cui, 2013
Cui, 2013
总结：
1、有限差分方法给出的数值解的精度取决于所用的差分形式（向 前、向后、中心）。
2、偏微分方程的显式有限差分格式通常是有条件稳定的，为了保 证得到精确的数值解，最关键的是需要根据稳定性条件选取正确的 空间和时间步长。
解，称该差分格式是收敛的。
如果
h l0i,m t 0(Uim ,j,kuim ,j,k)0
则称该差分格式是收敛的。
收敛性描述的是当差分网格无限细化时，差分方程的解是否具有无限逼 近偏微分方程的解的能力
Lax等价定理（Lax equivalence theorem）：如果逼近一个给定问题的差 分格式是相容的，那么该差分格式的收敛性与稳定性互为充分必要条件。
a(121,106)=-0.25; a(135,134)=-0.25; a(135,120)=-0.25; a(15,14)=-0.25; a(15,30)=-0.25; for i=2:14
a(i,i-1)=-0.25; a(i,i+1)=-0.25; a(i,i+15)=-0.25; end for i=122:134 a(i,i-1)=-0.25; a(i,i+1)=-0.25; a(i,i-15)=-0.25; end for i=1:7 for j=2:14; a(15*i+j,15*i+j-1)=-0.25; a(15*i+j,15*i+j+1)=-0.25; a(15*i+j,15*i+j+15)=-0.25; a(15*i+j,15*i+j-15)=-0.25; end end
1
2
3
4
5
6
边界条件：
初始条件：
Ac(n1) Bc(n)
1 0
0
0
0
s
2 2s
s
0
0
0 s 2 2s s 0
A
0
0
s 2 2s s
0
0
0
0
c (n1)
T T
n 1
n 2
1 1
T T
n 3
n 4
1 1
0 0 0 0
0 s 2 2s s
0
0
0 1
T T
n 5
n 6
B(1 ,1 ) = 1; B(nx,nx) = 1;
T T
n 1
n 2
c (n )
T T
n 3
n 4
T T
n 5
n 6
例子：牛顿冷却定律：温度高于周围环境的物体向周围媒质传递热量 逐渐冷却时所遵循的规律。当物体表面与周围存在温度差时，单位时 间从单位面积散失的热量与温度差成正比。
subplot(1,2,1),mesh(s) axis([0,17,0,11,0,100]) subplot(1,2,2),contour(s,32)
11
10
100
9
80
8
60
7
40
6
20
5
4 0
10
3
15
5
2 10
5 00
1
5
10
15
2.5 应用实例
南加州一次未来大地震的强地面运动的数值模拟
盆地效应
回顾
1. 有限差分法基础 2. 差分格式 3. 差分方程 4. 边界条件的处理 5. 相容性、稳定性和收敛性
回顾
1. 有限差分法的相容性、稳定性和收敛性
相容性：针对差分格式而言，在时间步长和空间步长趋近于零的情况下， 如果差分格式的截断误差（差分格式与原有偏微分方程之差）的模趋近于零， 则该差分格式与原偏微分方程是相容的，或称该差分方程与原偏微分方程 具有相容性。
a(i,i)=1; end; for i=1:7
a(15*i+1,15*i+2)=-0.25; a(15*i+1,15*i+16)=-0.25; a(15*i+1,15*i-14)=-0.25; end for i=1:7 a(15*i+15,15*i+14)=-0.25; a(15*i+15,15*i+30)=-0.25; a(15*i+15,15*i)=-0.25; End a(1,2)=-0.25; a(1,16)=-0.25; a(121,122)=-0.25;
边界节点：
A(1 ,1 ) = 1; A(nx,nx) = 1;
1 0
0
0
0 0
sBiblioteka Baidu
2 2s
s
0
0
0
0 s 2 2s s
0 0
B
0
0
s
2 2s
s
0
0 0
0
s 2 2s s
0 0
0
0
0 1
内部节点： 边界节点：
B= sparse(nx,nx); for i=2:nx-1 B(i,i-1) = s; B(i,i ) = (2-2*s); B(i,i+1) = s; end
稳定性（stability）：如果偏微分方程的严格解析解有界，差分格式给出的 解也有界，称该差分格式是稳定的；如果差分格式给出的解是无界的，则 称该差分格式是不稳定的。
稳定性反映了差分格式在计算中控制误差传递的能力
收敛性（convergence）：如果当时间和空间步长趋于零时，FDE解趋于PDE
载荷项：
rhs = zeros(nx,1); rhs(2:nx-1) = Told(2:nx-1); rhs(1) = Tleft; rhs(nx) = Tright;
内部 边界
Crank-Nicolson 隐式差分格式的程序实现
s T i n 1 1 ( 2 2 s ) T i n 1 s T i n 1 1 s T i n 1 ( 2 2 s ) T i n s T i n 1
Crank-Nicolson 隐式差分格式
s T i n 1 1 ( 2 2 s ) T i n 1 s T i n 1 1 s T i n 1 ( 2 2 s ) T i n s T i n 1
一般差分格式
1/2
Forward-Time Central-Space method Backward -Time Central -Space method Crank-Nicolson 隐式差分格式
b=a^(-1); c=zeros(135,1); for i=121:135
c(i,1)=25;end d=b*c; s=zeros(11,17); for i=2:16
s(11,i)=100; end for i=1:9 for j=1:15; s(i+1,j+1)=d(15*(i-1)+j,1); end end
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
U ( i 1 , j ) U ( i 1 , j ) U ( i , j 1 ) U ( i , j 1 ) 4 U ( i , j ) 0