第六讲_2.3高阶有限差分(6)

回顾：显式差分方法隐式差分方法Crank-Nicolson隐式差分显式最简单隐式Crank-Nicolson隐式1111111(22)(22)nnnnnniiiiiisTsTsTsTsTsTxt(nt+1)=nt*dt;plot(xt,T_e,'b.-',xt,T_i,'g.-',xt,T_m,'m.-',xt,T_a,'r.-',);holdon%set(gca,'DataAspectRatio',[(max(xt)-min(xt))/(max(T_e)-min(T_e))/311]);xlabel('Time(s)','Fontname','timesnewroman','FontSize',14);ylabel('Temperature','Fontname','timesnewroman','FontSize',14);title('dt=0.01tau=0.7');%Malab-1Dclear;clc;figure('color','w');t0=1;%initialtemperaturetau=0.7;%timeconstantdt=0.01;%timeintervalt_total=10;nt=round(t_total/dt);%totaltimestepsT_e(1)=t0;T_i(1)=t0;T_m(1)=t0;fori=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);%mixT_a(i)=t0*exp(-xt(i)/tau);%analyticalresultsend牛顿冷却问题不同差分格式的matlab程序混合差分格式精度最高！不同差分格式计算结果对比混合差分格式精度最高！不同差分格式计算结果对比混合差分格式精度最高！不同差分格式计算结果对比混合差分格式精度最高！不同差分格式计算结果对比混合差分格式精度最高！不同差分格式计算结果对比高阶差分算子主讲人：胡才博中国科学院大学地球科学学院中国科学院计算地球动力学重点实验室为什么需要高阶算子以上的差分形式又称为微分算子(derivativeoperator)。在实际应用中，基于一阶和二阶差分形式的差分格式的精度不能满足要求，在以上差分形式中，我们只用到相邻网格节点上的值。是否可以利用更多的信息提高微分算子的精度？用待定系数法构造精度较高的差分格式。Taylor算子微分微分算子长度函数f(x)在某一位置x的近似值以及导数（一阶、二阶、高阶等）如何计算？办法：充分利用更多相邻点的信息。问题：如何计算相邻点的权重（系数）？ii+1i+2i+3i-3i-2i-1Taylor算子在不同相邻点进行Taylor级数展开f(x)f(x+dx)f(x-dx)()()fxdxfxfxdx()()fxdxfxfxdx构造如下形式的表达式：()()1,()()iijjjlfxwfx一阶情况：f(x)f(x+dx)f(x-dx)afafafdxbfbfbfdx两点一阶情况：两式相加：()()afbfabfabfdx注意：上式对任意a，b都成立！插值：10abab令0.50.5ab插值权重：120.5,0.5ww微分：01/ababdx令11,22abdxdx插值权重：1211,22wwdxdx12xx+dxx-dx12x+2dxx-2dx34abcd四点三阶xx+dxx-dx12x+2dxx-2dx34abcd1221()(2)()()(2)6336fxfxdxfxdxfxdxfxdx四点插值1()((2)8()8()(2))12fxfxdxfxdxfxdxfxdxdx四点一阶导数31()((2)2()2()(2))2()fxfxdxfxdxfxdxfxdxdx四点三阶导数x-dx12x+2dxx-2dx34abcdxx+dx5e21()((2)16()30()16()(2))12()fxfxdxfxdxfxfxdxfxdxdx五点二阶导数单侧算子小结•1、高阶精度的有限差分算子可以利用Taylor级数展开推导出来。•2、对于双侧算子，权重系数随着算子长度的增加迅速减小。•3、对于单侧算子，权重系数随着算子长度的增加迅速增大。•4、以上给出的高阶算子提高了空间微分近似的精度，我们还需要给出提高时间域外推的精度的方法。实例：董良国等，2000，地球物理学报一阶弹性波方程交错网格高阶差分解法高阶Taylor外推高阶外推forn=0,1,2,3,…,N-1xn+1=xn+dx;k1=dx*f(xn,yn);k2=dx*f(xn+1,yn+k1);yn+1=yn+1/2*(k1+k2)end预估校正1121112n=0,1,2,3,...,n-1;();1();2nnnnnnforttdtdtkTdtkTkTTkkEnd预估校正高阶外推显式Euler方法是最简单的单步法，它是一阶的，它可以看作Talylor展开后取前两项。因此，得到高阶方法的一个直接想法是用Talylor展开，如果能计算的高阶导数，则可写出p阶方法的计算方法(,)yfxyx显式结构隐式结构几类显式Runge-Kutta方法对于L=2，则其局部截断误差是令11123411223431(22)6(,)(,)22(,)22(,)jjjjjjjjjjTTkkkkkdtftTkdtkdtftTkdtkdtftTkdtftdtTk四阶R-K方法是最常见的高精度方法：(,)yfxyxc=1,得到Euler公式1(),(,)2nnnnnnhyyfxhyhfxy1123411223431(22)6(,)(,)22(,)22(,)jjjjjjjjjjTTkkkkkdtftTkdtkdtftTkdtkdtftTkdtftdtTkT_rk(1)=t0;fori=1:nt;xt(i)=(i-1)*dt;k1=-dt*T_rk(i)/tau;k2=-dt*(T_rk(i)+k1/2)/tau;k3=-dt*(T_rk(i)+k2/2)/tau;k4=-dt*(T_rk(i)+k3)/tau;T_rk(i+1)=T_rk(i)+1/6*(k1+2*k2+2*k3+k4);endplot(xt,T_rk,'k.-');holdon%R_K_4differencescheme(,)dTfTtdt(,)TfTt四阶R-K有限差分解牛顿冷却问题1123411223431(22)6(,)(,)22(,)22(,)jjjjjjjjjjTTkkkkkdtftTkdtkdtftTkdtkdtftTkdtftdtTk牛顿冷却问题不同方法结果的对比四阶R-K方法精度最高！牛顿冷却问题不同方法结果的对比四阶R-K方法精度最高！