本科毕业设计--求解热传导方程的高精度隐式差分格式

三维非稳态导热问题的高效稳定数值解法介绍随着计算机技术的快速发展，数值模拟已经成为研究热传导问题的一种重要方法。

在工程领域，三维非稳态导热问题的数值解法应用广泛。

本文将详细探讨三维非稳态导热问题的高效稳定数值解法。

理论基础1. 热传导方程热传导方程是描述热量在物质中传递的方程，对于三维非稳态导热问题，可以表示为：∂u=α∇2u∂t其中，u是温度场，t是时间，α是热扩散系数，∇2是拉普拉斯算子。

2. 数值方法为了求解三维非稳态导热问题的数值解，有多种数值方法可供选择，例如有限差分法、有限元法和边界元法等。

本文将重点介绍有限差分法。

有限差分法是一种常用的数值方法，用于离散化偏微分方程。

对于三维非稳态导热问题，可以将空间和时间分别离散化，从而得到离散化的方程组。

通过迭代求解，可以得到数值解。

三维非稳态导热问题的数值解法1. 网格划分在使用有限差分法求解三维非稳态导热问题时，需要将计算区域进行网格划分。

通常采用正交网格的方式，将计算区域划分为多个小立方体。

2. 初始条件和边界条件为了求解三维非稳态导热问题，需要指定初始条件和边界条件。

初始条件用于确定初始温度分布，边界条件用于描述物体表面与外界的热交换过程。

3. 显式差分格式显式差分格式是一种简单但稳定性较差的差分格式，用于求解三维非稳态导热问题。

该方法基于前向差分公式，通过迭代计算得到数值解。

4. 隐式差分格式隐式差分格式是一种复杂但稳定性较好的差分格式，用于求解三维非稳态导热问题。

该方法基于后向差分公式，通过迭代求解线性方程组得到数值解。

数值实验为了验证所提出的数值方法的有效性和稳定性，进行了一系列数值实验。

以下是数值实验的步骤：1.确定计算区域和网格划分方式。

2.指定初始条件和边界条件。

3.选择显式差分格式或隐式差分格式进行数值求解。

4.通过迭代计算得到数值解。

5.与解析解进行比较，评估数值方法的精确度和稳定性。

数值实验结果表明，所提出的方法在求解三维非稳态导热问题方面具有较高的精确度和稳定性。

热传导方程的差分格式汇总1.显式差分格式：显式差分格式是最简单的一种方法，通过将导热方程时间和空间上的导数进行近似，引入差分算子，将方程转化为差分格式。

其中最常见的差分格式有：a. 前向差分法（Forward Difference Method）：利用当前节点和其相邻节点的温度值进行计算。

例如，在一维离散情况下，可以使用公式：u(i,j+1)=u(i,j)+α(u(i+1,j)-2u(i,j)+u(i-1,j))b. 后向差分法（Backward Difference Method）：利用当前节点和其相邻节点的温度值进行计算。

例如，在一维离散情况下，可以使用公式：u(i,j+1)=u(i,j)+α(u(i+1,j+1)-2u(i,j+1)+u(i-1,j+1))c. 中心差分法（Central Difference Method）：利用当前和其相邻节点的温度值进行计算。

例如，在一维离散情况下，可以使用公式：u(i,j+1)=u(i,j)+α(u(i+1,j)-2u(i,j)+u(i-1,j))+β(u(i+1,j)-u(i-1,j))其中α和β是时间和空间步长的比例因子。

2.隐式差分格式：显式差分格式具有较大的稳定性限制。

为了克服这个问题，可以使用隐式差分格式，其中使用下一个时间步长的温度值来求解当前时间步长。

常见的隐式差分格式有：a. C-N差分法（Crank-Nicolson Method）：利用前后两个时间步长的温度值进行计算。

例如，在一维离散情况下，可以使用公式：u(i,j+1)=u(i,j)+0.5α(u(i+1,j+1)-2u(i,j+1)+u(i-1,j+1))+0.5α(u(i+1,j)-2u(i,j)+u(i-1,j))b. 力学模拟法（Finite Element Method）：将空间离散化后，通过引入有限元方法，将热传导问题转化为线性方程组，再通过求解线性方程组得到温度分布。

