有限差分和Matlabpde求解一维稳态传热问题.(优选)

有限差分和pde 函数求解一维定态热传导方程分别用有限差分方法和pde 函数求解一维定态热传导方程，初始条件和边界条件，热扩散系数α=0.00001，22T T t x α∂∂=∂∂ （1） 求解过程：1. 用Tylaor 展开法推导出FTCS 格式的差分方程首先对T 进行泰勒展开得到如下两式子：23123123...232!3!23...232!3!nnnn n j j jj j nnnn n j jjjjttT T t x x TTx T T T t t t T T T x x x ++∆∆=+∆+++∆∆=+∆+++⎛⎫⎛⎫∂∂∂⎛⎫⎪ ⎪⎪∂∂∂⎝⎭⎝⎭⎝⎭⎛⎫⎛⎫∂∂∂⎛⎫⎪⎪ ⎪∂∂∂⎝⎭⎝⎭⎝⎭上述两个方程变换得：()1122323...23nnn n n n n j j j j j j jT T T T T t T t T o t t t t t t ++--⎛⎫⎛⎫∂∆∂∆∂⎛⎫=--=+∆ ⎪ ⎪ ⎪∂∆∂∂∆⎝⎭⎝⎭⎝⎭（2）223123...23nnn n n j j jj j T T T x T x T x x x x --⎛⎫⎛⎫∂∆∂∆∂⎛⎫=-- ⎪ ⎪ ⎪∂∆∂∂⎝⎭⎝⎭⎝⎭ ()1232422342222...3!4!nnnn n n j j j j j j T T T T x T x T x x x x x x +-⎛⎫⎛⎫⎛⎫∂∂∆∂∆∂⎛⎫=--- ⎪ ⎪ ⎪ ⎪∂∆∆∂∂∂⎝⎭⎝⎭⎝⎭⎝⎭()()21122222-n n n j j j T T T To x x x+--+⎛⎫∂=+∆ ⎪∆∂⎝⎭（3）将上述式子（2）（3）代入（1）得：12112222()nnn nn n n j j j j j j jT T T T T T T O t x t x t x αα+-+--+⎛⎫∂∂⎛⎫-=-+∆∆ ⎪ ⎪∂∂∆∆⎝⎭⎝⎭，（4）2.方程的相容性和稳定性讨论：上述方程截项为：2233242334()...4...23!3!4!n nn n j j j jt T t T x T x T O t x t t x x α⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫∆∂∆∂∆∂∆∂ ⎪∆∆=--++ ⎪ ⎪ ⎪ ⎪ ⎪∂∂∂∂⎝⎭⎝⎭⎝⎭⎝⎭⎝⎭， 由于(),0,0lim x t o t x ∆∆→∆∆=所以方程有相容性其经过傅里叶变换后，只有当,t x ∆∆满足下列条件时，方程具有较好的稳定性：2220sin 12m k xt x α∆∆≤≤∆ 其中m Nk Lπ=N 为节点个数，L 为边界长度由于：22sin sin 122m k x N x L π∆∆⎛⎫== ⎪⎝⎭所以当20.5ts x α∆=≤∆时方程具有稳定性3. 说明该方程的类型和定解条件，如何在程序中实现这些定解条件。

由式（4）得：1111122(12)n n n j j j n n n n n jjj j j T T T TT tsT s T sT x α-++-+-+≈+∆=+-+∆要想知道下一个时刻1n t +下的温度，则必须知道在n t 时刻下该点以及相邻两点的温度，如下图所示：所以需要知道的是在各个时刻x 轴两侧的温度，以及初始时刻x 轴上各点温度。

下述各题求解过程中，均采用0max 0.000010(1,2...1)j T j j α===-，0100n T =，作为边界条件。

边界条件是由外界给定，本身对过程求解，精度无影响4. 编写M 文件求解上述方程 M 文件的编写见附录所示。

