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

小结
• 1、高阶精度的有限差分算子可以利用Taylor级数展 开推导出来。 • 2、对于双侧算子，权重系数随着算子长度的增加 迅速减小。 • 3、对于单侧算子，权重系数随着算子长度的增加 迅速增大。 • 4、以上给出的高阶算子提高了空间微分近似的精 度，我们还需要给出提高时间域外推的精度的方法。
四点一阶导数
f ( x)
1 ( f ( x 2dx) 8 f ( x dx) 8 f ( x dx) f ( x 2dx)) 12dx
四点三阶导数
f ( x)
1 ( f ( x 2dx) 2 f ( x dx) 2 f ( x dx) f ( x 2dx)) 3 2(dx)
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');
为什么需要高阶算子
以上的差分形式又称为微分算子(derivative operator)。在实际应用中，基于一阶 和二阶差分形式的差分格式的精度不能满足要求，在以上差分形式中，我们只 用到相邻网格节点上的值。 是否可以利用更多的信息提高微分算子的精度？ 用待定系数法构造精度较高的差分格式。
Taylor算子
回顾：
显式差分方法
隐式差分方法
n1 n1 n1 n n n sTi sTi 1 (2 2s)T i 1 sT i 1 (2 2s)T i sT i 1
Crank-Nicolson隐式差分
显式
最简单隐式
Crank-Nicolson 隐式
牛顿冷却问题不同差分格式的matlab程序
两式相加：
af bf (a b) f (a b) f dx
注意：上式对任意a，b都成立！
插值：
令
a b 1 a b 0
a 0.5 微分： a b 0 令 a b 1/ dx b 0.5
插值权重：
1 1 a ,b 2dx 2dx
牛顿冷却问题不同方法结果的对比
四阶R-K方法精度最高！
牛顿冷却问题不同方法结果的对比
四阶R-K方法精度最高！
预估
校正
for n=0,1,2,3,...,n-1 tn 1 tn dt ; k1 k2 dt
预估

dt
Tn (Tn k 1);
校正

1 Tn 1 Tn (k 1 k 2); 2 End
高阶外推
y f ( x, y ) x
显式Euler方法是最简单的单步法，它是一阶的，它可以看作Talylor展 开后取前两项。因此，得到高阶方法的一个直接想法是用Talylor展开， 如果能计算的高阶导数，则可写出p阶方法的计算方法
j 1,l
f(x-dx)
f(x)
f(x+dx)
一阶情况：
f ( x dx) f x f ( x)dx
f ( x dx) f x f ( x)dx
f(x-dx)
f(x)
f(x+dx)
1 两点一阶情况：
2
af af af dx
bf bf bf dx
1 1 , w2 2dx 2dx
插值权重：
w1 0.5, w2 0.5
w1
x-2dx x-dx 1 a 2 b
x
x+dx 3 c 4 d
x+2dx
四点三阶
x-2dx x-dx 1 a 2 b
x
x+dx 3 c 4 d
x+2dx
四点插值
f ( x)
1 2 2 1 f ( x 2dx) f ( x dx) f ( x dx) f ( x 2dx) 6 3 3 6
不同差分格式计算结果对比
混合差分格式精度最高！
不同差分格式计算结果对比
混合差分格式精度最高！
不同差分格式计算结果对比
混合差分格式精度最高！
不同差分格式计算结果对比
混合差分格式精度最高！
不同差分格式计算结果对比
混合差分格式精度最高！
高阶差分算子
主讲人：胡才博 中国科学院大学地球科学学院 中国科学院计算地球动力学重点实验室
显式结构
隐式结构பைடு நூலகம்
几类显式Runge-Kutta方法
对于L=2，则
其局部截断误差是
令
1
四阶R-K方法是最常见的高精度方法：
1 T j 1 T j (k1 2k2 2k3 k4 ) 6 k1 dtf (t j , T j ) k dt ,Tj 1 ) 2 2 k2 dt k3 dtf (t j , T j ) 2 2 k4 dtf (t j dt , T j k3 ) k2 dtf (t j
算子长度
微分
微分 函数f(x)在某一位置x的近似值以及导数（一阶、二阶、高阶等）如何计算？ 办法：充分利用更多相邻点的信息。 问题：如何计算相邻点的权重（系数）？
i-3
i-2
i-1
i
i+1
i+2
i+3
Taylor算子
构造如下形式的表达式：
在不同相邻点进行Taylor级数展开
f (i ) ( x) w(ji ) f ( x j )
董良国等，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
%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 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)*(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
1 T j 1 T j (k1 2k2 2k3 k4 ) 6 k1 dtf (t j , T j ) k dt ,Tj 1 ) 2 2 k2 dt k3 dtf (t j , T j ) 2 2 k4 dtf (t j dt , T j k3 ) k2 dtf (t j
dT f (T , t ) dt
四阶R-K有限差分解牛顿冷却问题
f (T , t )
T

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); end plot(xt,T_rk,'k.-');hold on % R_K_4 difference scheme
x-2dx x-dx 1 a 2 b
x 3 c
x+dx 4 d
x+2dx 5 e
五点二阶导数
f ( x)
1 ( f ( x 2dx) 16 f ( x dx) 30 f ( x) 16 f ( x dx) f ( x 2dx)) 2 12(dx)
单侧算子
y f ( x, y ) x
c=1, 得到Euler公式
h yn 1 yn f ( xn h), yn hf ( xn , yn ) 2
1 T j 1 T j (k1 2k2 2k3 k4 ) 6 k1 dtf (t j , T j ) k dt ,Tj 1 ) 2 2 k dt k3 dtf (t j , T j 2 ) 2 2 k4 dtf (t j dt , T j k3 ) k2 dtf (t j