有限差分和Matlabpde求解一维稳态传热问题

有限差分和pde函数求解一维定态热传导方程分别用有限差分方法和pde函数求解一维定态热传导方程，初始条件和边界条件，热扩散系数α=0.00001，22TTtx（1）求解过程：1.用Tylaor展开法推导出FTCS格式的差分方程首先对T进行泰勒展开得到如下两式子：23123123...232!3!23...232!3!nnnnnjjjjjnnnnnjjjjjttTTtxxTTxTTTtttTTTxxx上述两个方程变换得：1122323...23nnnnnnnjjjjjjjTTTTTtTtTotttttt（2）223123...23nnnnnjjjjjTTTxTxTxxxx1232422342222...3!4!nnnnnnjjjjjjTTTTxTxTxxxxxx21122222-nnnjjjTTTToxxx（3）将上述式子（2）（3）代入（1）得：12112222()nnnnnnnjjjjjjjTTTTTTTOtxtxtx，（4）2.方程的相容性和稳定性讨论：上述方程截项为：2233242334()...4...23!3!4!nnnnjjjjtTtTxTxTOtxttxx，由于,0,0limxtotx所以方程有相容性其经过傅里叶变换后，只有当,tx满足下列条件时，方程具有较好的稳定性：2220sin12mkxtx其中mNkLN为节点个数，L为边界长度由于：22sinsin122mkxNxL所以当20.5tsx时方程具有稳定性3.说明该方程的类型和定解条件，如何在程序中实现这些定解条件。由式（4）得：1111122(12)nnnjjjnnnnnjjjjjTTTTTtsTsTsTx要想知道下一个时刻1nt下的温度，则必须知道在nt时刻下该点以及相邻两点的温度，如下图所示：所以需要知道的是在各个时刻x轴两侧的温度，以及初始时刻x轴上各点温度。下述各题求解过程中，均采用0max0.000010(1,2...1)jTjj，0100nT，作为边界条件。边界条件是由外界给定，本身对过程求解，精度无影响4.编写M文件求解上述方程M文件的编写见附录所示。在matlab中调用heat_conduct(0.00001,15000,500,0.1)，得到如下结果：rmserrorisrms=0.7959说明：1.上述结果采用的,xt分别为0.1和500，时间最终给定为15000当,xt为（0.01,5）时的结果附录1，传热方程有限差分的M文件编写functionheat_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;%计算各个点的温度fori=1:(Tmax-1)forj=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时刻的温度分布')%画整体分布图fori=1:JMAXX(:,i)=t;end;forj=1:TmaxY(j,:)=z;endsubplot(1,2,2)mesh(X,Y,T);view([-111]);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所要求的函数的描写格式：,,,,,,,,mmTTTcxtxxfxtusxtuxtxxx本题中，在描述@pdefun时，,,1Tcxtx；10mxm；,,,TTfxtuxx；,,,0sxtux@pdeic为初值条件，其描述的标准格式为：00,0Tx@pdebc为边界条件，本题要求边界温度为100℃，其描述的标注格式为：,,,,*,,,0TpxtTqxtTfxtux所以左右边界的,,pxtT与,,qxtT分别表示为：,,100pxtTT边界边界；,,=0qxtT边界上述程序编写见附录2heat_conduct1(0.1,500,15000)得到如下结果所示：附录2:MatlabPDE求解上述一维传热方程%带求解函数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;%初值设定为0functionT0=pdeic(x)T0=0;调用求解函数functionheat_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','PDEDemo——byMatlabsky')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