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

1、回顾,有限差分法基础差分格式差分方程边界条件的处理相容性、稳定性和收敛性,回顾,1. 有限差分法的相容性、稳定性和收敛性,相容性：针对差分格式而言，在时间步长和空间步长趋近于零的情况下， 如果差分格式的截断误差（差分格式与原有偏微分方程之差）的模趋近于零， 则该差分格式与原偏微分方程是相容的，或称该差分方程与原偏微分方程 具有相容性。,稳定性（stability）：如果偏微分方程的严格解析解有界，差分格式给出的解也有界，称该差分格式是稳定的；如果差分格式给出的解是无界的，则称该差分格式是不稳定的。,稳定性反映了差分格式在计算中控制误差传递的能力,收敛性（convergence）：如果当时间和空

2、间步长趋于零时，FDE解趋于PDE解，称该差分格式是收敛的。 如果则称该差分格式是收敛的。,收敛性描述的是当差分网格无限细化时，差分方程的解是否具有无限逼近偏微分方程的解的能力,Lax等价定理（Lax equivalence theorem）：如果逼近一个给定问题的差分格式是相容的，那么该差分格式的收敛性与稳定性互为充分必要条件。,相容性是比较容易满足的。在此基础上，如果满足了稳定性条件，差分格式的收敛性就自动满足。,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

