有限差分法

有限差分法 一、单变量函数： 用中心差分法（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 隐式差分法 小结：显式 格式 pdf格式笔记格式下载页码格式下载公文格式下载简报格式下载 仅当时格式是稳定的。（其中称为网格比） 隐式格式从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')