热传导方程有限差分法的MATLAB实现

《热传导方程有限差分法的MATLAB实现》由会员上传分享，免费在线阅读，更多相关内容在应用文档-天天文库。

1、万方数据2009年7月第24卷第4期咸阳师范学院学报JournalofXianyaagNormalUniversityJul．2009V01．24No．4[理论物理与应用物理学研究]热传导方程有限差分法的MATLAB实现史策(西安建筑科技大学理学院，陕西西安710055)摘要：对于有界热传导齐次方程的混合问题．用分离变量法求解往往很复杂。为了更好地理解热传导方程的解，使用MATLAB软件将方程的解用图像表示出来。通过区域转换的思想，利用MATLAB编程实现一定区域内热传导方程的有限差分方法，数值

2、表明了方法的可行性和稳定性。关键词：热传导方程；有限差分；MATLAB中图分类号：0552文献标识码：A文章编号：1672—2914(2009)04—0027-03近些年来。求解热传导方程的数值方法【ll取得进展，特别是有限差分区域分解算法121，此类算法的特点是在内边界处设计不同于整体的格式．将全局的隐式计算化为局部的分段隐式计算。使人从感觉上认为这样得到的解会比全局隐式得到的解的精度差，但大量的数值实验表明事实正好相反。用区域分解算法求得的解的精度更好。MATLAB具有强大的图形绘制功能【3

3、1．为科学计算和图形处理提供了很大的方便。用户只须指定绘图方式，并提供充足的绘图数据，用很少的程序指令就可得到直观、形象的图形结果。因此。近些年来。越来越多的人开始使用MATLAB来求解热传导方程14’5】。借助MATLAB的数值计算和图形处理技术[61。我们可以绘制出热传导方程数值解的二维、三维图形．从而可以更好地理解热传导方程解的意义。一维热传导方程婴=舻霎譬，Otd髫。是最简单的偏微分方程之一。其定解问题的数值解法主要有有限元法和有限差分法等。对于有限元法来说，适用处理复杂区域、精度可选；

4、缺点在于内存和计算量巨大。不易编程实现。对于有限差分法来说，虽然比较直观、理论也比较成熟、精度可选；但是不规则区域处理繁琐．网格生成可以使有限差分方法[71(FDM)应用于不规则区域。但是对区域的连续性等要求较严。适用FDM的好处在易于编程，易于并行。鉴于以上情况．本文考虑以下边界值问题：等彳詈问舢tt[矿O，uI，产0，t>0u驴sin(平)，o<矧利用区域转化的思想通过极坐标对求解区域进行转化，求解区域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以Taylor级数展开等方法．

5、把控制方程中的导数用网格节点上的函数值的差商代替进行离散。从而建立以网格节点上的值为未知数的代数方程组。1求解热传导方程的基本思想基本思想是把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似：把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题

6、在整个区域上的近似解【8期。下面是有限差分法数值计算的基本步骤：(1)区域的离散或子区域的划分；(2)插值函数的选择；(3)方程组的建立；(4)方程组的求解。收稿日期：2009—04—20作者简介：吏策(1986-)，男。陕西兴平市人，西安建筑科技大学理学院硕士研究生，研究方向为微分方程数值解法。万方数据·28·咸阳师范学院学报第24卷2热传导方程的离散分析r庸+lnn一1AxlAt(_，，曲T(j-1一‘＼j．／’(j+l，，1)(j,n-1)图1热传导方程隐格式网格划分通过已知方程，建立一个

7、关于时间和步长(戈一t坐标)的函数，设步长为Ax=I／M，每次运行的时间为At=T／N．这样就把初始的矩形区域划分成了一个长方形的网格．本文就是通过对原方程建立的差分格式．以及对初始条件和边界条件建立相应的差分近似进行计算，即把原方程离散到各个节点上从而进行数值近似解的计算。利用Ⅱ以‘)关于f的向后差商㈣：等角兰著L，关于￡的二阶中心差商[101．．雾j兰蔓悉#纽，对方程进行离散。离散后的方程为：npln^nn．孚--d气誊堑。令：产器Ax)，即&(缸)2。一‘(2’”7(1+2r)#-r·吱。

8、一r·吩n一。=譬一。可化为矩阵形式：l+2r_rO-rl+2r_r0O—r0⋯0，卜lnⅡI+r。‰pl％a-It‰■0000-rl+2rL㈣⋯I●l+2r．rL．．，⋯．．。、3热传导方程的MATLAB编程％热传导方程问题的数值解clear；clc；formatshortea=input(’请输入系数a的值：’)；l=input(’请输入长度l的值：’)；M=input(’请输入将Af．1[0，l】等分的个数M：’)；ot=input(’请输入时间增量ot的值：’)；n=input(’请输人