有限差分法.doc

有限差分法

一、单变量函数：

用中心差分法（matlab程序见附录）计算结果如下：

图1 中心差分法

图2 绝对误差

表1 数据对比

　x
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
精确解y
0
0.070467
0.14264
0.21824
0.29903
0.38682
0.48348
0.59099
0.71141
0.84696
1
数值解y
0
0.070489
0.14268
0.2183
0.29911
0.3869

0.48357
0.59107
0.71148
0.847
1
绝对误差
0
2.2E-05
4E-05
6E-05
8E-05
8E-05
9E-05
8E-05
7E-05
4E-05
0

小结：可见由中心差分法解本算例，其绝对误差数量级为

二、一维热传导：

在此取φ(x)=0，g1(t)= g2(t)=100-100*exp(-t)

问题描述：
已知厚度为l的无限大平板，初温0度 ，初始瞬间将其放于温度为100度的流体中，流体与板面间的表面传热系数为一常数。
试确定在非稳态过程中板内的温度分布。
(1) 显式差分法：

图3 显式差分法

(2) 隐式差分法：

图4 隐式差分法

小结：显式格式仅当时格式是稳定的。（其中称为网格比）
隐式格式从k层求k+1层时，需要求解一个阶方程组。而且隐式格式的稳定性对网格比没有要求，即为绝对稳定的。

三、Possion方程：

取f=1,R=1

图5差分法

图6 误差

小结：观察误差曲面，其绝对误差数量级为
附Matlab程序：

第1题：
%===========================Boundary Value Problem 1
clear;clc;
A=[-2.01 1 0 0 0 0 0 0 0;
1 -2.01 1 0 0 0 0 0 0;
0 1 -2.01 1 0 0 0 0 0;
0 0 1 -2.01 1 0 0 0 0;
0 0 0 1 -2.01 1 0 0 0;
0 0 0 0 1 -2.01 1 0 0;
0 0 0 0 0 1 -2.01 1 0;
0 0 0 0 0 0 1 -2.01 1;
0 0 0 0 0 0 0 1 -2.01;];
c1=[0.1;0.2;0.3;0.4;0.5;0.6;0.7;0.8;0.9];
C=0.01*c1-1*[00;0;0;0;0;0;0;1];
y=A\C;
x=0:0.1:1;
yn=[0;y;1];
ye=2*(exp(x)-exp(-x))/(exp(1)-exp(-1))-x;
figure(1);
plot(x,yn,'*',x,ye);
legend('numerical solution','exact solution')
xlabel('x','fontsize',20);
ylabel('y','fontsize',20);
set(gca,'fontsize',18);

figure(2);
err=abs(ye'-yn);
plot(x,err);
legend('error')
xlabel('x','fontsize',20);
ylabel('y','fontsize',20);
set(gca,'fontsize',18);

第2题：
%========================Boundary Value Problem 1_Explicit
%显式
clear;clc
l=20;%板厚
h=1;%步长
J=l/h;

T=50;%时间
tao=2.5;%步长
N=T/tao;
%下面组合A矩阵
a=0.2;
lamda=tao/(h^2);
zhu=1-2*a*lamda;
ce=a*lamda;
a00=ones(1,J-1);
a0=diag(a00);
A0=zhu*a0;
a01=ones(1,J-2);
a1=diag(a01,1);
A1=ce*a1;
a2=diag(a01,-1);
A2=ce*a2;
A=A0+A1+A2;
u(:,1)=0; %板的初始温度
for i=2:N+1
u(1,i)=100-100*exp(-(i-1)*tao); %边界条件
u(J+1,i)=100-100*exp(-(i-1)*tao); %边界条件
end
% g01=u(:,1);
% g02=u(:,J);
for k=1:N
% g01=ce*g1(1,k);
% g02=ce*g2(1,k);
oo=zeros(J-3,1);
g(:,k)=[ce*u(1,k);oo;ce*u(J+1,k)];
u(2:end-1,k+1)=A*u(2:end-1,k)+g(:,k);
end
t=0:h:l;
x=0:tao:T;
mesh(x,t,u)
xlabel('t','fontsize',20);
ylabel('x','fontsize',20);
zlabel('T','fontsize',20);
set(gca,'fontsize',18);

%========================Boundary Value Problem 1_2Implicit
%隐式
clear;clc
l=20;%板厚
h=1;%步长
J=l/h;

T=50;%时间
tao=2.5;%步长
N=T/tao;
%下面组合A矩阵
a=0.2;
lamda=tao/(h^2);
zhu=1+2*a*lamda;
ce=-a*lamda;
a00=ones(1,J-1);
a0=diag(a00);
A0=zhu*a0;
a01=ones(1,J-2);
a1diag(a01,1);
A1=ce*a1;
a2=diag(a01,-1);
A2=ce*a2;
A=A0+A1+A2;
u(:,1)=0; %板的初始温度
for i=2:N+1
u(1,i)=100-100*exp(-(i-1)*tao); %边界条件
u(J+1,i)=100-100*exp(-(i-1)*tao); %边界条件
end
% g01=u(:,1);
% g02=u(:,J);
for k=1:N
% g01=ce*g1(1,k);
% g02=ce*g2(1,k);
oo=zeros(J-3,1);
g(:,k)=[ce*u(1,k);oo;ce*u(J+1,k)];
u(2:end-1,k+1)=inv(A)*u(2:end-1,k)-inv(A)*g(:,k);
end
t=0:h:l;
x=0:tao:T;
mesh(x,t,u)
xlabel('t','fontsize',20);
ylabel('x','fontsize',20);
zlabel('T','fontsize',20);
set(gca,'fontsize',18);

第3题：
%=============================used by centered difference
clear;

function pdemodel
[pde_fig,ax]=pdeinit;
pdetool('appl_cb',1);
set(ax,'DataAspectRatio',[1 1 1]);
set(ax,'PlotBoxAspectRatio',[1.5 1 1]);
set(ax,'XLim',[-1.5 1.5]);
set(ax,'YLim',[-1 1]);
set(ax,'XTickMode','auto');
set(ax,'YTickMode','auto');

% Geometry description:
pdecirc(0,0,1,'C1');
set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','C1')

% Boundary conditions:
pdetool('changemode',0)
pdesetbd(4,...
'dir',...
1,...
'1',...
'0')
pdesetbd(3,...
'dir',...
1,...
'1',...
'0')
pdesetbd(2,...
'dir',...
1,...
'1',...
'0')

pdesetbd(1,...
'dir',...
1,...
'1',...
'0')

% Mesh generation:
setappdata(pde_fig,'Hgrad',1.3);
setappdata(pde_fig,'refinemethod','regular');
pdetool('initmesh')
pdetool('refine')

% PDE coefficients:
pdeseteq(1,...
'1.0',...
'0.0',...
'1',...
'1.0',...
'0:10',...
'0.0',...
'0.0',...
'[0 100]')
setappdata(pde_fig,'currparam',...
['1.0';...
'0.0';...
'1 ';...
'1.0'])

% Solve parameters:
setappdata(pde_fig,'solveparam',...
str2mat('0','1524','10','pdeadworst',...
'0.5','longest','0','1E-4','','fixed','Inf'))

% Plotflags and user data strings:
setappdata(pde_fig,'plotflags',[1 1 1 1 1 1 1 1 0 0 0 1 0 0 1 0 0 1]);
setappdata(pde_fig,'colstring','');
setappdata(pde_fig,'arrowstring','');
setappdata(pde_fig,'deformstring','');
setappdata(pde_fig,'heightstring','');

% Solve PDE:
pdetool('solve')