文献综述信息与计算科学热传导方程差分格式的收敛性和稳定性在实际研究物理问题过程中, 往往能给出问题相应的数学表达式, 但是由于实际物理问题的复杂性, 它的解却一般不容易求出. 由此计算物理应运而生, 计算物理是以计算机为工具, 应用数学的方法解决物理问题的一门应用性学科, 是物理、数学和计算机三者结合的交叉性学科. 它产生于二战期间美国对核武器的研究, 伴随着计算机的发展而发展.计算物理的目的不仅仅是计算, 而是要通过计算来解释和发现新的物理规律. 这一点它与传统的实验物理和理论物理并无差别, 所不同的只是使用的工具和方法. 计算物理早已与实验物理和理论物理形成三足鼎立之势, 甚至有人提出它将成为现代物理大厦的“栋梁”.在一个物理问题中一个数值解往往比一个式子更直观, 更有价值. 在实际求解方程时, 除了一些特殊的情况下可以方便地求得其精确解外, 在一般情况下, 当方程或定解条件具有比较复杂的形式, 或求解区域具有比较复杂的形状时, 往往求不到, 或不易求到其精确解. 这就需要我们去寻找方程的近似解, 特别是数值近似解, 简称数值解. 这里主要研究的是热传导方程.有限差分法是微分方程和积分微分方程数值解的方法. 其基本思想是把连续的定解区域用有限个离散点构成的网格来代替, 这些离散点称作网格的节点;把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似;把原方程和定解条件中的微商用差商来近似, 积分用积分和来近似, 于是原微分方程和定解条件就近似地代之以代数方程组, 即有限差分方程组, 解此方程组就可以得到原问题在离散点上的近似解. 然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解.热传导的差分法是求解热传导方程的重要方法之一. 对于差分格式的的求解, 我们首先要关注差分格式的收敛性和稳定性. 对于一个微分方程建立的各种差分格式, 为了有实用意义, 一个基本要求是它们能够任意逼近微分方程, 即相容性要求. 一个差分格式是否有用, 就要看差分方程的精确解能否任意逼近微分方程的解, 即收敛性的概念. 此外, 还有一个重要的概念必须考虑, 即差分格式的稳定性. 因为差分格式的计算过程是逐层推进的, 在计算第n ＋1层的近似值时要用到第n 层的近似值 , 直到与初始值有关. 前面各层若有舍入误差, 必然影响到后面各层的值, 如果误差的影响越来越大, 以致差分格式的精确解的面貌完全被掩盖, 这种格式是不稳定的, 相反如果误差的传播是可以控制的, 就认为格式是稳定的. 只有在这种情形, 差分格式在实际计算中的近似解才可能任意逼近差分方程的精确解. 由Lax 等价定理告诉我们, 对于各适定的线性的初值问题, 对相容性的差分逼近来说, 稳定性则是差分方程的解收敛于微分方程的解的充分必要条件. 收敛是差分方程的本质要求, 稳定是差分方程的基本特性, 对于计算的问题来说, 数值稳定性事差分格式必须要具备的条件, 一个不稳定的差分格式, 即使其他方面有很多的优点, 也是不能用来计算的. 可见由于收敛性和稳定性的重要性, 对于他们的研究是非常具有价值的.热传导方程： 2222222.u u u u a t x y z ⎛⎫∂∂∂∂=++ ⎪∂∂∂∂⎝⎭ 一维热传导方程的初边值问题：22200120(0,0),()(0),(),()(0).t x x l u u a x l t t x u x x l u t u t t ϕμμ===⎧∂∂==<<>⎪∂∂⎪⎪ =<<⎨⎪⎪⎪ = =>⎩用, , 及分别表示初边值问题的解及其偏导数及n j u n j u t ∂⎛⎫ ⎪∂⎝⎭22nj u x ⎛⎫∂ ⎪∂⎝⎭(,)u x t (,)u x t t ∂∂在点之值, 表示求解区域内网格节点. 当初边值问题的解在22(,)u x t x ∂∂(,)j n x t (,)j n x t 区域内部适当光滑时, 对任一区域内部的节点利用泰勒展开公式, 然后化简得(,)j n x t 到显示差分格式：1112200220,()()(1,,1),(),()(0,1,2,).n n nn n j j j j j j n n J U U U U U a t x U j x j J U n t U n t n ϕμμ++-⎧--+-=⎪∆∆⎪⎪=∆=⋅⋅⋅-⎨⎪⎪⎪=∆=∆=⋅⋅⋅⎩这里由于差分方程的解与原初边值问题的解一般是不同的, 故用不同的记号表示.U u 明显的用上式近似热传导方程的初边值问题, 所忽略掉的项, 即截断误差是. 记 2()(())O t O x ∆+∆22()t a x λ∆=∆ 其隐式格式： 111110012(12),()(1,,1),(),()(0,1,2,).n n n n j j j j j n n J U U U U U j x j J U n t U n t n λλλϕμμ+++-+⎧-++-=⎪⎪=∆=⋅⋅⋅-⎨⎪=∆=∆=⋅⋅⋅⎪⎩ 其中. 22()t a x λ∆=∆参考文献[1] 谷超豪, 李大潜, 陈恕行等. 数学物理方程[M ]. 北京: 高等教育出版社, 2002.[2] 刘盾. 实用数学物理方程[M ]. 重庆: 重庆大学出版社, 1996.[3] 张锁春. 抛物型方程定解问题的有限差分数值计算[M ]. 北京: 科学出版社, 2010.[4] (美)哈伯曼. 实用偏微分方程[M ]. 北京: 机械工业出版社, 2007.[5] 陆金甫, 关治. 偏微分方程数值解法[M ]. 北京: 清华大学出版社, 2003.[6] K. W. Morton, D. F. Mayers. 偏微分方程数值解[M ]. 北京: 人民邮电出版社, 2006.[7] 戴嘉尊, 邱建贤. 微分方程数值解法[M ]. 南京: 东南大学出版社, 2002.[8] 徐琛梅. 一类非线性偏微分方程差分格式的稳定性分析[J ]. 江西科学, 2008,27(3) :227~230.[9] 张天德, 张希华, 王玮. 偏微分方程差分格式的构造[J]. 山东工业大学学报, 1997,26(2) :245~246.[10] P. Darania and A. Ebadian. A method for the numerical solution of integrodifferentialequations [J]. Applied Mathematics and Computation , 2007, 188(1): 657~668.[11] Yang Zhang. A finite difference method for fractional partial differential equation [J].Journal of Computational and Applied Mathematics, 2009, 215(2):524~529.。

第1章前言1.1问题背景在史策教授的《一维热传导方程有限差分法的MATLAB实现》和曹刚教授的《一维偏微分方程的基本解》中，对偏微分方程的解得MATLAB实现问题进行过研究，但只停留在一维中，而实际中二维和三维的应用更加广泛。

诸如粒子扩散或神经细胞的动作电位。

也可以作为某些金融现象的模型，诸如布莱克-斯科尔斯模型与Ornstein-uhlenbeck过程。

热方程及其非线性的推广形式也被应用与影响分析。

在科学和技术发展过程中，科学的理论和科学的实验一直是两种重要的科学方法和手段。

