热传导方程几种差分格式的MATLAB数值解法比较

《热传导方程几种差分格式的MATLAB数值解法比较》由会员上传分享，免费在线阅读，更多相关内容在应用文档-天天文库。

1、第25卷第2期沈阳化工大学学报Vol．25No．22011．06JOURNALOFSHENYANGUNIVERSITYOFCHEMICALTECHNOLOGYJun．2011文章编号:2095－2198(2011)02－0179－04热传导方程几种差分格式的MATLAB数值解法比较冯立伟(沈阳化工大学数理系，辽宁沈阳110142)摘要:对求解热传导方程的几种差分格式进行分析，并讨论使用MATLAB编程求解偏微分方程的方法．编制几种差分格式的MATLAB程序，使用算例进行数值实验，在不同网格比情况下，比较几种算法的优劣．

2、关键词:热传导方程;MATLAB;Crank-Nicolson离散中图分类号:O241.82文献标识码:A许多工程问题需要研究热量在物体内部的为空间步长与时间步长，用2族平行直线xj=传导情况或某种物质在液体中的扩散情况，因此jh，tk=kτ将矩形区域［0，T］×［0，L］分割成矩研究热传导问题特别是非稳态热传导问题十分形网格．［1－2］重要．目前热传导方程已有多种求解格式．显式格式:MATLAB是目前最流行、应用最广泛的科学和工k+1kUj－Uj=程计算软件．MATLAB基于矩阵运算，具有强大τkkk的数值运算能力和

3、图形可视化能力，是方便实Uj+1－2Uj+Uj－1k［3］a2+fJ．(2)用、功能强大的数学软件．用MATLAB求解常h［4－6］微分方程已有大量的研究．王飞等介绍了如隐式格式:k+1k何使用MATLAB实现有限差分法求解微分方U－Ujj［7］=程．高理平等给出了对两点边值问题有限元τ［8］．李灿等对热传导问题的MAT-k+1k+1k+1方法的程序Uj+1－2Uj+Uj－1k［9］a2+fJ．(3)LAB数值计算进行了讨论．本文讨论求解一h维热传导方程几种不同差分格式的MATLAB编Crank-Nicolson格式:

4、程方法，并使用算例进行检验和对结果进行分Uk+1－Ukk+1k+1k+1jjUj+1－2Uj+Uj－1=a+析．τ(2h2uuUk－2Uk+Uk=a+f(x，t)(1)j+1jj－1kt2)+fj．(4)x2hu(x，t)表示在t时刻物体内部坐标为x处的温DuFortFrankel格式:度，a是热传导系数，f(x，t)为热源．k+1k－1Uj－Uj=2τ1热传导方程差分格式kk+1k－1kUj+1－Uj－Uj+Uj－1ka2+fj．(5)h1.1差分格式的建立由Taylor公式容易得出:它们都与一维热传导2首先

5、对x-t平面进行网格剖分．分别取h，τ方程相容，其截断误差分别为O(τ+h)，O(τ+收稿日期:2010－08－30作者简介:冯立伟(1980－)，男，辽宁沈阳人，助教，硕士，主要从事偏微分方程数值解的研究．180沈阳化工大学学报2011年22222［1］1h)，O(τ+h)和O(τ+h)．可算出u在第一层各个节点处的近似值uj．重k复使用此式，可以逐层计算出所有的uj．1.2初、边值条件的处理隐式格式:将(7)式与离散化的初边值条件对定解条件进行离散化．由初始条件及第一联立，得差分方程组:k+1k+1k+1kk类边界

6、条件，可直接得到:ì－rUj－1+(1+2r)Uj－rUj+1=Uj+τfj0ïUj=u(xj，0)=j，(j=0，1，…，M);ï(k=1，2，…，N－1，j=0，1，…，M－1)Uk=u(0，t)=g，í0(11)0k1kïUj=j(k=0，1，…，N)kïUL=u(L，tk)=g2k，(k=1，…，N)．îUk=g，Uk=g(j=0，1，…，M)01kL2k将上述方程组改写成矩阵形式1.3稳定性分析aτ令r=2为网格比，显式格式公式(2)变为:hk+1kkUJ=rUj+1+(1－2r)Uj+kkrUj－1+τ

7、fj(6)使用Fourier方法可知，当r≤1/2时显式格k+1k式稳定．éU1ùéU1+rg1j+1ùêUk+1úêUkú隐式格式变为:22êúêú－rUk+1+(1+2r)Uk+1－êú=êú(12)j－1jêk+1úêkúk+1kkUUrUj+1=Uj+τfj．(7)êM－2úêM－2úêúêú由Fourier方法可得隐式格式恒稳定．ëUk+1ûëUk+rgûM－1M－12j+1Crank-Nicolson格式变为:此方程组是三对角方程组，且系数矩阵严格对角rk+1k+1占优，故解存在唯一．通过在每一时间层上

8、求解－Uj－1+(1+r)Uj－2一个这样的线性方程组得到在各个时刻各个网rUk+1=rUk+格点上的函数值．j+1j+122Crank-Nicolson格式:将(8)式与离散化的krkk(1－r)Uj+Uj－1+τfj．(8)初边值条件联立并整理，得差分方程组:2ì－rUk+1+(1+r)Uk+1－rUk+1=对于r＞0恒有增长