差分.docx
《差分.docx》由会员分享，可在线阅读，更多相关《差分.docx（14页珍藏版）》请在冰点文库上搜索。
差分
问题
（1）利用中心差分格式近似导数d^2y/dx^2，数值求解常微分方程的边值问题：
空间步长取Δx=0.01
我们容易得到方程的解析解为
可以将有限差分的结果和解析解进行对比，检查计算误差。
用有限差分法求解的步骤：
（A）对应微分方程的中心差分格式为
边界条件为
U0=0；UN=1；
根据上述差分格式和边界条件，可构建关于未知变量Uj（j=1,2,..N-1）的线性方程组。
方程组系数为（N-1）*（N-1）的三对角阵。
（B）用追赶法解三对角线性方程组
有限差分实现常微分方程的MATLAB程序：
dx=0.01， 差分网格空间步长
x=dx:
dx:
1-dx;网格节点坐标
N=length（x）数据点长度
a三对角矩阵的对角线元素
b三对角矩阵对角线上的元素
c三对角矩阵对角线下的元素
u解向量
f节点上的自由项函数值（线性方程组Ax=f的右项）
%%%%--------------------------------------------
主程序：
%%-------------------
dx=0.01;
x=dx:
dx:
1-dx;
%%-------------------
N=length（x）;
a=-2*ones（N,1）;
b=ones（N-1,1）;;
c=ones（N-1,1）;
u=zeros（N,1）;
f=dx^2*sin（2*x）;
f
（1）=f
（1）;
f（N）=f（N）-1;
%%%%%%%%%%%%%%%
u=my_chasing（a,b,c,f）;
%%%%%%%%%%%%
y_an=-sin（2*x）/4+（1+sin
（2）/4）*x;%theanalysissolutionofODE
subplot（2,1,1）
plot（x,u,x,y_an,'.'）
title（['\deltax=',num2str（dx）]）
subplot（2,1,2）
plot（x,u-y_an'）
mean_error=mean（abs（u-y_an'））;
title（['mean_-error=',num2str（mean_error）]）
问题
（2）利用中心差分格式离散
，分别调用FTCS，BTCS和C-N格式。
求解偏微分方程的混合定解问题：
（一维热传导方程）
T（x,0）=0;（0<=x<=1）
空间步长取Δx=0.01
偏微分方程定解问题的解析解为
FTCS
步骤：
（A）写出对应微分方程的FTCS差分格式
（B）构建递推公式,递推得到每一层的结果
根据上面的双层差分格式，由第（n-1）层的结果递推n曾的结果。
对所有的1<=n<=N，改写上式
（1）为
其中λ=τ/h^2为网格比网格比
由初始值（t=0）便可递推得到任何时刻的温度值
FTCS差分格式的MATLAB程序：
程序中主要符号：
dx=0.01...，空间步长
x=dx:
dx:
1-dx;空间网格节点坐标
J=length（x）空间数据点长度
dt=0.25E-4 时间步长
N=2000 时间数据点长度
t=（0：
N-1）*dt时间节点坐标
ubc_l左侧边界条件
ubc_r 右侧边界条件
ubc_initial初始条件
lamda网格比
u解矩阵
f节点上的自由项函数值
%%%%--------------------------------------------
%%%%---------
dx=0.01;
x=dx:
dx:
1-dx;
J1=length（x）;
dt=0.25e-4;
N=2000;
t=（0:
N-1）*dt;
%%%%---------
%%Initialandboundarycondiction
ubc_l=100*ones（N,1）;
ubc_r=100*ones（N,1）;
u_initial=zeros（1,J1）;
%%%%---------
lamda=dt/dx^2;
u=zeros（N,J1）;
u（1,:
）=u_initial;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
forn1=2:
N
j1=1;
u（n1,j1）=u（n1-1,j1）+lamda*（u（n1-1,j1+1）-2*u（n1-1,j1）+ubc_l（n1-1））;
forj1=2:
J1-1
u（n1,j1）=u（n1-1,j1）+lamda*（u（n1-1,j1+1）-2*u（n1-1,j1）+u（n1-1,j1-1））;
end
j1=J1;
u（n1,j1）=u（n1-1,j1）+lamda*（ubc_r（n1-1）-2*u（n1-1,j1）+u（n1-1,j1-1））;
end
[xx,tt]=meshgrid（x,t（1:
10:
end））;
subplot（2,1,1）
contourf（xx,tt,u（1:
10:
end,:
））
colorbar
title（['\lambda=',num2str（lamda）]）
%%TheanalysissolutionofPDE
k3=find（t==max（t））;
T0=100;
fork2=1:
100
Dn=-200/k2/pi*（1-（-1）^k2）;
temp3=Dn*sin（k2*pi*x）*exp（-k2^2*pi^2*t（k3））;
T0=T0+temp3;
end
subplot（2,1,2）
plot（x',T0,x',u（k3,:
）,'*'）
mean_error=mean（abs（T0-u（k3,:
）））;
title（['t=',num2str（t（k3））,' mean_-error=',num2str（mean_error）]）
%%%%--------------------------------------------
上图中上侧图横坐标是x，纵坐标是时间，表示随着时间的延长，高温逐渐向内部扩散。
下侧图是最后时间步解析解与差分解的比较。
横坐标表示x。
可见平均误差为0.0054。
注意：
FTCS是条件稳定的差分格式，稳定条件为网格比λ<=0.5。
故在Δx=0.01的条件下，Δt应小于0.5E-4
BTCS差分格式（C-N差分格式与之类似，其步骤不另写）：
步骤：
（A）写出对应微分方程的BTCS差分格式
（B）构建递推公式,由第n-1层结果递推得到第n层的结果
递推公式为
差分格式为隐式格式，故对每一层，需利用边界条件，上一层结果求解。
求解时需用到追赶法。
BTCS差分格式的MATLAB程序：
程序中主要符号：
dx=0.01...，空间步长
x=dx:
dx:
1-dx;空间网格节点坐标
J=length（x）空间数据点长度
dt=0.25E-3 时间步长
N=500 时间数据点长度
t=（0：
N-1）*dt时间节点坐标
ubc_l左侧边界条件
ubc_r 右侧边界条件
ubc_initial初始条件
lamda网格比
u 解矩阵
f 节点上的自由项函数值
%%%%--------------------------------------------
主程序：
%%%%---------
dx=0.01;
x=dx:
dx:
1-dx;
J1=length（x）;
dt=0.25e-3;
N=500;
t=（0:
N-1）*dt;
%%%%---------
%%Initialandboundarycondiction
ubc_l=100*ones（N,1）;
ubc_r=100*ones（N,1）;
u_initial=zeros（1,J1）;
%%%%---------
lamda=dt/dx^2;
u=zeros（N,J1）;
u（1,:
）=u_initial;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%a:
themaindiagnalvectorofcoefficientmatrixA;Au=f
%b:
theupperdiagnalvectorofA
%c:
thelowerdiagnalvectorofA
%f:
Au=f
a=ones（J1,1）*（2*lamda+1）;
b=-lamda*ones（J1-1,1）;
c=-lamda*ones（J1-1,1）;
f=zeros（J1,1）;
%%%%%%%%%
forn1=2:
N
f=u（n1-1,:
）;
f
（1）=f
（1）+lamda*ubc_l（n1）;
f（J1）=f（J1）+lamda*ubc_r（n1）;
u（n1,:
）=my_chasing（a,b,c,f）;
end
[xx,tt]=meshgrid（x,t（1:
10:
end））;
subplot（2,1,1）
contourf（xx,tt,u（1:
10:
end,:
））
colorbar
title（['\lambda=',num2str（lamda）]）
%%TheanalysissolutionofPDE
k3=find（t==max（t））;
T0=100;
fork2=1:
100
Dn=-200/k2/pi*（1-（-1）^k2）;
temp3=Dn*sin（k2*pi*x）*exp（-k2^2*pi^2*t（k3））;
T0=T0+temp3;
end
subplot（2,1,2）
plot（x',T0,x',u（k3,:
）,'*'）
mean_error=mean（abs（T0-u（k3,:
）））;
title（['t=',num2str（t（k3））,' mean_-error=',num2str（mean_error）]）
求解三对角方程组的追赶法程序请见另一篇博文
解线性方程组的直接法及matlab程序
（2）
%%%%--------------------------------------------
下图为差分方程计算结果。
BTCS差分格式为无条件稳定的差分格式，从下图可见，网格比为2.5时，计算结果仍然稳定收敛。
另外，BTCS差分格式误差比FTCS差分格式误差大。
C-N差分格式的MATLAB程序：
对应微分方程的的C-N差分格式为：
其递推公式为
根据上面的双层差分格式，由第（n-1）层的结果递推n层的结果。
程序中主要符号：
与BTCS差分格式相同。
%%%%--------------------------------------------
主程序：
%%%%---------
dx=0.01;
x=dx:
dx:
1-dx;
J1=length（x）;
dt=0.25e-3;
N=500;
t=（0:
N-1）*dt;
%%%%---------
%%Initialandboundarycondiction
ubc_l=100*ones（N,1）;
ubc_r=100*ones（N,1）;
u_initial=zeros（1,J1）;
%%%%---------
lamda=dt/dx^2;
u=zeros（N,J1）;
u（1,:
）=u_initial;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%a:
themaindiagnalvectorofcoefficientmatrixA;Au=f
%b:
theupperdiagnalvectorofA
%c:
thelowerdiagnalvectorofA
%f:
Au=f
a=ones（J1,1）*（lamda+1）;
b=-lamda/2*ones（J1-1,1）;
c=-lamda/2*ones（J1-1,1）;
f=zeros（J1,1）;
%%%%%%%%%
forn1=2:
N
temp1=lamda/2*ubc_l（n1-1）+（1-lamda）*u（n1-1,1）+lamda/2*u（n1-1,2）;
f
（1）=temp1+lamda/2*ubc_l（n1）;
f（2:
J1-1）=lamda/2*u（n1-1,1:
J1-2）+（1-lamda）*u（n1-1,2:
J1-1）+lamda/2*u（n1-1,3:
J1）;
f（J1）=lamda/2*u（n1-1,J1-1）+（1-lamda）*u（n1-1,J1）+lamda/2*ubc_r（n1-1）+lamda/2*ubc_r（n1）;
u（n1,:
）=my_chasing（a,b,c,f）;
end
[xx,tt]=meshgrid（x,t（1:
10:
end））;
subplot（2,1,1）
contourf（xx,tt,u（1:
10:
end,:
））
colorbar
title（['\lambda=',num2str（lamda）]）
%%TheanalysissolutionofPDE
k3=find（t==max（t））;
T0=100;
fork2=1:
100
Dn=-200/k2/pi*（1-（-1）^k2）;
temp3=Dn*sin（k2*pi*x）*exp（-k2^2*pi^2*t（k3））;
T0=T0+temp3;
end
subplot（2,1,2）
plot（x',T0,x',u（k3,:
）,'*'）
mean_error=mean（abs（T0-u（k3,:
）））;
title（['t=',num2str（t（k3））,' mean_-error=',num2str（mean_error）]）
求解三对角方程组的追赶法程序请见另一篇博文
解线性方程组的直接法及matlab程序
（2）
%%%%--------------------------------------------
下图为C-N差分格式的数值结果。
C-N为无条件稳定的差分格式。
从下图可见，C-N差分格式的精度远高于BTCS差分格式。
且为无条件稳定的差分格式。
升级会员 