3、(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,内部节点：,边界节点：,A矩阵非零系数减少， 同时引入第一类边界， 方程右端项B向量出现 非零元素。,局部节点编号,总体节点编号,组建A和B矩阵，求解线性方程组得到X,%Matlab 2D clear; clc; figure(color,w);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.2

4、5;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(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+1

5、5)=-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,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:

6、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 endsubplot(1,2,1),mesh(s) axis(0,17,0,11,0,100) subplot(1,2,2),contour(s,32),2.5 应用实例,南加州一次未来大地震的强地面运动的数值模拟,盆地效应,Cui, 2013,Cui, 2013,Cui, 2013,Cui, 2013,总结：1、有限差分方法给出的数值解的精度取决于所用的差分形式（向 前、向后、中心）。2、偏微分方程的显式有限差分格式通常是有条件稳定的，为了保 证

7、得到精确的数值解，最关键的是需要根据稳定性条件选取正确的 空间和时间步长。,显式与隐式差分格式,主讲人：胡才博 中国科学院大学地球科学学院 中国科学院计算地球动力学重点实验室,显式差分格式（explicit difference scheme）差分方法中可逐层逐点分别求解的格式。特点 1. 不联立解方程； 2.时间步长和空间步长的选择受限制。通常要求时间步长足够小。,隐式差分格式（implicit difference scheme）特点 时间步长和空间步长的选择不受限制； 需要联立解方程组,显式和隐式：求解问题与时间相关,例子：,1. 显式差分格式：,左端：n+1时刻的值； 右端：n时刻的值

8、。,特点：结构简洁，直接求解，求解速度快。,但是，时间步长需满足：,显式差分格式才能得到稳定的数值解，否则，数值解将会不稳定而振荡。,显示差分格式示意图,2. 隐式差分格式：,时间一阶精度 空间二阶精度,隐式有限差分格式,Crank-Nicolson 隐式差分格式,Crank-Nicolson 隐式差分格式,Forward-Time Central-Space method,Backward -Time Central -Space method,Crank-Nicolson 隐式差分格式,一般差分格式,求解区域：,边界条件：,初始条件：,一种隐式差分格式的程序实现,A = sparse(nx

9、,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;,内部节点：,边界节点：,载荷项：,内部,边界,边界条件：,初始条件：,Crank-Nicolson 隐式差分格式的程序实现,A = sparse(nx,nx); for i=2:nx-1 A(i,i-1) = -s; A(i,i ) =

10、(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,B(1 ,1 ) = 1; B(nx,nx) = 1;,内部节点：,边界节点：,例子：牛顿冷却定律：温度高于周围环境的物体向周围媒质传递热量 逐渐冷却时所遵循的规律。当物体表面与周围存在温度差时，单位时 间从单位面积散失的热量与温度差成正比。,Tair,一阶常微分方程的数值解,首先对时间和温度进行离

11、散：,利用向前差分形式：,得到以下的显式差分格式：,Tcap,利用向前差分格式：,现在改用向后差分形式进行近似，得到隐式差分格式：,可以验证，当时间步长趋近于零时，以上近似解趋于解析解。因此， 该格式收敛。 稳定性条件：T单调减小的条件,显式差分格式,隐式差分格式,当dt=1.25，tau=0.7时，显式差分格式不稳定，结果振荡； 隐式差分格式稳定，结果不精确。,隐式差分格式,无条件稳定,重点考察差分格式的收敛性,当dt=1，tau=0.7时，显式差分格式不稳定，结果振荡； 隐式差分格式稳定，结果不精确。,当dt=0.5，tau=0.7时，显式差分格式稳定， 隐式差分格式稳定，结果不精确，两者

12、都不精确。,当dt=0.1，tau=0.7时，显式差分格式稳定； 隐式差分格式稳定；结果都比较精确。,当dt=0.01，tau=0.7时，显式差分格式稳定； 隐式差分格式稳定；结果都相当精确。,当dt和tau都大于零时，该式无条件满足，因此混合差分格式无条件稳定。,xt(nt+1)=nt*dt;plot(xt,T_e,b.-, xt,T_i,g.-, xt,T_m,m.-, xt,T_a,r.-,);hold on % set(gca,DataAspectRatio,(max(xt)-min(xt)/(max(T_e)-min(T_e)/3 1 1);xlabel(Time (s),Fontn

13、ame,times new roman,FontSize,14);ylabel(Temperature,Fontname,times new roman,FontSize,14);title(dt=0.01 tau=0.7);,%Malab-1D clear; clc; figure(color,w);t0=1; % initial temperature tau=0.7; % time constant dt=0.01; % time interval t_total=10; nt=round(t_total/dt); % total time stepsT_e(1)=t0; T_i(1)=

14、t0; T_m(1)=t0; for i=1:nt;xt(i)=(i-1)*dt; T_e(i+1)=T_e(i)*(1-dt/tau); %explicit T_i(i+1)=T_i(i)/(1+dt/tau); % implicit T_m(i+1)=T_m(i)*(1-dt/2/tau)/(1+dt/2/tau); %mix T_a(i)=t0*exp(-xt(i)/tau); %analytical results end,不同差分格式的matlab程序,混合差分格式精度最高！,不同差分格式计算结果对比,混合差分格式精度最高！,不同差分格式计算结果对比,混合差分格式精度最高！,不同差分

15、格式计算结果对比,混合差分格式精度最高！,不同差分格式计算结果对比,混合差分格式精度最高！,不同差分格式计算结果对比,显式差分格式,1.对步长有要求； 2.无需解方程,二阶精度,%Matlab 2D clear; clc; figure(color,w);lx=17;ly=11; % v1=zeros(ly,lx); % for j=2:lx-1 v1(ly,j)=100; end % v2=v1;maxt=1;t=0; k=0;,while(maxt1e-6) % k=k+1 maxt=0; for i=2:ly-1 for j=2:lx-1; v2(i,j)=(v1(i,j+1)+v1(i

16、+1,j)+v1(i-1,j)+v1(i,j-1)/4; t=abs(v2(i,j)-v1(i,j); if(tmaxt) maxt=t;end end end v1=v2; end %,subplot(1,2,1),mesh(v1) axis(0,17,0,11,0,100) subplot(1,2,2),contour(v1,32),迭代解法:,K:迭代步数,K=419,总结：显式格式算法简单、易于编程，可以从给定的初始条件开始，在时间上逐层前进求解。一些与时间有关的偏微分方程的求解，需要用到隐式差分格式，在时间上计算数值解的传播时，需要求解线性方程组。通常在计算的每一个时间步，需要求矩阵

17、的逆矩阵。因此，隐式格式算法相对于显式格式更复杂，编程更困难。显式格式通常比隐式格式的稳定性差，如果时间步长取得过大，可能会给出物理上不正确的结果。某些隐式格式的优点是其无条件稳定性，因此时间步长可以取得大一些，但是并不能保证精度很高。可以利用显式-隐式混合格式（如： Crank-Nicholson scheme ），它们无条件稳定，而且精度高于相应的显式和隐式格式。,总结有限差分法的主要内容包括：1、根据问题的特点将定解区域作网格剖分； 在所有网格节点上用有限差分格式对导数求近似，对函数、初始和边界条件求近似； 把原方程离散化为代数方程组，即有限差分方程组. 2、有限差分模型形态的理论分析，

18、以保证计算过程可行及计算结果正确： 解的相容性：对于一个微分方程建立的各种差分格式，为了有实用意义，一个基本要求是它们能够任意逼近微分方程，这就是相容性要求。 解的稳定性：因为差分格式的计算过程是逐层推进的，在计算第m1层的近似值时要用到第m层的近似值，直到与初始值有关。前面各层若有舍入误差，必然影响到后面各层的值，如果误差的影响越来越大，以致差分格式的精确解的面貌完全被掩盖，这种格式是不稳定的，相反如果误差的传播是可以控制的，就认为格式是稳定的。只有在这种情形下，差分格式在实际计算中的近似解才可能任意逼近差分方程的精确解。 解的收敛性：一个差分格式是否有用，最终要看差分方程的精确解能否任意逼

19、近微分方程的解，这就是收敛性的概念。,总结：3、差分格式：对于有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式。考虑时间因子的影响，差分格式还可以分为显式、隐式、显隐交替格式等。目前常见的差分格式，主要是上述几种形式的组合，不同的组合构成不同的差分格式。通过对时间和空间这几种不同差分格式的组合，可以组合成不同的差分计算格式。显式格式算法简单、易于编程，可以从给定的初始条件开始，在时间上、逐层前进求解。隐式格式需要联立求解线性方程组。,总结： 4、差分格式的构造：构造差分的方法有多种形式，目前主要采用的是泰勒级数展开方法。其基本的差分表达式主要有三种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等， 其中前两种格式为一阶计算精度，后两种格式为二阶计算精度。 5、差分网格的步长一般根据稳定性条件来决定。,分别用显式、最简单隐式以及Crank-Nicolson隐式差分格式求解如下定解问题（统一选取： ，计算100步）：,作业二,要求：1. 编制相关计算程序（建议使用matlab）；2. 绘制第100步时不同方法得到Tx曲线（画在一张图上）；3. 分别绘制第2、3、4、5点不同方法得到的T-t对比曲线（2、3、4、5点各一张图）,