第五讲——显式差分和隐式差分(5)(格式整齐)

高级材料
16
显式与隐式差分格式
主讲人：胡才博 中国科学院大学地球科学学院 中国科学院计算地球动力学重点实验室
高级材料
17
高级材料
10
盆地效应
高级材料
11
高级材料
Cui, 2013
12
高级材料
Cui, 2013
13
高级材料
Cui, 2013
14
高级材料
Cui, 2013
15
总结：
1、有限差分方法给出的数值解的精度取决于所用的差分形式（向 前、向后、中心）。
2、偏微分方程的显式有限差分格式通常是有条件稳定的，为了保 证得到精确的数值解，最关键的是需要根据稳定性条件选取正确的 空间和时间步长。
解，称该差分格式是收敛的。
如果
lim (U
h0,t 0
m i, j
,k

um i, j
,k
)

0
则称该差分格式是收敛的。
收敛性描述的是当差分网格无限细化时，差分方程的解是否具有无限逼 近偏微分方程的解的能力
Lax等价定理（Lax equivalence theorem）：如果逼近一个给定问题的差 分格式是相容的，那么该差分格式的收敛性与稳定性互为充分必要条件。
j (i-1,j)
j-1 (i-1,j-1)
3
h2 0
(hi3,j)
ห้องสมุดไป่ตู้h1
h4 4
(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) 4U (i, j) 0
end
end
高级材料
边界节点：
A矩阵非零系数减少， 同时引入第一类边界，
方程右端项B向量出现 非零元素。
6
AX B
组建A和B矩阵，求解线性方程组得到X
A A(135,135)
X X (135,1)
B B(135,1)
高级材料
7
%Matlab 2D clear; clc; figure('color','w');
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) 4U (i, j) 0
高级材料
5
j+1 (i-1,j+1) (i,j+1) 2
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
相容性是比较容易满足的。在此基础上，如果满足了稳定性条件，差分格式的 收敛性就自动满足。
高级材料
4
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)
高级材料
2
稳定性（stability）：如果偏微分方程的严格解析解有界，差分格式给出的 解也有界，称该差分格式是稳定的；如果差分格式给出的解是无界的，则 称该差分格式是不稳定的。
稳定性反映了差分格式在计算中控制误差传递的能力
高级材料
3
收敛性（convergence）：如果当时间和空间步长趋于零时，FDE解趋于PDE
回顾
1. 有限差分法基础 2. 差分格式 3. 差分方程 4. 边界条件的处理 5. 相容性、稳定性和收敛性
高级材料
1
回顾
1. 有限差分法的相容性、稳定性和收敛性
相容性：针对差分格式而言，在时间步长和空间步长趋近于零的情况下， 如果差分格式的截断误差（差分格式与原有偏微分方程之差）的模趋近于零， 则该差分格式与原偏微分方程是相容的，或称该差分方程与原偏微分方程 具有相容性。
subplot(1,2,1),mesh(s) axis([0,17,0,11,0,100]) subplot(1,2,2),contour(s,32)
8
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
高级材料
9
2.5 应用实例
南加州一次未来大地震的强地面运动的数值模拟
内部节点：
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;
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
高级材料
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
a=zeros(135,135); for i=1:135
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;