有限差分法一、单变量函数：用中心差分法（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 1clear;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;clcl=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+1u(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);endt=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;clcl=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+1u(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);endt=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')。