第五讲——显式差分和隐式差分.pptx

回顾1.有限差分法基础2.差分格式3.差分方程4.边界条件的处理5.相容性、稳定性和收敛性回顾1. 有限差分法的相容性、稳定性和收敛性相容性 ： 针对差分格式而言 ，在时间步长和空间步长趋近于零的情况下，如果差分格式的截断误差（差分格式与原有偏微分方程之差）的模趋近于零，则该差分格式与原偏微分方程是相容的，或称该差分方程与原偏微分方程具有相容性稳定性（ stability）：如果偏微分方程的 严格解析解有界 ， 差分格式给出的解也有界 ，称该差分格式是稳定的；如果差分格式给出的解是无界的，则称该差分格式是不稳定的稳定性 反映了差分格式在计算中控制误差传递的能力收敛性（ convergence） ：如果当时间和空间步长趋于零时， FDE解趋于 PDE解，称该差分格式是收敛的如果则称该差分格式是收敛的收敛性 描述的是当差分网格无限细化时，差分方程的解是否具有无限逼近偏微分方程的解的能力Lax等价定理（ Lax equivalence theorem）：如果逼近一个给定问题的差分格式是相容的，那么该差分格式的收敛性与稳定性互为充分必要条件相容性是比较容易满足的在此基础上，如果满足了稳定性条件，差分格式的收敛性就自动满足。

U=0U=0U=100U=021 43 65 87 109 1211 1413 152.5 有限差分法实例(i,j) (i+1,j)(i-1,j)0 1234(i,j) (i+1,j-1)(i-1,j-1)(i,j+1) (i+1,j+1)(i-1,j+1)i-1 i i+1j-1jj+1h1h3h2h4(i,j) (i+1,j)(i-1,j)0 1234(i,j) (i+1,j-1)(i-1,j-1)(i,j+1) (i+1,j+1)(i-1,j+1)i-1 i i+1j-1jj+1h1h3h2h4for j=2:n-1for 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; endend内部节点： 边界节点：A矩阵非零系数减少，同时引入第一类边界，方程右端项 B向量出现非零元素局部节点编号总体节点编号组建 A和 B矩阵，求解线性方程组得到 X%Matlab 2Dclear;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.25;a(15*i+1,15*i-14)=-0.25;endfor 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; Enda(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+15)=-0.25; endfor i=122:134 a(i,i-1)=-0.25;a(i,i+1)=-0.25; a(i,i-15)=-0.25; endfor i=1:7for 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;endendb=a^(-1);c=zeros(135,1);for i=121:135 c(i,1)=25;endd=b*c;s=zeros(11,17);for i=2:16 s(11,i)=100; endfor i=1:9for j=1:15;s(i+1,j+1)=d(15*(i-1)+j,1);endendsubplot(1,2,1),mesh(s)axis([0,17,0,11,0,100])subplot(1,2,2),contour(s,32)2.5 应用实例南加州一次未来大地震的强地面运动的数值模拟盆地效应Cui, 2013Cui, 2013Cui, 2013Cui, 2013总结：1、有限差分方法给出的数值解的精度取决于所用的差分形式（向前、向后、中心）。