虽然这两种科学方法都有十分重要的作用，但是一些研究对象往往由于他们的特性（例如太大或太小，太快或太慢）不能精确的用理论描述或用实验手段来实现。

自从计算机出现和发展以来，模拟那些不容易观察到的现象，得到实际应用所需要的数值结果，解释各种现象的规律和基本性质。

科学计算在各门自然科学和技术科学与工程科学中其越来越大的作用，在很多重要领域中成为不可缺少的重要工具。

而科学与工程计算中最重要的内容就是求解科学研究和工程技术中出现的各种各样的偏微分方程或方程组。

解偏微分方程已经成为科学与工程计算的核心内容，包括一些大型的计算和很多已经成为常规的计算。

为什么它在当代能发挥这样大的作用呢？第一是计算机本身有了很大的发展；第二是数值求解方程的计算法有了很大的发展，这两者对人们计算能力的发展都是十分重要的。

1.2问题现状近三十年来，解偏微分方程的理论和方法有了很大的发展，而且在各个学科技术的领域中应用也愈来愈广泛，在我国，偏微分方程数值解法作为一门课程，不但在计算数学专业，而且也在其他理工科专业的研究生的大学生中开设。

同时，求解热传导方程的数值算法也取得巨大进展，特别是有限差分法方面，此算法的特点是在内边界处设计不同于整体的格式，将全局的隐式计算化为局部的分段隐式计算。

而且精度上更好。

目前，在欧美各国MATLAB的使用十分普及。

在大学的数学、工程和科学系科，MATLAB苏佳园：二维热传导方程有限差分法的MATLAB实现被用作许多课程的辅助教学手段，MATLAB也成为大学生们必不可少的计算工具，甚至是一项必须掌握的基本技能。

解高维热传导方程的一族高精度的显式差分
格式
1 热传导方程及其差分格式
热传导方程是传统数学物理中最基础和最重要的方程之一，它可以描述物体温度随时间、空间变化的过程。

该方程最早出现在18格仑偏微分方程当中。

由于它与现实生活息息相关，自20世纪以来，它发展成为热传导理论的基础，以及热传导问题的基本处理方法和工具。

同时也是热科学及工程中最重要的模拟问题之一。

高维热传导方程有分量形式和平均值形式，它关系到很多跨越学科的问题，是普通微分方程解的典型应用。

但是，通常的数值方法很难满足它的解的准确性要求，尤其是分量形式的高维热传导方程，计算它的精度更为重要。

为了解决高维热传导方程的精度问题，高精度的显式差分格式发展出来，它利用了正交网格，并用空间参数指数外推算法求解热传导方程。

首先，把分量形式简化为差分表达式，格式化为矩阵形式，采用插值方程构成差分法，然后把位置和时间进行外推；最后对比解答解，得出传热率的数值。

该差分格式提供了解高维热传导方程的精准而可靠的工具，可以有效提高高维热传导的研究的质量与速度。

综上所述，高维热传导方程解的准确性极其重要，而高精度的显式差分格式则为此提供了有力的工具，极大地提升了对高维热传导方程的研究的可能性。

热传导方程的一个高精度8点隐式差分格式谢安来;邱淑芳;黄何露【期刊名称】《江西科学》【年(卷),期】2016(000)001【摘要】In this paper,we mainly consider the finite difference method for an initial-boundary prob-lem of standard heat conduction equation inone-dimensional space. An implicit difference scheme of two layers with eight points for the initial-boundary problem was constructed by using Taylor expan-sion and undetermined parameters method. The local truncation error of the proposed difference scheme is Q(τ3 +h6),and the stability condition is 0<r≤0. 126 or 0. 203≤r≤1. 352. Finally,two numerical examples were given to verify the effect of the proposed difference scheme.%主要考虑一维标准的热传导方程初边值问题的有限差分解法。

利用Taylor展开与待定系数的方法，构造出一个2层8点隐式差分格式，并得到该差分格式的局部截断误差为Q(τ3+h6)和稳定条件为0＜r≤0.126或0.203≤r≤1.352。

最后，给出了2个数值算例以验证所得差分格式的计算效果。

【总页数】5页(P10-14)【作者】谢安来;邱淑芳;黄何露【作者单位】东华理工大学理学院，330013，南昌;东华理工大学理学院，330013，南昌;东华理工大学理学院，330013，南昌【正文语种】中文【中图分类】O241.82【相关文献】1.解抛物型方程的一个高精度隐式差分格式 [J], 马明书;王同科2.求解热传导方程的一个高精度格式 [J], 詹涌强3.解热传导方程的一族高精度隐式差分格式 [J], 谭志明; 罗森月4.解高维热传导方程的一个高精度ADI格式 [J], 马明书; 付立志; 朱霖霖; 谷淑敏5.解高维热传导方程的一个新的高精度显格式 [J], 马明书;马小宁因版权原因，仅展示原文概要，查看原文内容请购买。

新疆大学毕业论文(设计)题目:求解热传导方程的高精度隐式差分格式所属院系：数学与系统科学学院专业：信息与计算科学声明本人郑重声明该毕业论文（设计）是本人在开依沙尔老师指导下独立完成的，本人拥有自主知识产权，没有抄袭、剽窃他人成果，由此造成的知识产权纠纷由本人负责。

声明人（签名）：年月日亚库甫江.买买提同学在指导老师的指导下，按照任务书的内容，独立完成了该毕业论文（设计），指导教师已经详细审阅该毕业论文（设计）。

