偏微分一维热传导问题

偏微分大作业一维热传导方程问题——运用隐式格式求解数值解目录问题描述.....................................................................................................31解析解——分离变量法.........................................................................32数值解——隐式格式.............................................................................53证明隐式格式的相容性与稳定性.........................................................54数值解——分析与Matlab实现...........................................................65数值解与解析解的比较.........................................................................96随时间变化的细杆上的温度分布情况...............................................117稳定后细杆上的温度分布情况............................................................12参考文献...................................................................................................13附录...........................................................................................................14有限长杆的一维热传导问题问题描述一根单位长度的细杆放入100℃的沸水中，当细杆的温度达到100℃时取出。假设细杆四周绝热；在时间t=0时，细杆两端浸入0℃的冰水中。一维热传导方程：20txxuau，现在令21a，从而可知本题：0txxuu。现在要求细杆温度分布：(,)uxt。1解析解——分离变量法热传导偏微分方程：0txxuu(1)(0,t)(1,t)0uu(0)()uxx，其中，001xx，或()x100(0,1)x，首先令：(,)()()uxtXxTt(2)将(2)式带入(1)式得：()T()()()0XxtTtXx于是可得：T()()()()tXxTtXx可以得到两个微分方程：T()()0tTt()()0XxXx先求解空间项：当0时，()xxXxAeBe由于(0,t)(1,t)0,.uut可知:由于解的收敛性，0B(0)=(1)00XXAAeA则此时是平庸解。当0时，()XxABx(0)=(1)00,0XXAABAB则此时是平庸解。当0时，()cossinXxAkxBkx，其中k。(0)00XAA(1)sin0,1,2,3...XBkknn所以，()sin()nXxBnx，1,2,3...n因为22n所以，22()ntnTtCe，1,2,3...n则，221(,)sin()ntnnuxtDenx初始条件：(0)()uxx，1(0)sin()()nnuxDnxx，102()sin()nDxnxdx12100sin()1200()cos(1)cosnxdxnnn2000lim=(1cos)nDnn当时，最终，221200(,)(1(1))sin()nntnuxtenxn，1,2,3...n2数值解——隐式格式目前，研究热传导问题特别是非稳态热传导问题十分重要。这里使用隐式格式1。利用(,)uxt，关于t进行向前差商：1kkjjUUt；关于x进行二阶中心差：1111122()kkkjjjUUUx；代入偏微分方程可以得到隐式差分格式：11111122()kkkkkjjjjjUUUUUtx（1）3证明隐式格式的相容性与稳定性（1）相容性2122222122331222212233122Taylor1=+()211+()()2211+()()22kkjjkkjjkkjjUUUUtttttUUUUUUttxxtxttxxUUUUUUttxxtxttxx根据展开：代入隐式格式得：2222221()=+()()2UUUtttxttx(2)将（2）与原微分方程相减，得到截断误差：221=-()02Utxt所以此隐式格式与原微分方程相容。（2）稳定性令网格比为2rtx，则可以将（1）式改写得到：11111(12)kkkkjjjjrUrUrUU（3）首先令：11(-111111(+1+1kkIjjkkIjjkkIjjkkIjjUUeUUeUUeUUe））（4）将（4）代入（3）式，根据欧拉公式化简得：11+22cos)kkUrrU（（5）故得放大因子是：1111+2(1cos)kkUGUr所以根据Fourier方法，隐式格式恒稳定。4数值解——分析与Matlab实现（1）边值与初值离散化将边值与初值离散化，与式（3）联立得差分线性方程组：11111(12)kkkkjjjjrUrUrUU，(0,1,2,,M-1)j(0,1,2,,N-1)k0=(),jjUx(0,1,2,,M)j00,kU(0,1,2,,N)k0kLU，(0,1,2,,N)k再将方程组改写成AUB的形式：111101221331221111(1)(M1)1+2+12121212kkkkkkkkkMMkkkMMMMrrUUrUrrrUUrrUUrrrUUrrUUrU本题的边界条件均为零。所以可以将上式改写。111122133122111(1)(M1)1+212121212kkkkkkkkMMkkMMMrrUUrrrUUrrUUrrrUUrrUU（2）Matlab的实现杆长1米，时间2秒。设计空间步长h=0.1和时间步长t=0.01，网格比是2trh。从而得到划分的空间网格点数是M1+1，时间网格点数是M2+1。先设初始的温度矩阵U(M2+1,M1+1)。再将边界条件和初始条件编写到表示温度分布的矩阵中。具体代码可见最后附录。编写矩阵A核心代码：对角线：A(i,i)=1+2r对角线的右方和下方：A(i,i+1)=-r;A(i+1,i)=-r;下面就要运用*(1,)(,)AUkjUkj进行迭代。当k=1时，A*U(2,j)=U(1,j)当k=2时，A*U(3,j)=U(2,j)当k=3时，A*U(4,j)=U(3,j)以此迭代下去直到k=M2。就可以得到整个温度随时间和空间的分布矩阵U。数值解画图，如图1(a)和图1(b)所示。图1(a)数值解的温度分布图现在将着色平稳过渡。图1(b)着色平稳过渡的数值解的温度分布图5数值解与解析解的比较首先，我们需要将解析解离散化，解析解中有一项22nte，当n越来越大时，会快速趋于0，故我们可以取n=8000。现在来证明可行性，在matlab里的工作空间运算。将解析解的温度分布画出来，数值解画图2，如图2所示。图2解析解的温度分布图将数值解与解析解相减，得到误差图。如图3(a)和图3(b)，我们从图3(a)上可以看出空间上的误差，在边界处误差比较大。图3（a）数值解与解析解空间误差我们从图3(a)上可以看出时间的误差，在时间的最开始，处误差最大，然后又有一个小的波动，最后就误差渐渐变小，最后趋于0。图3（b）数值解与解析解时间误差6随时间变化的细杆上的温度分布情况从数值解的温度分布三维图，如图4(a)和图4(b)可以看出随着时间的增加，细杆温度下降最后趋于0℃。从物理角度来说：细杆的温度会不断地向两端扩散，热量会慢慢散失，最终随着时间的增加，细杆的温度会趋于0℃。图4(a)细杆温度随时间的变化图现取细杆中心处一点，观看它随时间的温度变化情况。图4(b)细杆中央（x=0.5）温度随时间的变化图7稳定后细杆上的温度分布情况从图像上可以看出，最后稳定的情况下，细杆的温度是0℃。参考文献[1]冯立伟．热传导方程几种差分格式的MATLAB数值解法的比较[J]．沈阳化工大学，辽宁沈阳．2011(6)．[2]一维热传导方程数值解法及Matlab实现[EB/OL]．2014-11-20附录代码：%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%此程序用于解决一维热传导方程:ut-a^2uxx=0%%边界条件：u(0,t)=u(L,t)=0%%初始条件：u(x,0)=100,x!=0和L%%u(0,0)=0%%u(L,0)=0%%其中,a^2=1,L=1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clc;clearall;%区域及划分网格L=1;%单位长度的细杆T=2;%时间h=0.1;%%%%空间的划分%%%%t=0.01;%%%%时间的划分%%%%r=t/(h*h);%网格比%设计步长M1=L/h;M2=T/t;%构造边界条件%构造的矩阵：U(时间，空间）U=zeros(M2+1,M1+1);%编程包含边值，如U(k,1)=u(0,t)fork=1:M2+1%时间划分了M2份，有M2+1个节点U(k,1)=0;%两个边界处温度恒为零U(k,M1+1)=0;end;%构造初始条件forj=2:M1%位置划分了M1份，有M1+1个节点U(1,j)=100;end;U(1,1)=0;U(1,M1+1)=0;%差分格式的矩阵形式A*U(k+1,j)=U(k,j)%构造矩阵AA=zeros(M1-1);fori=1:M1-1A(i,i)=1+2*r;end;fori=1:M1-2A(i,i+1)=-r;A(i+1,i)=-r;end;%构造AU=B中的B%本题边值的特殊，矩阵B大大简化了B=zeros(M1-1,1);fork=1:M2j=2:M1;B(j-1,1)=U(k,j);x=A\B;forj=2:M1U(k+1,j)=x(j-1);%k+1时刻的不同位置的温度分布end;end;%作图x=0:h:1;y=0:t:2;[xx,yy]=meshgrid(x,y);figure(1);surf(xx,yy,U);shadingflattitle('一维热传导方程--数值解--温度分布图');xlabel('位置x');ylabel('时间t');zlabel('温度T');figure(2)s=0;fori=1:8000s=s+(200*(1-(-1)^i))/(i*pi)*sin(i*pi*xx).*exp(-i^2*pi^2*yy);end;surf(xx,yy,s);t