在matlab 中调用heat_conduct(0.00001,15000,500,0.1)，得到如下结果： rms error is rms = 0.7959∆∆分别为0.1和500，时间最终给定为15000 说明：1.上述结果采用的,x t当(),x t∆∆为（0.01,5）时的结果附录1，传热方程有限差分的M文件编写function heat_conduct(alpha,Timemax,dt,dx)%方程需要输入的参数有，最大时刻数，alpha，dt，dx,实际Tmax=Timemax/dt+1;%时刻数JMAX=1/dx+1;%x轴上的差分点个数J=JMAX-1;MAXEX=J;T=zeros(Tmax,JMAX);%构造温度矩阵TE=zeros(Tmax,JMAX);%实际温度矩阵%设定边界条件T(:,1)=100;T(:,JMAX)=100;s=alpha*dt/dx/dx;%计算各个点的温度for i=1:(Tmax-1)for j=2:JT(i+1,j)=s*T(i,j+1)+(1-2*s)*T(i,j)+s*T(i,j-1);endend%计算实际值%画图z=0:dx:1;t=0:(Tmax-1);a=3000/dt+1;b=9000/dt+1;subplot(1,2,1)plot(z,T(a,:),'g*',z,T(a,:),'g',z,T(b,:),'o',z,T(b,:),'r',z,T(Tmax,:),'b')xlabel('x方向')ylabel('温度T')title('3000,9000,15000时刻的温度分布')%画整体分布图for i=1:JMAXX(:,i)=t; end; for j=1:Tmax Y(j,:)=z; endsubplot(1,2,2) mesh(X,Y,T); view([-1 1 1]); xlabel('时刻t'); ylabel('方向x'); zlabel('温度T'); title('整体分布')二、调用MATLAB 函数完成上述计算1. 编写M 文件求解上述方程，并用适当的文字对程序做出说明。

Matlab 中自带的pdepe 函数可直接求解一般偏微分方程。

其调用格式为： T=pdepe(m,@pdefun,@pdeic,@pdebc,x,t);其中T 为输出结果，在本题中，其输出格式为一个二维数组，分别表示一维空间和时间上的温度；x ，t ，分别为所需输入的时间和空间数组；为达到该pdepe 函数求解要求，将上述方程转化为如下函数：即@pdefun 所要求的函数的描写格式：,,,,,,,,m mT T T c x t x x f x t u s x t u x t x x x -∂∂∂⎡∂⎤∂⎛⎫⎛⎫⎛⎫=+ ⎪ ⎪ ⎪⎢⎥∂∂∂∂∂⎝⎭⎝⎭⎝⎭⎣⎦本题中，在描述@pdefun 时，,,1T c x t x ∂⎛⎫= ⎪∂⎝⎭；10mx m -=⇒=；,,,T T f x t u x x α∂∂⎛⎫= ⎪∂∂⎝⎭；,,,0s x t u x ∂⎛⎫= ⎪∂⎝⎭@pdeic 为初值条件，其描述的标准格式为：()00,0T x =@pdebc 为边界条件，本题要求边界温度为100℃，其描述的标注格式为：()(),,,,*,,,0T p x t T q x t T f x t u x ∂⎛⎫+= ⎪∂⎝⎭所以左右边界的(),,p x t T 与(),,q x t T 分别表示为：(),,100p x t T T =-边界边界；()q x t T,,=0边界上述程序编写见附录2heat_conduct1(0.1,500,15000)得到如下结果所示：附录2:Matlab PDE求解上述一维传热方程%带求解函数function [c,f,s]=pdefun(x,t,T,dT)c=1;f=0.00001*dT;s=0;%边值条件，a为左边界，b为右边界，都为100℃function [pa,qa,pb,qb]=pdebc(xa,Ta,xb,Tb,t)pa=100-Ta;qa=0;pb=100-Tb;qb=0;%初值设定为0function T0=pdeic(x)T0=0;调用求解函数function heat_conduct1(dx,dt,Timemax)clcx=0:dx:1;t=0:dt:Timemax;m=0;T=pdepe(m,@pdefun,@pdeic,@pdebc,x,t);figure('numbertitle','off','name','PDE Demo——by Matlabsky')subplot(1,2,2)surf(x,t,T(:,:))title('solution')xlabel('X')ylabel('t')zlabel('T')subplot(1,2,1)plot(x,T(7,:),'b',x,T(19,:),'g*',x,T(19,:),'g',x,T(31,:),'ro',x,T(31,:),'r')title('3000,9000,15000时刻温度')xlabel('X')ylabel('T')end最新文件---------------- 仅供参考--------------------已改成-----------word文本--------------------- 方便更改赠人玫瑰，手留余香。