指导教师（签名）：年月日新疆大学毕业论文（设计）任务书班级：信计07-2 姓名：亚库甫江.买买提论文（设计）题目：求解热传导方程的高精度隐式差分格式专题：毕业设计论文（设计）来源：教师自拟要求完成的内容：学习和掌握一维热传导方程已有的各种差分格式的基础上，扩散方程对空间变量应用紧致格式离散，对时间变量应用梯形方法，构造热传导方程的精度为()24τ+数值格式，O h讨论格式的稳定性，最后数值例子来验证。

发题日期：2012 年12月25日完成日期：2012 年5月28 日实习实训单位：数学学院地点：数学学院论文页数：19页；图纸张数：4指导教师：开依沙尔老师教研室主任院长（系主任）摘要本文首先对热传导方程经典差分格式进行复习和讨论，然后热传导方程对空间变量四阶紧致格式进行离散,时间变量保持不变，把一维热传导方程转化为常微分方程组的初值问题, 再利用梯形方法构造热传导方程方程的时间二阶空间四阶精度的一种差分格式，并稳定性进行分析，数值结果与Crank-Nicholson 格式进行比较，数值结果表明, 该方法是有效求解热传导方程的数值计算.关键词: 热传导方程,高精度紧致格式; 梯形方法；两层隐格式; Crank-Nicolson格式ABSTRACTThis paper first study on some classical finite difference for the heat conduction equation, secondely secondely we apply compact finite difference approximation of fourth order for discretizing spatial derivatives but leave the time variable Continuous. This approach results in a system of ODEs, which can then be used trapezodial formula derived fourth order in space and second order in time unconditionally stable implicit scheme .the stability and local truncation error of the obtained method are analysied. Numerical experiments shows that this method Useful, efficient method for solving diffusion equationKeywords: Heat conduction eqution;Higher- oder compact scheme; Trapezodial formula ;Two- level implict scheme; Crank- Nicolson scheme目录引言 (1)预备知识 (2)1.扩散方程的经典差分格式 (3)1.1 显式差分格 (3)1.1.1 显式的截断误差................ . (4)1.1.2 显式差分格式的稳定性 (4)1.2 隐式差分格式 (5)1.2.1 隐式差分格式的截断误差 (5)1.2.2 隐式差分格式的稳定性 (6)1.3 Crank-Nicolson格式 (6)1.3.1 Crank-Nicolson差分格式的截断误差 (7)1.3.2 Crank-Nicolson差分格式的稳定性 (8)2.高精度格式的构造 (9)2.1梯形方法 (9)2.2本文格式的构造 (10)2.3 稳定性分析 (11)3.数值实验 (13)结论 (17)致谢 (18)参考文献 (19)引言热传导方程是一类描述物理量随时间的扩散和衰减规律的抛物型微分方程．自然环境、工程设备及生物机体中的许多物理现象，诸如气体的扩散、液体的渗透、热的传导、以及半导体材料中杂质的扩散等都可用热传导方程来描述．由于物理问题本身的复杂性，其精确解往往不容易求得，因此研究其数值求解方法无疑具有非常重要的理论意义和工程应用价值【1】．求解该问题的数值方法主要有 差分法、有限元法、边界元法等，其中有限差分方法数值求解扩散方程的应用广泛的有效地方法之一。

毕业论文（设计) 题目: 热传导方程初边值问题的差分解法院（系): 数学与计算机科学学院 _专业年级： 2008级数学与应用数学系姓名: XXX __ _学号: 200808101134 __ _指导教师： XXX______ _2012年5月摘要文章目的是为了探讨热传导方程初边值问题的差分解法。

本文包括以下两部分主要内容：第一部分即是对比传统热传导方程初边值问题的变量分离法的差分解法；第二部分即是热传导方程初边值问题差分解法的具体例子。

其中主要涉及到的方法有热传导方程初边值问题的分离变量法和有限差分法.那么先具体介绍有限差分法。

基本思想是把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点;把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似,于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解.然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解. 在采用数值计算方法求解偏微分方程时，若将每一处导数由有限差分近似公式替代，从而把求解偏微分方程的问题转换成求解代数方程的问题，即所谓的有限差分法。

有限差分法求解偏微分方程的步骤如下:1。

区域离散化，即把所给偏微分方程的求解区域细分成由有限个格点组成的网格;2.近似替代,即采用有限差分公式替代每一个格点的导数；3.逼近求解。

换而言之，这一过程可以看作是用一个插值多项式及其微分来代替偏微分方程的解的过程。

对比与分离变量法，有限差分法有着其特性,方便且更精确的特性.经过下面的一番比较,我们有理由相信有限差分法是大有用途的.关键词: 差分格式步长网络节点截断误差AbstractThe article aims to explore the heat conduction equation initial boundary value problem of the finite difference method。

一维热传导方程的差分法一维热传导方程是描述材料内部温度分布随时间变化的数学模型。

它在许多实际工程问题中起着重要的作用，比如热传导、材料加工、建筑设计等。

差分法是一种用于数值求解偏微分方程的常用方法，其原理是将偏微分方程中的导数项用差分近似代替，然后将求解区域划分为离散点，最终得到一个代数方程组。

本文将介绍一维热传导方程的差分法求解过程。

一维热传导方程可以写成如下形式：\[\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}\]\(u(x, t)\)表示材料内部温度分布，\(x\)是空间坐标，\(t\)是时间，\(\alpha\)是热扩散系数。

为了使用差分法求解该方程，我们需要对空间和时间进行离散化。

假设求解区域为\(0 \leq x \leq L\)，时间区间为\(0 \leq t \leq T\)，将空间和时间分别划分成\(N_x\)和\(N_t\)个小区间，步长分别为\(\Delta x = \frac{L}{N_x}\)和\(\Delta t = \frac{T}{N_t}\)。

接下来，我们将使用显式差分格式对一维热传导方程进行离散化。