2、偏微分方程的显式有限差分格式通常是有条件稳定的，为了保证得到精确的数值解，最关键的是需要根据稳定性条件选取正确的空间和时间步长显式与隐式差分格式主讲人：胡才博中国科学院大学地球科学学院中国科学院计算地球动力学重点实验室• 显式差分格式 （ explicit difference scheme）差分方法中可逐层逐点分别求解的格式• 特点• 1. 不联立解方程；• 2.时间步长和空间步长的选择受限制通常要求时间步长足够小隐式差分格式 （ implicit difference scheme）特点1. 时间步长和空间步长的选择不受限制；2. 需要联立解方程组显式 和隐式：求解问题与时间相关例子：1. 显式差分格式：左端： n+1时刻的值；右端： n时刻的值 特点：结构简洁，直接求解，求解速度快但是，时间步长需满足：显式差分格式才能得到稳定的数值解，否则，数值解将会不稳定而振荡显示差分格式示意图2. 隐式差分格式：时间一阶精度空间二阶精度隐式有限差分格式Crank-Nicolson 隐式差分格式Crank-Nicolson 隐式差分格式Forward-Time Central-Space methodBackward -Time Central -Space methodCrank-Nicolson 隐式差分格式一般差分格式31 2 54 6求解区域： 边界条件：初始条件：一种隐式差分格式的程序实现A = sparse(nx,nx);for i=2:nx-1A(i,i-1) = -s;A(i,i ) = (1+2*s);A(i,i+1) = -s;endA(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;内部节点：边界节点：载荷项：内部边界31 2 54 6边界条件： 初始条件：Crank-Nicolson 隐式差分格式的程序实现A = sparse(nx,nx);for i=2:nx-1A(i,i-1) = -s;A(i,i ) = (2+2*s);A(i,i+1) = -s;endA(1 ,1 ) = 1;A(nx,nx) = 1;内部节点：边界节点：B= sparse(nx,nx);for i=2:nx-1B(i,i-1) = s;B(i,i ) = (2-2*s);B(i,i+1) = s;endB(1 ,1 ) = 1;B(nx,nx) = 1;内部节点：边界节点：例子：牛顿冷却定律：温度高于周围环境的物体向周围媒质传递热量逐渐冷却时所遵循的规律。

当物体表面与周围存在温度差时，单位时间从单位面积散失的热量与温度差成正比Tair一阶常微分方程的数值解首先对时间和温度进行离散：利用向前差分形式：得到以下的显式差分格式：Tcap利用向前差分格式：现在改用向后差分形式进行近似，得到隐式差分格式：可以验证，当时间步长趋近于零时，以上近似解趋于解析解因此，该格式收敛稳定性条件： T单调减小的条件显式差分格式隐式差分格式当 dt=1.25， tau=0.7时，显式差分格式不稳定，结果振荡；隐式差分格式稳定，结果不精确隐式差分格式无条件稳定重点考察差分格式的收敛性当 dt=1， tau=0.7时，显式差分格式不稳定，结果振荡；隐式差分格式稳定，结果不精确当 dt=0.5， tau=0.7时，显式差分格式稳定，隐式差分格式稳定，结果不精确，两者都不精确当 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)','Fontname','times new roman','FontSize',14);ylabel('Temperature','Fontname','times new roman','FontSize',14);title('dt=0.01 tau=0.7');%Malab-1Dclear;clc;figure('color','w');t0=1; % initial temperaturetau=0.7; % time constantdt=0.01; % time intervalt_total=10;nt=round(t_total/dt); % total time stepsT_e(1)=t0; T_i(1)=t0; T_m(1)=t0; for i=1:nt;xt(i)=(i-1)*dt;T_e(i+1)=T_e(i)*(1-dt/tau); %explicitT_i(i+1)=T_i(i)/(1+dt/tau); % implicitT_m(i+1)=T_m(i)*(1-dt/2/tau)/(1+dt/2/tau); %mix T_a(i)=t0*exp(-xt(i)/tau); %analytical resultsend不同差分格式的 matlab程序混合差分格式精度最高！不同差分格式计算结果对比混合差分格式精度最高！不同差分格式计算结果对比混合差分格式精度最高！不同差分格式计算结果对比混合差分格式精度最高！不同差分格式计算结果对比混合差分格式精度最高！不同差分格式计算结果对比显式差分格式1.对步长有要求；2.无需解方程二阶精度%Matlab 2Dclear;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(maxt>1e-6) %k=k+1maxt=0;for i=2:ly-1for j=2:lx-1;v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v1(i-1,j)+v1(i,j-1))/4; t=abs(v2(i,j)-v1(i,j));if(t>ma。