差分.docx
- 文档编号:6743332
- 上传时间:2023-01-09
- 格式:DOCX
- 页数:12
- 大小:167.19KB
差分.docx
《差分.docx》由会员分享,可在线阅读,更多相关《差分.docx(12页珍藏版)》请在冰豆网上搜索。
差分
问题
(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差分格式。
且为无条件稳定的差分格式。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 差分.docx