我们定义离散点\(u_i^n = u(i\Delta x, n\Delta t)\)，用\(u_i^n\)表示时间\(n\)、空间\(i\)处的温度。

那么热传导方程可以用差分格式表示为：\[\frac{u_i^{n+1} - u_i^n}{\Delta t} = \alpha \frac{u_{i+1}^n - 2u_i^n +u_{i-1}^n}{\Delta x^2}\]为了进行数值求解，我们需要给定初始条件和边界条件。

初始条件可以表示为：\[u_i^0 = f(i\Delta x)\]边界条件可以是温度固定或热传导定律，比如：\[u_0^n = g_1(t), u_{N_x}^n = g_2(t)\]或者\[\frac{\partial u}{\partial x}(0, t) = 0, \frac{\partial u}{\partial x}(L, t) = 0\]接下来，我们可以通过迭代计算离散点的温度值来求解一维热传导方程。

热传导方程的差分解法物理学中对热传导现象和扩散现象等物理过程的描述, 通常采用二阶偏微分方程, 统称为热传导方程.9.1. 热传导方程概述一般而言, 在介质内部传导的热量与传热时间、传热截面及温度梯度成正比. 设t 时刻, 点(),,x y z 处的温度为(),,,u x y z t , 则t ∆时间内通过S ∆横截面积传导的热量为(),,,uQ k x y z t t S n∆∆∆∂=-∂其中(),,,0k x y z t >, 是介质的热传导系数. un ∂∂是温度沿S ∆面的法向微商, 即温度梯度的法向分量. 为讨论热传导的规律, 设在介质中任取一小区域V , 其边界面S 为一封闭曲面. 现讨论自1t 至2t 时间内, 小体积V 内热量变化的情况. 首先, 通过包面S 传入V 的热量为()211,,,t t S u Q dt k x y z t ds n ∂=∂⎰⎰⎰ 由矢量积分定理可得()211,,,t t VQ dt k x y z t u dV =∇⋅∇⎡⎤⎣⎦⎰⎰⎰⎰ 其中∇是哈密顿算子.设介质的比热容为c , 密度为ρ, 则V 内温度变化所消耗的热量为212t t V u Q dt c dV tρ∂=∂⎰⎰⎰⎰设体积V 内部热源密度为(),,,F x y z t , 其物理意义是, t 时刻, 点(),,x y z 处, 单位体积热源在单位时间内产生的热量. 所有内部热源产生的热量为()213,,,t t VQ dt F x y z t dV =⎰⎰⎰⎰由能量守恒定律, 即213Q Q Q =+可得()2110t t Vu Q dt c k u F dV t ρ∂⎡⎤=-∇⋅∇-=⎢⎥∂⎣⎦⎰⎰⎰⎰因为体积和时间都是任取的, 所以有 ()u c k u F t ρ∂=∇⋅∇+∂ (9.1) 式(9.1)称为各向同性介质有热源的热传导方程, 也叫做三维非齐次热传导方程. 为简单起见, 设介质是均匀的, 即c 、ρ和k 都是常量. 再设体积V 内无热源, 即(),,,0F x y z t =, 则有u c k u t ρ∂=∆∂ (9.2) 式(9.2)称为各向同性介质无热源的热传导方程, 也叫做三维齐次热做传导方程. 其中∆是拉普拉斯算子. 式(9.2)也可表示为2222222u u u u a t xy z ⎛⎫∂∂∂∂=++ ⎪∂∂∂∂⎝⎭ (9.3)其中2k a c ρ=. 9.2. 一维热传导方程的差分解法各向同性介质中无热源的一维热传导方程为22220,0u u a a t T t x ∂∂=><≤∂∂ (9.4) 其中T 表明时间的有限范围. 要求解方程(9.4), 需要一定的初始条件和边界条件, 统称为定解条件.9.2.1 初值问题()(),0u x x x ϕ=<+∞ (9.5)即初始时刻空间各点的温度颁布函数.9.2.2 初、边值混合问题初始条件为()(),00u x x x l ϕ=≤≤ (9.6)0x =和x l =两端的边界条件有三种情况:第一类边界条件()()()()120,0,u t g t t u l t g t =⎧⎪≥⎨=⎪⎩(9.7) 第二类边界条件()()()()120,0,u t g t xt T u l t g t x∂⎧=⎪⎪∂≤≤⎨∂⎪=⎪∂⎩ (9.8) 其中()1g t 、()2g t 为给定函数.第三类边界条件()()()()()()()()11220,0,0,,u t t u t g t xt T u l t t u l t g t xλλ∂⎧-=⎪⎪∂≤≤⎨∂⎪-=⎪∂⎩ (9.9) 其中1λ、2λ、()1g t 、()2g t 为给定函数, 其中10λ≥, 20λ≥, 且不同时为零.用差分方法求解偏微分方程式(9.4), 首先要建立差分格式. 通常取空间步长和时间步长均为常量. 设空间步长为h , 时间步长为τ, 计算时的步序号空间用i 表示, 时间用k 表示.定义一阶向前商近似为1kk k i i i u u u xh ++∂-=∂一阶向后差商近似为1k k k i i i u u u xh--∂-=∂二阶中心差商作为二阶微商近似为21122,2k k k i i i i k u u u ux h +--+∂=∂ (9.10) 对时间的一阶差分近似为1,k k i i i k u u ut τ+-∂=∂ (9.11) 将(9.10)和式(9.11)代入(9.4), 并令22a h τα=(9.12)即可得一维热传导方程的差分格式为()111121,2,,10,1,,k k k k i i i i u u u u i N k Mααα++-=+-+=-= (9.13)其中,l T N M h τ⎡⎤⎡⎤==⎢⎥⎢⎥⎣⎦⎣⎦, “[]”表示取整.定解条件为()()()()()()001211,2,,11,11,,1i kk Nu i h i N u g k u g k k M ϕττ=-=+=-=-=+差分公式(9.13)为显式格式, 可由初始条件和边界条件逐次计算出任一时刻各点的温度. 习惯上把时刻计算的各点称为一层, 而计算则是一层一层进行的. 计算过程中层间各点的关系如图9.1所示. 从图中可直观地看出, 1k +时刻第i 个点的值是由k 时刻1i -, i 和1i +三点的值算出来.由于初始条件和边界条件的误差及其计算中的舍入误差, 用式(9.13)计算出的值并非该式的精确解k i u . 设计算值与其精确之间的误差为k i ε, 若当k 增加时, k i ε有减小的趋势, 或至少不增加, 则称其差分格式为稳定差分格式. 可以证明, 对于一维热传导方程, 差分格式(9.13)为稳定差分格式的充分条件是2212a h τα=≤(9.14) 差分格式(9.13)计算的具体步骤如下: 1. 给定2,,,,a l h T α2. 计算初始值: ,l T N M h τ⎡⎤⎡⎤==⎢⎥⎢⎥⎣⎦⎣⎦, 计算22h aατ=3. 计算初始值:()()011,2,,1i u i h i N ϕ=-=+ ;计算边界值:()()()()0121,11,,1k k N u g k u g k k M ττ=-=-=+ 4. 用差分格式(9.13)计算1k i u +. 泛定方程2201,0u ux t t x ∂∂=<<<∂∂初始条件()(),04101u x x x x =-≤≤边界条件()()0,001,0u t t u t =⎧⎪≤⎨=⎪⎩程序设计: clear %设置参数 lambda=1; alpha=1/6; L=1; h=0.01; T=0.6;tao=alpha*h^2/lambda; N=fix(L/h); M=fix(T/tao);%设置u 矩阵及x 的值 I=N+1; K=M+1; for i=1:I x(i)=(i-1)*h; endu(I,K)=zeros; %设置初始条件 u(:,1)=4.*x.*(1-x);%设置左端第一类边界条件 u(1,:)=0;%设置右端第一类边界条件 u(I,:)=0; %计算矩阵u for k=1:K for i=2:I-1u(i,k+1)=1/6*u(i+1,k)+2/3*u(i,k)+1/6*u(i-1,k); end end %u ； for k=1:1000:Kplot(x,u(:,k),'-k','LineWidth',2) hold on endx/cmT C Oaxis([0,1,0,1])xlabel ('\fontsize{14}\bfx/cm') ylabel ('\fontsize{14}\bfTC^O') grid on8.6 一维扩散方程的有限差分格式8.6.1 隐式六点差分格式(C —N 格式)以下介绍一维扩散方程或热传导方程的有限差分解法, 考虑一维扩散方程的定解问题()()()()()()()22max 2002111222,,0,0,0t u x t u x t a x l t t t x u x t u x k a u c a u b c x n ua ubc x l n ρ=⎧∂∂=≤≤<<⎪∂∂⎪=⎪⎛⎫⎪=⎨ ⎪∂⎝⎭+==⎪∂⎪⎪∂+==⎪∂⎩ (8.62) 取,x h t ∆∆τ==进行离散化, 如图8.12所示, 结点坐标为()()()()11,2,11,2,i kx i h i N t k k K τ=-=⎧⎪⎨=-=⎪⎩ (8.63) 结点处的函数为(),ki k i u x t u =. 在(),12i k +点, u t∂∂用中心差商,22ux∂∂用(),i k 和(),1i k +两点的中心差商的平均值代替, 则(8.62)中的偏微分方程变为()()()1111111121222k k k k k k k k i i i i i i i i u u u u u u u u hλτ+++++-+-⎡⎤-=-++-+⎣⎦(8.64) 引入212211,1,1a P P P h P Pτ==+=-, 将上式中的含()1k u +项移至等号左边, 将含()ku 项移至等号右边, 式(8.64)变为11111112122k k k k k k i i i i i i u Pu u u Pu u ++++-+--+-=++ (8.65) 上式表明由k 时的值可求得1k +时的u 值, 但要解联立方程组, 所以这种差分格式是隐式的. 整个方程涉及到六个点处的u 值, 所以称为隐式差分格式, 又称为Crank_Nicolson 格式, 简称C_N 格式, 误差为()()22O O h τ+, 是无条件稳定的.8.6.2 边界条件的差分格式由式(8.62)知, 一维扩散方程的边界条件为()()11122200u a u b c x n u a u b c x a n ∂⎧+==⎪⎪∂⎨∂⎪+==⎪∂⎩(8.66)在x 轴上设置两个虚格点0i =和1i N =+(见图8.13). 用中心差商代替.66)中的un∂∂, 则得()()1110212211222N N N b a u u u c hb a u u uc h +-⎧+-=⎪⎪⎨⎪+-=⎪⎩(8.67) 由式(8.67)解出011111222u hc b ha u b u =-+tk +k 图8.12和12222122N N N u hc b ha u b u +-=-+,代入1i =的式(8.65), 有()()11112111111122211111121111121111112121211111222222222242k k k k k k k k k k k k k k k k u Pu hc ha u u u P u hc ha u b u b ub Puha ub ub u b P u hc ha u b u ++++++++-+--+=++-+-++-=++-+整理得到()()1111111212111212k k k kb P ha u b u b P ha u b u hc +++-=-++ (8.68(a)) 同理, 代入i N =的式(8.65), 得到 ()()()()1111111211111222211122221211111222121212212222222222222222k k k k k k N N N N N N k k k k k k k kN N N N N N N N k k k k k k k N N N N N N N u Pu u u P u u hc b ha u b u Pu u hc ha u b u P u u hc ha u b u b Pu b u hc ha u u b P u b ++++-+-++++----++++----+-=++--++-=-+++--++-=-+++21k N u - 整理得()()11212122122222k k k kN N N N b u b P ha u b u b P ha u hc ++---++=+-+ (8.68(b))8.6.3 差分方程组及其求解把式(8.65)和式(8.68(a))和(8.68(b))结合起来, 构成差分方程组, 其形式为AU R = (8.69)其中, ()12,,N U u u u = 是未知量组成的矢量. 系数矩阵A 是三对角的, 而R 是由前一时刻的u 值组成的矢量()12,,N R R R R = . 该方程可利用MA TLAB 求解. 由式(8.65)和式(8.68(a))和(8.68(b))可知()()11211121121212222222k kk k ki i i i k k NN N R b P ha u b u hc R u P u u R b u b P ha u hc-+-⎧=-++⎪=++⎨⎪=+-+⎩ (8.70) 11111112213121121121b P ha b P P A P b b P ha +-⎛⎫ ⎪-- ⎪ ⎪-- ⎪= ⎪ ⎪ ⎪-- ⎪ ⎪-+⎝⎭ (8.70)8.6.4 计算实例研究细杆导热问题. 杆的初始温度是均匀的0u , 保持杆的一端的温度为不变的0u ,至于另一端则有强度为恒定的0q 的热流进入. (解析解见数理方法P214)杆上温度(),u x t 满足下列泛定方程和定解条件(数理方法P214)()20t xx u a u a k c ρ-== (8.71)00x x x lu u u q k ==⎧=⎪⎨=⎪⎩ (8.72) ()000t u u x l ==<< (8.73)边界条件不是齐次的, 首先要处理这个问题. 取一个既满足边界条件(8.72)又满足泛定方程(8.71)的函数(),v x t ,()00,q v x t u x k=+(8.74)计算程序: clear%设置边界条件参数 u0=0; q0k=10; D=1; a1=1.0; b1=0.0; a2=0.0; b2=1.0; c1=u0; c2=q0k;%设置u 矩阵及计算解方程系数 I=101; K=101; h=0.1; tao=0.1; P=tao*D/h^2; P1=1/P+1; P2=1/P-1;for i=1:I x(i)=(i-1)*h; end for k=1:K t=(k-1)*tao; endu(I,K)=zeros; %设置初始条件 u(:,1)=u0;%设置左端第一类边界条件 u(1,:)=u0; %设置系数矩阵A A(I,K)=zeros; A(1,1)=b1*P1+h*a1;x/cmu /u 0A(1,2)=-b1;A(I,K-1)=-b2;A(I,K)=b2*P1+h*a2;for i=2:K-1A(i,i)=2*P1;A(i,i-1)=-1;A(i,i+1)=-1;end%解方程for k=1:K-1R(1,1)=(b1*P2-h*a1)*u(1,k)+b1*u(2,k)+2*h*c1;R(I,1)=b2*u(I-1,k)+(b2*P2-h*a2)*u(I,k)+2*h*c2;for i=2:I-1R(i,1)=u(i-1,k)+2*P2*u(i,k)+u(i+1,k);endc=rank(A)==rank([A R]);u(:,k+1)=A\R;end%作图程序for k=10:20:100plot(x,u(:,k),'-k','LineWidth',2)hold onendaxis([0,10,0,35])xlabel ('\fontsize{14}\bfx/cm')ylabel ('\fontsize{14}\bfu/u_0')grid on%理论结果作图程序clearu0=0;q0k=10;I=101;h=0.1;D=1;for i=1:Ix(i)=(i-1)*h;endl=10;a=sqrt(D);for k=10:20:100t=0.1*k;u=0;for n=1:1000u=u+2*q0k*l/pi^2*(-1).^(n)./(n-1/2)^2*exp(-(n-1/2).^2*pi^2*a^2.*t/l^2).*sin((n-1/2).*pi.*x/l);endU=u+u0+q0k.*x;plot(x,U,':r','LineWidth',2)hold onendaxis([0,10,0,35])grid on例：clear%设置边界条件参数u0=0;q0k=10;D=1;a1=1.0;b1=0.0;a2=0.0;b2=1.0;c1=u0;c2= 0;%设置u矩阵及计算解方程系数I=101;K=101;h=0.1;tao=0.1;P=tao*D/h^2;P1=1/P+1;P2=1/P-1;for i=1:Ix(i)=(i-1)*h;endfor k=1:Kt=(k-1)*tao;endu(I,K)=zeros;%设置初始条件u(:,1)=-q0k.*x;%设置左端第一类边界条件u(1,:)=u0;%设置系数矩阵AA(I,K)=zeros;A(1,1)=b1*P1+h*a1;A(1,2)=-b1;A(I,K-1)=-b2;A(I,K)=b2*P1+h*a2;for i=2:K-1A(i,i)=2*P1;A(i,i-1)=-1;A(i,i+1)=-1;end%解方程for k=1:K-1R(1,1)=(b1*P2-h*a1)*u(1,k)+b1*u(2,k)+2*h*c1;R(I,1)=b2*u(I-1,k)+(b2*P2-h*a2)*u(I,k)+2*h*c2;for i=2:I-1R(i,1)=u(i-1,k)+2*P2*u(i,k)+u(i+1,k);endc=rank(A)==rank([A R]);u(:,k+1)=A\R;end%作图程序for k=1:10:101plot(x,u(:,k),'-k','LineWidth',2)hold onend%axis([0,10,0,35])xlabel ('\fontsize{14}\bfx/cm')ylabel ('\fontsize{14}\bfu/u_0')grid on%理论结果作图程序clearu0=0;q0k=10;I=101;h=0.1;D=1;for i=1:Ix(i)=(i-1)*h;endl=10;a=sqrt(D);for k=1:10:101t=0.1*(k-1);u=0;for n=1:10000u=u+2*q0k*l/pi^2*(-1).^(n)./(n-1/2)^2*exp(-(n-1/2).^2*pi^2*a^2.*t/l^2).*sin((n-1/2).*pi.*x/l);endU=u;plot(x,U,':r','LineWidth',2)hold onend%axis([0,10,0,35])grid onx/cmu /u 0热传导方程的混合问题泛定方程2201,0u u x t t x ∂∂=<<<∂∂初始条件 ()(),04101u x x x x =-≤≤边界条件()()0,001,0u t t u t =⎧⎪≤⎨=⎪⎩ 问题的数值解.clear%设置边界条件参数u0=0;D=1;a1=1.0;b1=0.0;a2=1.0;b2=0.0;c1=u0;c2=u0;%设置u 矩阵及计算解方程系数I=101;K=101;h=0.01;tao=0.01;P=tao*D/h^2;P1=1/P+1;P2=1/P-1;for i=1:Ix(i)=(i-1)*h;endfor k=1:Kt=(k-1)*tao;endu(I,K)=zeros;%设置初始条件u(:,1)=4.*x.*(1-x);%设置左端第一类边界条件u(1,:)=u0;%设置右端第一类边界条件u(101,:)=u0;x/cmu /u 0%设置系数矩阵AA(I,K)=zeros;A(1,1)=b1*P1+h*a1;A(1,2)=-b1;A(I,K-1)=-b2;A(I,K)=b2*P1+h*a2;for i=2:K-1A(i,i)=2*P1;A(i,i-1)=-1;A(i,i+1)=-1;end%解方程for k=1:K-1R(1,1)=(b1*P2-h*a1)*u(1,k)+b1*u(2,k)+2*h*c1;R(I,1)=b2*u(I-1,k)+(b2*P2-h*a2)*u(I,k)+2*h*c2;for i=2:I-1R(i,1)=u(i-1,k)+2*P2*u(i,k)+u(i+1,k);endc=rank(A)==rank([A R]);u(:,k+1)=A\R;end%作图程序for k=1:5:101plot(x,u(:,k),'-k','LineWidth',2)hold onend%axis([0,10,0,35])xlabel ('\fontsize{14}\bfx/cm')ylabel ('\fontsize{14}\bfu/u_0')grid on泛定方程2201,0u u x t t x ∂∂=<<<∂∂初始条件 ()()(),0sin 4101u x x x x =-≤≤边界条件 ()()0,001,0u t t u t =⎧⎪≤⎨=⎪⎩问题的数值解.clear%设置参数lambda=1;alpha=1/6;L=1;h=0.01;T=0.6;tao=alpha*h^2/lambda;N=fix(L/h);M=fix(T/tao);%设置u 矩阵及x 的值I=N+1;K=M+1;for i=1:Ix(i)=(i-1)*h;endu(I,K)=zeros;%设置初始条件u(:,1)=sin(4*pi.*x.*(1-x));%设置左端第一类边界条件u(1,:)=0;%设置右端第一类边界条件u(I,:)=0;%计算矩阵ufor k=1:Kfor i=2:I-1 u(i,k+1)=1/6*u(i+1,k)+2/3*u(i,k)+1/6*u(i-1,k); endendu;for k=1:100:1000plot(x,u(:,k),'-k','LineWidth',2)hold onend%axis([0,1,0,1])xlabel ('\fontsize{14}\bfx/cm') ylabel ('\fontsize{14}\bfTC^O') grid on。

