【毕业设计(论文)】二维热传导方程有限差分法的matlab实现管理资料

第1章前言在史策教授的《一维热传导方程有限差分法的MATLAB实现》和曹刚教授的《一维偏微分方程的基本解》中，对偏微分方程的解得MATLAB实现问题进行过研究，但只停留在一维中，而实际中二维和三维的应用更加广泛。诸如粒子扩散或神经细胞的动作电位。也可以作为某些金融现象的模型，诸如布莱克-斯科尔斯模型与Ornstein-uhlenbeck过程。热方程及其非线性的推广形式也被应用与影响分析。在科学和技术发展过程中，科学的理论和科学的实验一直是两种重要的科学方法和手段。虽然这两种科学方法都有十分重要的作用，但是一些研究对象往往由于他们的特性（例如太大或太小，太快或太慢）不能精确的用理论描述或用实验手段来实现。自从计算机出现和发展以来，模拟那些不容易观察到的现象，得到实际应用所需要的数值结果，解释各种现象的规律和基本性质。科学计算在各门自然科学和技术科学与工程科学中其越来越大的作用，在很多重要领域中成为不可缺少的重要工具。而科学与工程计算中最重要的内容就是求解科学研究和工程技术中出现的各种各样的偏微分方程或方程组。解偏微分方程已经成为科学与工程计算的核心内容，包括一些大型的计算和很多已经成为常规的计算。为什么它在当代能发挥这样大的作用呢？第一是计算机本身有了很大的发展；第二是数值求解方程的计算法有了很大的发展，这两者对人们计算能力的发展都是十分重要的。近三十年来，解偏微分方程的理论和方法有了很大的发展，而且在各个学科技术的领域中应用也愈来愈广泛，在我国，偏微分方程数值解法作为一门课程，不但在计算数学专业，而且也在其他理工科专业的研究生的大学生中开设。同时，求解热传导方程的数值算法也取得巨大进展，特别是有限差分法方面，此算法的特点是在内边界处设计不同于整体的格式，将全局的隐式计算化为局部的分段隐式计算。而且精度上更好。目前，在欧美各国MATLAB的使用十分普及。在大学的数学、工程和科学系科，MATLAB被用作许多课程的辅助教学手段，MATLAB也成为大学生们必不可少的计算工具，甚至是一项必须掌握的基本技能。在我国，MATLAB在各大专院校的应用日益普遍，许多专业已把MATLAB作为基本计算工具。在科研机构和工业界，MATLAB正得到越来越广泛的应用。MATLAB具有强大的图形绘制功能，为科学计算和图形处理提供了很大的方便。我们只需制定的绘图方式，再提供绘图数据，有程序指令就可以得到形象、直观的图形结果。因此，近些年越来越多的人开始使用MATLAB来求解数值计算和图形处理技术，我们也可以绘制出热传导方程数值解的二维、三维图形，从而可以更好的理解热传导方程的意义。问题解决目前，对于求解偏微分方程有很多方法，但差分法和有限元离散法式主要解决问题的两种方法。一般来说，用差分法来接偏微分方程，解得得结果就是方程的准确解函数再借点上的近似值。而用变分近似的方法求解，是将近似解表示成有限维子空间中基函数的线性组合。有限元法也是基于变分原理，由于选择了特殊的基函数，使它能适用于一般的区域。这种基函数是与区域的剖分有关的，近似解表示为基函数的线性组合，二线性组合中的系数，又是剖分节点上或其导数的近似值。有关一维热传导方程的有限差分法求解的MATLAB实现，西安建筑科技大学的史策教授已经解决，本文借鉴史老师的求解思想，对二维热传导方程进行转换，再对解法编程实现，从而进一步对热传导方程进行探讨。二维热传导方程求解在现实生活中的应用也更加广泛，所以有很好的现实意义。第2章预备知识[8]含有未知函数的偏导数的方程称为偏微分方程。[8]方程 称热传导方程（或扩散方程）。其中，是固体的传热过程中在处、时刻的温度。系数称为热传导系数，当时，方程为其中，为维数。[8]在特定条件下求解方程的解。这样的条件成为定解条件。给出了方程和定结条件，就构成了定解问题。[1]一般说，边界条件有下列形式其中为边界的外法向导数。有如下几种特殊形式(1)Dirichlet（或第一类）条件：即值给定；(2)Neumann（或第二类）条件：.即的外法向导数给定；(3)Robbins（或第三类）条件：。定[8]只有出事条件而没有边界条件的定解问题。定[8]只有边界条件而没有初值条件的定解问题。定[8]既有边值条件又有初值条件的定解问题。定[8]定义在上的函数的一个关系式，设，有关系式以上变换称为Fourier变换。其中是虚数单位。[8]由第个时间层推进到第个时间层时差分方程提供了逐点直接计算的表达式，我们称次差分方程为显式格式。[8]有限差分格式在新的时间层上包含有多于一个的节点，这种有限差分格式称为隐式格式。[11]称为向前差分。[11]称为向后差分。[11]称为中心差分。[11]用微分方程的解代替差分方程的全部近似解，这样得到的方程两边的差就是截断误差。[8]给定一个适定的线性初值问题以及与其相容的差分格式，则差分格式的稳定性是差分格式收敛性的充要条件。第3章求解二维热传导方程的基本思想基本思想是把连续的定解区域用有限个离散点构成网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数的近似；把原方程和定解条件中的微商用差商来近似，积分用积分来近似，于是原微分方程和定解条件近似的代之以代数方程，即优先差分方程组，解此方程组就可以得到原问题在离散点上的近似解。下面是有限差分法数值计算的基本步骤：用有限差分方法求解偏微分方程问题必须把连续问题进行离散化。为此首先要对求解区域给出网格剖分，由于求解的问题不同，因此求解区域也不尽相同。下面用例子来说明不同区域的剖分离散。并引入一些常用术语。双曲型和抛物型方程的初值问题，求解区域是我们在的上半平面画出两族平行于坐标轴的直线，把上班平面分成矩形网格。其交点称为节点（或网格点）。可设距离，称其为空间步长，平行线的距离按具体问题而定。可设距离，称其为时间步长。这样两族网格线可以写作网格节点有时记为。双曲型和抛物型方程的初边值问题，设求解区域是这个区域的网格由平行于轴的直线族与平行于轴的直线族所构成，其中选择不同的插值函数对偏微分方程进行估计，可得到不同的差分方程，进而稳定性和精度会有所不同。用Taylor级数展开方法是最常用的方法，下面建立差分格式的同时引入一些基本概念及术语。我们主要从对流方程的初值问题()和扩散方程的初值问题()（其中）进行讨论。假定偏微分方程初值问题的解是充分官话的，由Taylor级数展开有()其中或用，表示看括号内的函数在节点处取的值。利用（）表达式中的第1式和第3式有.如果是满足偏微分方程（）的光滑解，则由此看一看出，偏微分方程（）在处可以近似的用下面的方程来代替()其中为的近似值。（）式称为逼近微分方程（）的有限差分方程或简称差分方程。差分方程再加上初始条件的离散型式就可以按时间逐层推进，算出各层的值。差分格式（）和初始条件的离散形式结合在一起构成了一个差分格式。将离散后的差分方程转化为方程组的形式，便于求解。利用矩阵的解法求解方程组，再用MATLAB对矩阵求解方法进行程序化，以便对以后类似的方程进行求解。隐式差分格式方程矩阵化后，得到的矩阵是严格的对角占优三对角矩阵，我们可以根据线性方程组的求解方法对其求解。其中这要应用的是追赶法，追赶法对于此类线性方程组的求解非常方便，用MATLAB对追赶法进行编程，就可以轻松实现矩阵的求解，进而解出差分方程的近似解。第4章二维热传导方程网格剖分在区域中，我们设二维热传导方程的初始值和边界条件如下：其中为正常数。通过已知方程，建立一个关于时间和步长的函数，这样就把初始区域划分为一个网格图。先将定义域剖分为网格其中为时间步长，分别为轴和轴的空间步长。稳定性分析利用有限差分格式进行计算时是按时间层逐层推进的。那么计算第层上的值时要用到第层上计算出来的结果值，而计算第层结果值时的舍入误差必然会影响到第层的值。从而就要分析这种误差传播的情况。希望误差不至于越来越大，以至掩盖差分格式的解的面貌，这便是稳定性问题。我们先考虑一维差分格式()的稳定性，其中为网格比，假设差分格式从初层开始计算，当初始数据存在误差时考察这个误差在以后计算中的在传播情况。为方便起见，不考虑计算过程中的舍入误差。及确定初始数据误差绝对值为，则差分格式在处的误差为于是，对于固定网格比及的情况，差分格式的解的误差随时间步长的步数的增加而增加。初始数据的误差将必定掩盖了差分格式的解的面貌，所以我们认为差分格式（）时不稳定的。差分格式的稳定性不仅与差分格式本身有关，而且还与网格比的大小有关。差分格式的稳定性在差分方法的研究中具有特别的意义，我们再做进一步的叙述。[8]为了度量误差及其他应用，引入范数设有一个误差，则就有误差。如果存在一个正常数,使得当时，一致的有则称差分格式是稳定的。差分格式一旦具有稳定性，就可以用差分格式计算出偏微分方程的近似解来。一维热传导方程的各种类型的差分格式可以推广到二维热传导方程，利用向前差分格式对（）式进行离散，引入记号（）其中为差分方程在节点的计算值。差分格式（）利用Taylor级数展开易得差分格式（）的截断误差为。为方便稳定性的判断，设，令，为网格比。即改写为：（）用Fourier方法来分析（）式的稳定性。令把此式带入（）式中有因此差分格式（）的增长因子是其中。如果，则有由此得出差分格式（）的稳定性条件是容易看出在二维情况下采用这样的显式格式是不适合的，为此我们再转向考虑隐式格式。用一维向后差分格式的直接推广是（）此格式的截断误差仍为，仍用Fourier方法来分析这个格式的稳定性，仿前可以得出其增长因子是：因此有，即差分格式（）是绝对稳定的。为了提高精度，对微分方程（）也可以用Crank-Nicolson型差分格式，这也是一维问题的直接推广。其格式可写为（）这也是二阶精度格式，其增长因子是因此，对任何都有，所以()式也是绝对稳定的。现在考虑一下隐式格式（）式和（）式的求解方法。我们知道，在一位格式形成的方程组是系数矩阵为三对角矩阵的线性代数方程组，因此用追赶法很容易求解。而对于（）式和（）式导出的系数矩阵不是三对角矩阵，因此求解就不容易了。我们对于显式格式和隐式格式的分析知道，在实际使用上都受到限制，因此构造每层计算量不大的绝对稳定的格式就成为一个具有现实意义并很有兴趣的问题。在一维中，隐式格式是绝对稳定并可用追赶法很容易求解。由此产生了下面将使用的交替方程隐式格式。它具有绝对稳定、容易求解和有相当精度的特点。建立方程组我们在构造微分方程（）的隐式格式和显式格式中，对和做了同样的处理，即同时在第层或第层取值。为了构造一维形式的隐式格式，对二阶导数用在第层上用未知的二阶中心差分来代替，而则用在第层上用已知的二阶中心差分来代替，这样得到的方程组在仅方向是隐式的。比较容易求解，用追赶法就可以了。同理，为对称起见，在下一时间层上重复上述步骤，即又仅在方向是隐式的，对方向是显式的。这样相邻的两个时间层合并起来构成一个差分格式。故用多次追赶法就可以解出了。我们解向后差分格式的方程。令，则（）式变为:（）（）在x方向是隐式时，代入（），变形为（）（）化为矩阵形式：（）（）在y方向是隐式时，代入（），可变形为（）（）化为矩阵形式：（）第5章MATLAB编程追赶法对于前面求出的（）和（）矩阵，我们称之为三对角矩阵。在一些实际问题中，例如解常微分方程的边值问题，接热传导方程一级船体数学放样中建立三次样条函数等，都会要求解系数矩阵为对角占优的三对角方程组，解三对角矩阵我们有一种特殊的方法叫追赶法。例如：求解线性方程组，其中为三对角矩阵。首先对矩阵做如下三角分解：其中为下三角矩阵，为上三角矩阵，且主对角线元素全为1.求解等价于求解两个三角形方程组。先求解的解，再求解的解。追赶法公式实际上就是把高斯消去法用到求解三对角方程组上去的结果。这是由于特别简单，因此使得求解的计算公式非常简单，而且计算量很小。另外追赶法的计算公式中不会出现中间结果数量级的巨大增长和舍入误差的严重积累。最后在MATLAB中编程实现。追赶法程序：见附录。用追赶法求解线性方程组解：在MATLAB命令窗口中输入求解程序：>>A=[2-100;-13-20;0-24-3;00-35]A=2-100-13-200-24-300-35 >>f=[61-21]’;>>x=chase(A,f)输出计算结果为：x=91210.即追赶法求出了方程组的解为：二维热传导方程有限差分解的MATLAB编程程序见附录。第6章数值举例分析算例为了便于观察和表示微分方程的解，我们将最终解用图像表示出来。求满足下列条件的二维热传导方程的解：在MATLAB命令窗口中输入求解程序：a=;u_xy0=inline('0','x','y');u_xyt=inline('x^2*sin(y)-y^2*sin(x)','x','y','t');D=[0,5,0,5];T=1000;Mx=50;My=50;N=50;[u,x,y,t]=sjy(a,D,T,u_xy0,u_xyt,Mx,My,N);mesh(x,y,u)xlabel('x')ylabel('y')zlabel('U')输出计算结果为：rx=+003ry=+003此解即为差分方程的近似解。求解如下二维热传导方程在MATLAB命令窗口中输入求解程序：a=10;u_xy0=inline('0','x','y');u_xyt=inline('x*sin(pi*y)-y*sin(pi*x)','x','y','t');D=[0,2,0,4];T=500;Mx=80;My=80;N=40;[u,x,y,t]=sjy(a,D,T,u_xy0,u_xyt,Mx,My,N);mesh(x,y,u)xlabel('x')ylabel('y')zlabel('t')输出计算结果为：rx=+005ry=+004运行结果分析利用MATLAB的工具箱PDETOOL求解方程结果为：对上述数值解和精确解进行比较，将两解得图画在一个坐标面内，可以更加直观的观察图形的区别，若两图形相差比较大则说明此算法还有不对或不完善的地方；若两图形重合度比较高，则说明此算法和编程过程比较好，能真实的反应偏微分方程的解的情况。如下图：从上图可以看出结果基本一致，重合度比较高，即说明计算差分方程的的结果比较真实，能够很好的反应真实解的情况。若数值解的误差在一定的范数下满足不等式（为初始层的误差），则热传导方程的数值解是稳定的，可从上图看出，两图几乎重合，说明数值解和精确解的误差非常小。第7章结束语本篇文章利用MATLAB数学软件来求解二维热传导方程，通过对区域的剖分，把微分方程进行离散化，在选择隐式差分格式对微分方程进行转化，得到隐式差分方程，再把差分方程转换成矩阵形式求解，由于矩阵形式是严格对角占优的三对角矩阵，通过追赶法轻松的求解。利用MATLAB对追赶法的思想进行编程实现。因此，可以快速的求解类似的二维热传导方程。而且误差非常小，可以控制在有效的范围内。以上算例表明了此方法的可行性和稳定性。相比一维热传导方程有限差分法的MATLAB实现，二维热传导方程的求解更加深入，实际应用性更强。但求解的基本思想不变。致谢本论文是在导师王其林老师和张学富老师的悉心指导下完成的。导师渊博的专业知识，严谨的治学态度，精益求精的工作作风，诲人不倦的高尚师德，严于律己，宽于待人的崇高风范和朴实无华、平易近人的人格魅力对我影响深远。不仅使我在学术上有了更进一步的深入、专业知识更加巩固，最重要的是在我以后的工作和学习中明白了对于待人接物和为人处事的道理，着实受益匪浅。这将会是我人生经历中一笔重要的财富。文章从选题到完成，每一步都有导师的津津指导，倾注了导师的大量心血。在此，谨向导师表示崇高的敬意和由衷的感谢。本文的顺利完成，离不开各位老师、同学、亲人和朋友的关心与帮助。非常感谢邹昌文老师的指导和帮助；感谢负责管理实验室的老师和同学的帮助；感谢室友的帮助和勉励；感谢哥哥苏佳斌的关心和编程方面的指导；感谢爸爸妈妈的多年教育及对我学业的大力支持；感谢好朋友孙斌的帮助和支持。在此，对以上几位以及其他在背后默默支持和帮助我的老师、同学和朋友表示深深地感谢。没有他们的帮助和支持是不可能完成我的学士论文的。愿他们好人一生平安、幸福。重庆交通大学理学院数学与应用数学06级2班苏佳园 2010年6月参考文献[1]陆金甫，关[M]，北京：清华大学出版社，2006。[2]史[J].西安建筑科技大学学报，2009，24(4)：27-29.[3]马明书，马菊意，[J]，河南师范大学学报(数学与信息科学学院)，2008，23(3)：446-452.[4]王正林，何倩，龚[M],北京：电子工业出版社，2009。[5]王飞，[J].新疆师范大学学报（自然科学版），2003，22(4)：21-27.[6][J],河南大学学报（自然科学版）,2004，34(3)：16-18.[7]万正苏，方春华，[J].湖南理工学院学报（自然科学版），2007，20(3)：12-14.[8]曹钢，王桂珍，[J].山东轻工业学院学报，2005，19(4)：76-80.[9]李傅山，[J]，曲阜师范大学数学科学学院，2006，77(2)：42-43.[10]阳[J]，湖南环境生物职业技术学院学报，2009，3(17)：17-19.[11]蔡国梁，李[J]，江苏大学理学院计算科学系学报，2008，29(3)：273-276.[12]张正林，[J]。上海大学理学院数学系学报，2008，22(1)：21-30.[13]刘相国，周宏宇，张海燕，[J],西安理工大学理学院，2007，27(2)：199-204.[14]谷建涛，佟玉霞，付景红.一阶偏微分方程的Mathematica和Matlab解法比较[J],河北理工大学学报，2008，30(2)：40-42.[15]LILi-kangetal.NumericalSolutionofDierentialEquations[M].Shanghai:FudanUniversity[16]SUNHong-lie.Aclassofhighaccuracyexplicitdifferenceschemeforsolvingheat-conductionequationsofmulti-dimension[J].AppliedMathematicsJournalofUniversities,Edition,1999,14(4):427-432.附录追赶法程序functionx=chase(A,f)n=rank(A);b=zeros(n,1);a=zeros(n-1,1);c=zeros(n-1,1);fori=1:n-1b(i,1)=A(i,i);a(i,1)=A(i+1,1);c(i,1)=A(i,i+1);endb(n,1)=A(i,i);fori=2:nb(i,1)=f(i,1)-(a(i-1,1)/b(i-1,1))*c(i-1,1);f(i,1)=f(i,1)-(a(i-1,1)/b(i-1,1))*f(i-1,1);endx(n)=f(n,1)-b(n,1);fori=(n-1):-1:1x(i)=(f(i,1)-c(i,1)*x(i+1))/b(i,1);end求解热传导方程的程序function[u,x,y,t]=sjy(a,D,T,u_xy0,u_xyt,Mx,My,N)%a为方程系数%D为在x轴和y轴上的边界值,D(1)<=x<=D(2),D(3)<=y<=D(4)%T为时间上限%u_xy0为t=0时的初值%u_xyt为边界取值函数%Mx为x轴的等分段数%My为y轴的等分段数%N为时间轴t的等分段数ox=(D(2)-D(1))/Mx;x=D(1)+[0:Mx]*ox;oy=(D(4)-D(3))/My;y=D(3)+[0:My]*oy;ot=T/N;t=[0:N]*ot;%初始化uforj=1:Mx+1fori=1:My+1u(i,j)=u_xy0(x(