有限差分法

有限差分法

一、单变量函数:

有限差分法

用中心差分法(matlab程序见附录)计算结果如下:

有限差分法

图1 中心差分法

有限差分法

表1 数据对比

有限差分法

二、一维热传导:

有限差分法

在此取φ(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*[0;0;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);

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)=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')

相关文档