热传导方程的差分格式第2页一维抛物方程的初边值问题分别用向前差分格式、向后差分格式、六点对称格式，求解下列问题：.:u ；：2ua 2，0 ：：： x ：： 1,.:t ；xu(x,0) =sin 二x, 0 ：： x :: 1u(O,t) =u(1,t) =0, t 02在t =0.05,0.1和0.2时刻的数值解，并与解析解u(x,t)=ef t s in (二x)进行比较1差分格式形式2设空间步长h =1/ N ,时间步长• • 0, T =M •，网比r = • / h .(1) 向前差分格式该问题是第二类初边值问题(混合问题)，我们要求出所需次数的偏微商的函数Eu c2uu(x,t)，满足方程—a—^, 0 ：：：x ：：：1,和初始条件u(x,0)= sin x , 0 x ：：： 1抚ex及边值条件u(0,t) =u(1,t) =0, t 0。

已知sin二x在相应区域光滑，并且在x =0,丨与边值相容，使问题有唯一充分光滑的解。

取空间步长h =1/ N，和时间步长-T /M，其中N,M都是正整数。

用两族平行直线x= j X =( j h0Hj 1 ，和tlNt k =X(k=0,1,|||,M) 将矩形域G二｛0 — x — 1,t — 0｝分割成矩形网络，网络格节点为(X j,t k)。

以G h表示网格内点集合，即位于矩形G的网点集合；G h表示闭矩形G的网格集合；丨h=G h-G h是网格界点的集合。

向前差分格式，即k 1 k k k ku, -u, u, 1 -2比■ u, 4- j二a 亠2亠t ( 1)£hk 1 kU j _Uj[ T2k 1 c k 1a U j 1 -2U jh 2k k kU j 1 _2U j U j 和]f j(3)0 U j=(X j ), U 0 = U N = 0.f i = f(x),k kU j j = (X j ), U o = U N = 0其中，j =1,2，…，N _1,k =1,2,…，M 一1.以r =a ./h 2表示网比。