差分Word文件下载.docx
《差分Word文件下载.docx》由会员分享，可在线阅读，更多相关《差分Word文件下载.docx（14页珍藏版）》请在冰点文库上搜索。
a三对角矩阵的对角线元素
b三对角矩阵对角线上的元素
c三对角矩阵对角线下的元素
u解向量
f节点上的自由项函数值（线性方程组Ax=f的右项）
%%%%--------------------------------------------
主程序：
%%-------------------
dx=0.01;
%%-------------------
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'
））;
mean_-error='
num2str（mean_error）]）
问题
（2）利用中心差分格式离散
，分别调用FTCS，BTCS和C-N格式。
求解偏微分方程的混合定解问题：
（一维热传导方程）
T（x,0）=0;
（0<
=x<
=1）
偏微分方程定解问题的解析解为
FTCS
步骤：
（A）写出对应微分方程的FTCS差分格式
（B）构建递推公式,递推得到每一层的结果
根据上面的双层差分格式，由第（n-1）层的结果递推n曾的结果。
对所有的1<
=n<
=N，改写上式
（1）为
其中λ=τ/h^2为网格比网格比
由初始值（t=0）便可递推得到任何时刻的温度值
FTCS差分格式的MATLAB程序：
程序中主要符号：
dx=0.01...，空间步长
空间网格节点坐标
J=length（x）空间数据点长度
dt=0.25E-4
时间步长
N=2000
时间数据点长度
t=（0：
N-1）*dt时间节点坐标
ubc_l左侧边界条件
ubc_r
右侧边界条件
ubc_initial初始条件
lamda网格比
u解矩阵
f节点上的自由项函数值
%%%%---------
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））;
contourf（xx,tt,u（1:
end,:
））
colorbar
\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;
plot（x'
T0,x'
u（k3,:
）,'
*'
mean_error=mean（abs（T0-u（k3,:
）））;
t='
num2str（t（k3））,'
上图中上侧图横坐标是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程序：
dt=0.25E-3
N=500
u
解矩阵
f
节点上的自由项函数值
%%%%--------------------------------------------
dt=0.25e-3;
N=500;
%a:
themaindiagnalvectorofcoefficientmatrixA;
Au=f
%b:
theupperdiagnalvectorofA
%c:
thelowerdiagnalvectorofA
%f:
a=ones（J1,1）*（2*lamda+1）;
b=-lamda*ones（J1-1,1）;
c=-lamda*ones（J1-1,1）;
f=zeros（J1,1）;
%%%%%%%%%
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）;
求解三对角方程组的追赶法程序请见另一篇博文
解线性方程组的直接法及matlab程序
（2）
下图为差分方程计算结果。
BTCS差分格式为无条件稳定的差分格式，从下图可见，网格比为2.5时，计算结果仍然稳定收敛。
另外，BTCS差分格式误差比FTCS差分格式误差大。
C-N差分格式的MATLAB程序：
对应微分方程的的C-N差分格式为：
其递推公式为
根据上面的双层差分格式，由第（n-1）层的结果递推n层的结果。
与BTCS差分格式相同。
a=ones（J1,1）*（lamda+1）;
b=-lamda/2*ones（J1-1,1）;
c=-lamda/2*ones（J1-1,1）;
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）;
下图为C-N差分格式的数值结果。
C-N为无条件稳定的差分格式。
从下图可见，C-N差分格式的精度远高于BTCS差分格式。
且为无条件稳定的差分格式。
升级会员 