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

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