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

《第六讲_2.3高阶有限差分(6).pptx》由会员分享，可在线阅读，更多相关《第六讲_2.3高阶有限差分(6).pptx（48页珍藏版）》请在一课资料网上搜索。

1、回顾回顾：显式差分方法隐式差分方法Crank-Nicolson隐式差分显式最简单隐式Crank-Nicolson 隐式1111111(22)(22)nnnnnniiiiiisTs TsTsTs TsT 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);ylabe

2、l(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 steps T_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)

3、*(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程序混合差分格式精度最高！不同差分格式计算结果对比混合差分格式精度最高！不同差分格式计算结果对比混合差分格式精度最高！不同差分格式计算结果对比混合差分格式精度最高！不同差分格式计算结果对比混合差分格式精度最高！不同差分格式计算结果对比高阶差分算子主讲人：胡才

4、博中国科学院大学地球科学学院中国科学院计算地球动力学重点实验室为什么需要高阶算子以上的差分形式又称为微分算子(derivative operator)。在实际应用中，基于一阶和二阶差分形式的差分格式的精度不能满足要求，在以上差分形式中，我们只用到相邻网格节点上的值。是否可以利用更多的信息提高微分算子的精度？用待定系数法构造精度较高的差分格式。Taylor算子微分微分算子长度函数f(x)在某一位置x的近似值以及导数（一阶、二阶、高阶等）如何计算？办法：充分利用更多相邻点的信息。问题：如何计算相邻点的权重（系数）？ii+1i+2i+3i-3i-2i-1Taylor算子在不同相邻点进行Taylor级

5、数展开f(x)f(x+dx)f(x-dx)()()f xdxf xfx dx()()f xdxf xfx dx构造如下形式的表达式：()()1,()()iijjjlfxwf x一阶情况：f(x)f(x+dx)f(x-dx)afafaf dxbfbfbf dx两点一阶情况：两式相加：()()afbfab fab f dx注意：上式对任意a，b都成立！插值：10abab令0.50.5ab插值权重：120.5,0.5ww微分：01/ababdx令11,22abdxdx 插值权重：1211,22wwdxdx 12xx+dxx-dx12x+2dxx-2dx34abcd四点三阶xx+dxx-dx12x+2

6、dxx-2dx34abcd1221()(2)()()(2)6336f xf xdxf xdxf xdxf xdx 四点插值1()(2)8()8()(2)12fxf xdxf xdxf xdxf xdxdx四点一阶导数31()(2)2()2()(2)2()fxf xdxf xdxf xdxf xdxdx四点三阶导数x-dx12x+2dxx-2dx34abcdxx+dx5e21()(2)16()30()16()(2)12()fxf xdxf xdxf xf xdxf xdxdx 五点二阶导数单侧算子单侧算子小结 1、高阶精度的有限差分算子可以利用Taylor级数展开推导出来。2、对于双侧算子，权重

7、系数随着算子长度的增加迅速减小。3、对于单侧算子，权重系数随着算子长度的增加迅速增大。4、以上给出的高阶算子提高了空间微分近似的精度，我们还需要给出提高时间域外推的精度的方法。实例：董良国等，2000，地球物理学报一阶弹性波方程交错网格高阶差分解法高阶Taylor外推高阶外推高阶外推for n=0,1,2,3,N-1 xn+1=xn+dx;k1=dx*f(xn,yn);k2=dx*f(xn+1,yn+k1);yn+1=yn+1/2*(k1+k2)end预估校正1121112 n=0,1,2,3,.,n-1 ;();1 ();2 nnnnnnforttdtdtkTdtkTkTTkkEnd 预估校

8、正高阶外推高阶外推显式Euler方法是最简单的单步法，它是一阶的，它可以看作方法是最简单的单步法，它是一阶的，它可以看作Talylor展展开后开后取前两项。因此，得到高阶方法的一个直接想法是用Talylor展开，展开，如果能计算如果能计算的高阶导数，则可写出p阶方法的计算方法阶方法的计算方法(,)yf x yx显式结构隐式结构几类显式Runge-Kutta方法方法对于L=2，则，则其局部截断误差是令11123411223431(22)6(,)(,)22(,)22(,)jjjjjjjjjjTTkkkkkdtf t Tkdtkdtf tTkdtkdtf tTkdtf tdt Tk四阶R-K方法是最

9、常见的高精度方法：(,)yf x yxc=1,得到Euler公式1(),(,)2nnnnnnhyyf xhyhf xy1123411223431(22)6(,)(,)22(,)22(,)jjjjjjjjjjTTkkkkkdtf t Tkdtkdtf tTkdtkdtf tTkdtf tdt Tk T_rk(1)=t0;for i=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.-);hold on%R_K_4 difference scheme(,)dTf T tdt(,)Tf T t 四阶R-K有限差分解牛顿冷却问题1123411223431(22)6(,)(,)22(,)22(,)jjjjjjjjjjTTkkkkkdtf t Tkdtkdtf tTkdtkdtf tTkdtf tdt Tk牛顿冷却问题不同方法结果的对比四阶R-K方法精度最高！牛顿冷却问题不同方法结果的对比四阶R-K方法精度最高！