二维三温热传导方程组的分数步隐式差分格式.
6第六讲 典型模型方程-热传导方的差分格式程
ADE Methods
同一时间步，同时左右扫描
6、Hopscotch Method
Comments
2、Richardaon’s Method: CTCS
3、Simple Implicit(Laasonen) Method
算子表示：
其中：
未知，每一个时间步长需要求解三对角方程组
放大因子：
4、Crank-Nicolson Method: famous
1 N 1 1 uN uN j 1 2u j j 1
where
G 1
3-D: ADI Methods
假设
∆t左移，右端合并即为C-N格式
其中 rx x 2 , ry y 2
at
at
3-D: ADI Methods
4、Splitting or Fractional-step Methods
5、ADE Methods
1-D:
或
Simple Explicit Method
差分格式的放大因子与精确解的放大因子比较
精度高
Simple Explicit Method
Simple Explicit Method: Example
Simple Explicit Method: Example
没有相位误差，幅值误差：1.88%
2-D: Crank-Nicolson Scheme
通常采用迭代方 法，需要比三对 角更多的计算资 源
3、2-D: ADI Methods
time
1/ 2 n 1 / 2 1/ 2 uin , uin 1, j , ui , j 1, j
1 n 1 n 1 uin, j 1 , ui , j , ui , j 1
二维热传导方程的可视化计算
二维热传导方程的可视化计算二维热传导方程是描述二维物体热传导过程的数学模型。
在工程领域中，通过求解二维热传导方程，可以预测物体内部的温度分布，进而进行热设计和优化。
热传导是指物体内部由高温区向低温区传递热量的过程。
二维热传导方程是基于热传导定律和能量守恒定律建立的，它可以描述物体内部温度的时空变化。
二维热传导方程的一般形式如下：∂²T/∂x² + ∂²T/∂y² = α ∂T/∂t其中，T是温度，x和y是空间坐标，t是时间，α是热扩散系数。
为了求解二维热传导方程，需要给定边界条件和初始条件。
边界条件是指在物体表面的温度分布情况，而初始条件是指在初始时刻物体内部各点的温度分布。
通常情况下，我们采用数值方法来求解二维热传导方程，其中最常用的方法是有限差分法。
有限差分法将连续的空间和时间离散化，将二维热传导方程转化为一组离散的代数方程。
在计算机中，可以使用计算软件来实现二维热传导方程的可视化计算。
首先，需要将物体的几何形状离散化为一个个小区域，然后对每个小区域进行温度计算。
在计算过程中，可以使用迭代方法来逐步求解离散方程，直到达到收敛条件。
通过迭代计算，可以得到物体在不同时间点的温度分布情况。
在可视化计算中，可以将温度用不同的颜色表示，从而直观地显示物体内部的温度分布。
通过观察温度分布的变化，可以了解物体的热传导特性，并对其进行优化设计。
除了温度分布的可视化，还可以计算物体的热流量、热传导速率等热学参数。
这些参数对于热设计和工程优化非常重要，可以帮助工程师在设计过程中做出准确的决策。
二维热传导方程的可视化计算在工程领域中具有重要的应用价值。
通过求解二维热传导方程，可以预测物体内部的温度分布，为工程设计提供参考依据。
同时，可视化计算也为工程师提供了直观的数据展示方式，帮助他们更好地理解和分析热传导过程。
第九章_热传导方程的差分解法_郑大昉
类似地,其偏微分用差分近似为: 类似地 其偏微分用差分近似为 近似为
∂ui, j,k ui, j,k+1 − ui, j,k = ∂t τ 2 ∂ ui, j,k ui+1, j,k − 2ui, j,k + ui−1, j,k = 2 ∂x h2 ∂2ui, j,k ui, j+1,k − 2ui, j,k + ui, j−1,k = 2 ∂y h2
∂ui,k ∂x ∂ui,k − ∂x + h
(9-18)
二阶中心差商可近似为 二阶中心差商可近似为: 可近似为
∂2ui,k ∂x
即:
2
=
−
(9-19)
ui+1,k − 2ui,k + ui−1,k ∂2u = 2 2 ∂x i,k h
(9-20)
时间的一阶差商近似为 近似为: 另, 对时间的一阶差商近似为
(9-27)
u(x, y,0) = ϕ(x, y)
(9-28)
其边界条件留待后面给出 边界条件留待后面给出. 留待后面给出
差分方法 仍设空间步长 h 仍设空间步长: 空间步长 时间步长: 时间步长 空间为: 网格. 空间为 N× M 网格
τ
则:
Nh = l,
M =s h
t = kτ , k = 0,1 2,... , x = ih, i = 0,1 N ,..., y = jh, j = 0,1,..., M
∆t
∂u ∆Q = −K(x, y, z, t)∆t∆S ∂n
(9-1)
t1
t2 t1
和
Q =∫ 1
∂u dt ∫∫ K(x, y, z, t) dS ∂n (S)
二维三温热传导方程组的二阶分数步差分格式
１ｄ ｉ Ｖ（ ｇ ａ ） ｒｄ 一 ． ７ ． ） （ 一 ，
ｃ 警＝ ｄ（ ｒ ７ ） 丢ｉ 坦ｄ） （ ， 一 ｖ ａ 一 １ Ｔ ，
式 中 、 和 ， 未 知 函 数 ， 别表 示 电子 温 度 、 子 温 度 和 光 子 温度 ； 、 分 别 为 电子 与 离 是 分 离 ％ 子、 电子 与 光 子 的 能 量 交 换 系数 ， Ｐ为 介 质 的 密 度 常 数 ， 、 和 Ｋ 为 扩 散 系 数 ， 为 已 知 函 Ｋ 均 数 。为 简 单 起 见 ， 文 讨 论 扩 散 系 数 与 未 知 函 数 无 关 的情 形 。假 设 问题 是 正 定 的 ， 本 即满 足
到离 散 日 数 的最 优 阶 先验 误 差估 计 及 稳定性 。 范 美蕾 词 中国 法分 类 号
０ 引
言
在惯 性 约束 聚 变 的 二 维 数 值 模 拟 中 ， 常 要 求 解 包 括 电子 热 传 导 ， 子 热 传 导 和 光 子 热 传 常 离 导 方 程 耦 合 的 辐 射 流 体 力 学 方 程 组 Ｉ 其 中 带 有 不 同 能 量 之 间 交换 的 三 温 热 传 导 方 程 组 的 求 １ ］ 解 占 有重 要 的位 置 。 尚武 等 ［ａ 究 了三 温 热 传 导 方程 组 的 能 适 应 各 种 二 维 拉 格 朗 日网 格 的 符 ２￣ －研
二 维 三 温 热 传 导 方 程 组 的
二 阶 分 数 步 差 分 格 式
谢耐森 薛 恒 李 春
（ 岛海 洋 大 学 数 学 系 ，青 岛 ．６ ０ ３ 青 ２６ ０ ）
摘
要
对 二 维三 温热 传 导 方 程 组提 出一 类 分 数步 有 限 差分 格 式 。利 用 变 分 形 式及 能 量 方 法 ． 得 三温 热 传 导方 程 组 ； 数步 有 限 差分 格 式 ｛ 敛性 ； 定性 分 收 稳 ０２ １ ８ ４ ．２ 文 章 编号 ］ ０ —８ ２ ２ ０ ） ３ ０ ９ ６ ０ １１ ６ （０ ２ ０ — ４ ５０
【文献综述】热传导方程差分格式的收敛性和稳定性
文献综述信息与计算科学热传导方程差分格式的收敛性和稳定性在实际研究物理问题过程中, 往往能给出问题相应的数学表达式, 但是由于实际物理问题的复杂性, 它的解却一般不容易求出. 由此计算物理应运而生, 计算物理是以计算机为工具, 应用数学的方法解决物理问题的一门应用性学科, 是物理、数学和计算机三者结合的交叉性学科. 它产生于二战期间美国对核武器的研究, 伴随着计算机的发展而发展.计算物理的目的不仅仅是计算, 而是要通过计算来解释和发现新的物理规律. 这一点它与传统的实验物理和理论物理并无差别, 所不同的只是使用的工具和方法. 计算物理早已与实验物理和理论物理形成三足鼎立之势, 甚至有人提出它将成为现代物理大厦的“栋梁”.在一个物理问题中一个数值解往往比一个式子更直观, 更有价值. 在实际求解方程时, 除了一些特殊的情况下可以方便地求得其精确解外, 在一般情况下, 当方程或定解条件具有比较复杂的形式, 或求解区域具有比较复杂的形状时, 往往求不到, 或不易求到其精确解. 这就需要我们去寻找方程的近似解, 特别是数值近似解, 简称数值解. 这里主要研究的是热传导方程.有限差分法是微分方程和积分微分方程数值解的方法. 其基本思想是把连续的定解区域用有限个离散点构成的网格来代替, 这些离散点称作网格的节点;把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似;把原方程和定解条件中的微商用差商来近似, 积分用积分和来近似, 于是原微分方程和定解条件就近似地代之以代数方程组, 即有限差分方程组, 解此方程组就可以得到原问题在离散点上的近似解. 然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解.热传导的差分法是求解热传导方程的重要方法之一. 对于差分格式的的求解, 我们首先要关注差分格式的收敛性和稳定性. 对于一个微分方程建立的各种差分格式, 为了有实用意义, 一个基本要求是它们能够任意逼近微分方程, 即相容性要求. 一个差分格式是否有用, 就要看差分方程的精确解能否任意逼近微分方程的解, 即收敛性的概念. 此外, 还有一个重要的概念必须考虑, 即差分格式的稳定性. 因为差分格式的计算过程是逐层推进的, 在计算第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.。
热传导方程三层并行差分格式初始条件的计算
１００１－２４６Ｘ （ ２０１１ ）０４－０４８８－０５热传导方程三层并行差分格式初始条件的计算左风丽１崔霞２袁光伟２（１．北京应用物理与计算数学研究所高性能计算中心，北京１０００８８；２．北京应用物理与计算数学研究所计算物理重点实验室，北京１０００８８）摘要：给出二维热传导问题的三层差分格式初始条件的一种显式计算方法，对于由此形成的内边界预估校正三层并行差分算法，证明稳定性和收敛性定理．并行数值试验表明，方法稳定，且与通常采用隐式格式计算初始条件的方法相比，易于程序实现；与已有的扰动算法相比，能大幅度减小误差．三层差分格式；初始条件；稳定显式计算方法 Ｏ２４１．８２；Ｏ２４６Ａ２０１０－０８－１２２０１１－０２－１８国防基础科研项目（ Ｂ１５２０１１００１１）、中国工程物理研究院科学技术基金（２０１０Ａ０２０２０１０）、计算物理实验室基金、国家自然科学基金（６０９７３１５１，６１０３３００９）、８６３课题（２００９ＡＡ０１Ａ１３４，２０１０ＡＡ０１２３０３）及９７３ 课题（２０１１ＣＢ３０９７０２）资助项目左风丽（１９７１ －），女，河南，副研究员，硕士，主要从事大规模科学与工程并行计算研究．第４期４９０＠＠［ １ ］ Ｔｈｏｍａｓ Ｊ Ｗ． Ｎｕｍｅｒｉｃａｌ ｐａｒｔｉａｌ ｄｉｆｆｅｒｅｎｔｉａｌ ｅｑｕａｔｉｏｎｓ ｆｉｎｉｔｅ ｄｉｆｆｅｒｅｎｃｅ ｍｅｔｈｏｄｓ ［ Ｍ］． Ｓｐｒｉｎｇｅｒ－Ｖｅｒｌａｇ， １９９７．＠＠［２］ Ｒｉｃｈｔｍｙｅｒ Ｒ Ｄ，Ｍｏｒｔｏｎ Ｋ Ｗ．初值问题的差分方法［Ｍ］．袁国兴，杜明笙，王汉强，译．广州：中山大学出版社，１９９２．＠＠［３］李德元，陈光南．抛物型方程差分方法引论［Ｍ］．北京：科学出版社，１９９５．＠＠［４］陆金甫，关治．偏微分方程数值解法［Ｍ］．第二版．北京：清华大学出版社，２００４．＠＠［ ５ ］ Ｙｕａｎ Ｇ Ｗ， Ｚｕｏ Ｆ Ｌ． Ｐａｒａｌｌｅｌ ｄｉｆｆｅｒｅｎｃｅ ｓｃｈｅｍｅｓ ｆｏｒ ｈｅａｔ ｃｏｎｄｕｃｔｉｏｎ ｅｑｕａｔｉｏｎ ［ Ｊ ］． Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｊｏｕｒｎａｌ ｏｆ Ｃｏｍｐｕｔｅｒ Ｍａｔｈｅｍａｔｉｃｓ， ２００３， ８０ （ ８ ） ： ９９５ － ９９９．＠＠［ ６ ］ Ｇüｎｔｅｒ Ｓ， Ｌａｃｋｎｅｒ Ｋ． Ａ ｍｉｘｅｄ ｉｍｐｌｉｃｉｔ－ｅｘｐｌｉｃｉｔ ｆｉｎｉｔｅ ｄｉｆｆｅｒｅｎｃｅ ｓｃｈｅｍｅ ｆｏｒ ｈｅａｔ ｔｒａｎｓｐｏｒｔ ｉｎ ｍａｇｎｅｔｉｚｅｄ ｐｌａｓｍａｓ ［ Ｊ ］． Ｊｏｕｒｎａｌ ｏｆ Ｃｏｍｐｕｔａｔｉｏｎａｌ Ｐｈｙｓｉｃｓ， ２００９， ２２８ ： ２８２ － ２９３．＠＠［ ７ ］ Ｙｕａｎ Ｇ Ｗ， Ｈａｎｇ Ｘ Ｄ， Ｓｈｅｎｇ Ｚ Ｑ． Ｐａｒａｌｌｅｌ ｄｉｆｆｅｒｅｎｃｅ ｓｃｈｅｍｅｓ ｗｉｔｈ ｉｎｔｅｒｆａｃｅ ｅｘｔｒａｐｏｌａｔｉｏｎ ｔｅｒｍｓ ｆｏｒ ｑｕａｓｉ－ｌｉｎｅａｒ ｐａｒａｂｏｌｉｃ ｓｙｓｔｅｍｓ ［Ｊ］． Ｓｃｉｅｎｃｅ ｉｎ Ｃｈｉｎａ Ｓｅｒｉｅｓ Ａ： Ｍａｔｈｅｍａｔｉｃｓ， ２００７， ５０（２） ： ２５３ －２７５．＠＠［ ８ ］ Ｓｈｅｎｇ Ｚ Ｑ， Ｙｕａｎ Ｇ Ｗ， Ｈａｎｇ Ｘ Ｄ． Ｕｎｃｏｎｄｉｔｉｏｎａｌ ｓｔａｂｉｌｉｔｙ ｏｆ ｐａｒａｌｌｅｌ ｄｉｆｆｅｒｅｎｃｅ ｓｃｈｅｍｅｓ ｗｉｔｈ ｓｅｃｏｎｄ ｏｒｄｅｒ ａｃｃｕｒａｃｙ ｆｏｒ ｐａｒａｂｏｌｉｃ ｅｑｕａｔｉｏｎ ［ Ｊ ］． Ａｐｐｌｉｅｄ Ｍａｔｈｅｍａｔｉｃｓ ａｎｄ Ｃｏｍｐｕｔａｔｉｏｎ， ２００７， １８４ ： １０１５ － １０３１．＠＠［ ９ ］ Ｙｕａｎ Ｇ Ｗ， Ｓｈｅｎｇ Ｚ Ｑ， Ｈａｎｇ Ｘ Ｄ． Ｔｈｅ ｕｎｃｏｎｄｉｔｉｏｎａｌ ｓｔａｂｉｌｉｔｙ ｏｆ ｐａｒａｌｌｅｌ ｄｉｆｆｅｒｅｎｃｅ ｓｃｈｅｍｅｓ ｗｉｔｈ ｓｅｃｏｎｄ ｏｒｄｅｒ ｃｏｎｖｅｒｇｅｎｃｅ ｆｏｒ ｎｏｎｌｉｎｅａｒ ｐａｒａｂｏｌｉｃ ｓｖｓｔｅｍ ［ Ｊ ］． Ｊｏｕｒｎａｌ ｏｆ Ｐａｒｔｉａｌ Ｄｉｆｆｅｒｅｎｔｉａｌ Ｅｑｕａｔｉｏｎｓ， ２００７， ２０ （ １ ） ： １ － ２０．＠＠［１０］许士良．计算机常用算法［Ｍ］．第二版．北京：清华大学出版社，１９９５．Ｉｎｉｔｉａｌｉｚａｔｉｏｎ Ｍｅｔｈｏｄ ｉｎ Ｔｈｒｅｅ－ｌａｙｅｒ Ｐａｒａｌｌｅｌ Ｄｉｆｆｅｒｅｎｃｅ Ｓｃｈｅｍｅ ｆｏｒ Ｈｅａｔ ＥｑｕａｔｉｏｎＺＵＯ ＦｅｎｇｌｉＣＵＩ Ｘｉａ ＹＵＡＮ Ｇｕａｎｇｗｅｉ。
常系数扩散方程三层十五点差分格式的稳定性
第28卷㊀第5期2023年10月㊀哈尔滨理工大学学报JOURNAL OF HARBIN UNIVERSITY OF SCIENCE AND TECHNOLOGY㊀Vol.28No.5Oct.2023㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀常系数扩散方程三层十五点差分格式的稳定性刘㊀莹,㊀毕㊀卉(哈尔滨理工大学理学院,哈尔滨150080)摘㊀要:基于向后差分格式和Crank-Nicolson 格式对二维扩散方程提出一种三层十五点隐式差分格式㊂采用泰勒展开求出截断误差,证明了该格式的相容性,接着用傅里叶变换和Von Neumann 条件证明了该格式的无条件稳定性㊂由于三层差分格式需要两层启动条件,在数值实验中,利用二维Saulᶄev 差分格式作为三层十五点隐式差分格式的启动格式㊂数值试验表明Saulᶄev 格式与三层十五点差分格式相结合误差小,精度高,并且网比的变化对误差的影响不大㊂关键词:三层十五点差分格式;二维扩散方程;稳定性;误差估计;Saulᶄev 格式DOI :10.15938/j.jhust.2023.05.018中图分类号:O241.3文献标志码:A文章编号:1007-2683(2023)05-0143-07Stability of Three-level Fifteen-point Difference Scheme for Diffusion Equations with Constant CoefficientsLIU Ying,㊀BI Hui(School of Sciences,Harbin University of Science and Technology,Harbin 150080,China)Abstract :Based on backward difference and Crank-Nicolson scheme,a three-level fifteen-point implicit difference scheme for two-dimensional diffusion equation is proposed.The truncation error is obtained by Taylor expansion,and the compatibility of the scheme is proved,then the unconditional stability of the scheme is proved by Fourier transform and Von Neumann condition.Since the three-level difference scheme needs two-level starting conditions,two-dimensional Saulᶄev difference scheme is used as the starting scheme of three-level fifteen-point implicit difference scheme in numerical experiments.Numerical experiments show that the combination of Saulᶄev scheme and three-level fifteen-point difference scheme has small error and high accuracy,and the change of network ratio has littleeffect on the error.Keywords :three-level fifteen-point difference scheme;two-dimensional diffusion equation;stability;error estimation;Saulᶄev scheme㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀收稿日期:2022-05-22基金项目:黑龙江省自然科学基金联合引导项目(LH2020A015).作者简介:刘㊀莹(1996 ),女,硕士研究生.通信作者:毕㊀卉(1982 ),女,博士,教授,博士研究生导师,E-mail:bihui@.0㊀引㊀言扩散方程是一类反映液体渗透㊁气体扩散㊁半导体材料杂质扩散等现象的数学模型,在化学㊁生物㊁物理等领域都有着非常多的应用㊂因此,扩散方程数值解法的研究具有重要科学价值㊂由于有限差分方法(finite difference method,简称FDM)简单灵活且具有很强的通用性,该方法成为一种求解扩散方程数值解的热门方法[1-9],越来越多的数值解法[10-11]也被广泛关注㊂基本的差分格式大多是一层和二层的,由于多层差分格式一般具有较高的精确度,本文旨在研究三层差分格式[12-17]㊂Zhan 等[18]用待定系数法构造了求解一维热传导方程的三层隐式格式,并研究了该方法的截断误差和稳定性条件㊂Amina 等[19]针对扩散方程提出了一种三层九点隐式差分格式,并证明该差分格式与扩散方程相容,二阶收敛且无条件稳定㊂本文基于三层九点差分格式提出一种三层十五点差分格式,并证明该差分格式与二维常系数扩散方程相容且无条件稳定㊂论文第1节给出了二维扩散方程的三层十五点差分格式;第2节计算了格式的截断误差并判断了相容性;第3节证明了差分格式的稳定性;第4节通过数值实验验证了理论结果;最后在第5节给出了结论㊂1㊀三层十五点差分格式二维常系数扩散方程[20]∂u∂t =a ∂2u ∂x 2+∂2u ∂y 2(),x ,y ɪR ,t >0(1)其中a 是一个正常数㊂初始条件为u (x ,y ,0)=g (x ,y ),x ,y ɪR首先,给出一个近似于微分方程(1)的三层十五点差分格式u n +1jk -u n -1jk 2τ-a 3h 2(δ2x u n +1jk +δ2x u n jk +δ2x u n -1jk +δ2y u n +1jk+δ2y u n jk +δ2y u n -1jk )=0(2)其中:τ为时间步长;h 为空间步长;x 和y 的步长相同,都为h ㊂δ2x u jk =u j +1,k -2u jk +u j -1,kδ2y u jk =u j ,k +1-2u jk +u j ,k -1{2㊀截断误差对式(2)中的各项进行泰勒展开,得到u n +1jk -u n -1jk2τ=∂u ∂t []n jk +13!∂3u ∂t 3[]n jk τ2+ (3)δ2xu n +1jk=∂2u ∂x 2[]n +1jk h 2+112∂4u∂x 4[]n +1jk h 4+ =∂2u ∂x 2[]n jk h 2+112∂4u ∂x 4[]njk h 4+∂∂t ∂2u∂x 2[]njk+[h 212∂4u ∂x 4[]njk ]τh 2+12!∂2∂t 2∂2u ∂x 2[]njk[+h 212∂4u∂x 4[]njk]τ2h 2+(4)δ2xu n jk=∂2u ∂x2[]n jkh 2+112∂4u∂x 4[]njkh 4+ (5)δ2xun -1jk=∂2u ∂x 2[]njk h 2+112∂4u ∂x 4[]njkh 4-∂∂t ∂2u ∂x 2[]n jk +[h 212∂4u ∂x 4[]n jk]τh 2+12!∂2∂t 2∂2u ∂x2[]njk+[h 212∂4u∂x 4[]njk]τ2h 2+ (6)δ2yun +1jk =∂2u ∂y 2[]n +1jkh 2+112∂4u ∂y 4[]n +1jkh 4+ =∂2u ∂y2[]njk h 2+112∂4u ∂y 4[]n jk h 4+∂∂t ∂2u∂y 2[]njk+[h 212∂4u ∂y 4[]njk]τh 2+12!∂2∂t 2∂2u ∂x 2[]njk[+h 212∂4u ∂x 4[]njk ]τ2h 2+(7)δ2y u n jk=∂2u ∂y2[]n jkh 2+112∂4u∂y 4[]njk h 4+ (8)δ2yu n -1jk=∂2u ∂y 2[]n jk h 2+112∂4u∂y 4[]njkh 4-∂∂t ∂2u ∂y 2[]n jk+[h 212∂4u∂y 4[]njk]τh 2+12!∂2∂t 2∂2u ∂x2[]njk+[h 212∂4u∂x 4[]njk]τ2h 2+ (9)把式(3)~式(9)代入差分格式(2),有u n +1jk -u n -1jk 2τ-a 3h2[δ2x u n +1jk +δ2x u n jk +δ2x u n -1jk+δ2yun +1jk +δ2yu n jk+δ2yun -1jk ]=∂u ∂t[][njk-a ∂2u ∂x 2[]njk-a ∂2u ∂y2[]njk]+13!∂3u∂t 3[]njk τ2-a 12∂4u ∂x 4[]njkh 2-a 12∂4u ∂y 4[]njkh 2-a 3∂2∂t 2∂2u ∂x 2[][n jk+∂2u ∂y 2[]njk+h 212∂4u ∂x 4[]njk+h 212∂4u ∂y 4[]n jk]τ2+假设方程(1)的解是光滑的,则有T (x ,y ,τ)=13!∂3u∂t 3[]n jk-a 3∂2∂t 2∂2u ∂x 2[][n jk+[∂2u ∂y 2[]n jk +h 212∂4u ∂x 4[]n jk+h 212∂4u ∂y 4[]n jk]]τ2-a 12∂4u ∂x 4[]njk+a 12∂4u ∂y 4[]njk[]h 2+因此差分格式(2)的截断误差具有二阶精度ο(τ2+h 2)㊂由τң0,h ң0得到T (x ,y ,τ)ң0,因此差分格式(2)与微分方程(1)相容㊂定理1 设u (x ,t )是定解问题的解,u n j 是差分格式的解,如果当时间步长τ和空间步长h 都趋于441哈㊀尔㊀滨㊀理㊀工㊀大㊀学㊀学㊀报㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀第28卷㊀零时有e n j =u (x j ,t n )-u nj ң0那么差分格式是收敛的㊂3㊀稳定性把差分格式(2)写成方便计算的形式12(u n +1jk -u n -1jk )=13aλ(δ2x u n +1jk +δ2x u n jk +δ2x u n -1jk +δ2y u n +1jk +δ2y u n jk +δ2y u n -1jk )其中λ=τh 2,整理得到1-23aλδ2x -23aλδ2y ()u n +1jk =1+23aλδ2x +(23aλδ2y )u n -1jk +23aλ(δ2x +δ2y )u n jk (10)为了证明稳定性,给出了等价于式(10)的两层方程:1-23aλ(δ2x +δ2y )()u n +1jk=1+23aλ(δ2x+δ2y)()vn jk +23aλ(δ2x +δ2y )u n jk ㊀㊀㊀㊀㊀㊀v n +1jk=u njk üþýïïïïïï(11)将δ2x u jk ,δ2y u jk (此处省略上标)代入上述两层差分方程(11),得到u n +1jk-23aλ(δ2x u n +1jk +δ2y u n +1jk )=v njk +23aλ(δ2x v n jk +δ2y v n jk )+23aλ(δ2x u n jk +δ2y u njk )㊀㊀㊀㊀㊀㊀v n +1jk =u njk üþýïïïïïï如果取U njk=(u n jk,v n jk)T,其中U n jk是一个两行一列的矩阵,然后把上面的方程写成向量形式,得到D 1U n +1j +1k +D 2U n +1jk +D 1U n +1j -1k +D 1U n +1jk +1+D 1U n +1jk -1=D 3U n j +1k +D 4U n jk +D 3U n j -1k +D 3U n jk +1+D 3U n jk -1(12)其中D 1=-23aλ000éëêêêùûúúú㊀D 2=1+83aλ001éëêêêùûúúúD 3=23aλ23aλ00éëêêêùûúúú㊀D 4=-83aλ1-83aλ10éëêêêùûúúú如果U n jk =V n e(i (fjh +gkh )),其中V n 与U n 形式相同,都是两行一列矩阵,那么由式(12)有D 1Vn +1eih (f (j +1)+gk )+D 2Vn +1eih (fj +gk )+D 1V n +1e ih (f (j -1)+gk )+D 1V n +1e ih (fj +(k +1)g )+D 1V n +1e ih (fj +(k -1)g )=D 3V n e ih (f (j +1)+gk )+D 4V n e ih (fj +gk )+D 3V n e ih (f (j -1)+gk )+D 3V n e ih (fj +(k +1)g )+D 3V n e ih (fj +(k -1)g )(13)为了简化矩阵,令γf =e ihf +e -ihf ,γg =e ihg +e -ihg ,ω=cos(hf )+cos(hg )将式(13)消去公因式得到1-φ001éëêêùûúúV n +1=φ1+φ10éëêêùûúúV n 其中φ=23aλγf +23aλγg -83aλ㊂由于e ihf =cos(hf )+i sin(hf )e -ihf =cos(hf )-i sin(hf ){,于是有-43aλω+1+83aλ001éëêêêùûúúúV n +1=43aλω-83aλ43aλω+1-83aλ1éëêêêùûúúúV n 令α=43aλω-83aλ=43aλ(cos(hf )+cos(hg ))-83aλ显然αɤ0,得到1-α001éëêêùûúúV n +1=α1+α10éëêêùûúúV n 于是得到增长因子G (τ,k )=1-α001éëêêùûúú-1α1+α10éëêêùûúú=α1-α1+α1-α10éëêêêùûúúúG (τ,k )的特征值函数可以被写成μ2-α1-αμ-1+α1-α=0(14)为了得到格式的稳定性,需要如下引理㊂引理1[20]㊀实系数二次方程aμ2-bμ-c =0的根按模小于等于1的充要条件是:|b |ɤ1-c ɤ2㊂已知αɤ0,根据式(14),有b =α1-α,c =1+α1-α541第5期刘㊀莹等:常系数扩散方程三层十五点差分格式的稳定性并且|b |ɤ1-c =1-1+α1-α=-2α1-α且1-c =-2α1-αɤ2㊂由引理1,可以得到|μi |ɤ1(i =1,2),所以|G |ɤ1㊂这样就满足了Von Neumann 条件,因此差分格式(2)无条件稳定㊂定理2㊀(Von Neumann 条件)差分格式稳定的必要条件是当τɤτ0,nτɤT ,对所有k ɪR 有|λj (G (τ,k ))|ɤ1+Mτ,j =1,2, ,p ,其中λj (G (τ,k ))表示G (τ,k )的特征值,M 为常数㊂4㊀数值算例已知扩散方程∂u∂t=4-2(u xx +u yy )(x ,y )ɪG =(0,1)ˑ(0,1),t >0u (0,y ,t )=u (1,y ,t )=0,0ɤy ɤ1,t >0u (x ,0,t )=u (x ,1,t )=0,0ɤx ɤ1,t >0u (x ,y ,0)=sinπx sinπyüþýïïï(15)通过变量分离可以得到方程的解析解u =sinπx sinπy exp -π28t()(x ,y )ɪG =(0,1)ˑ(0,1),t >0㊂离散方程(15)的定义域:令x j =jh ,y k =kh (j ,k =0,1, ,J ),t n =nτ(n =1,2, ),其中τ为时间步长,网格比λ=τh 2,重新排列(2),得到u n +1jk-23aλu n +1j +1k +43aλu n +1jk -23aλu n +1j -1k -23aλu n +1jk +1+43aλu n +1jk -23aλu n +1jk -1=u n -1jk +23aλu n -1j +1k -43aλu n -1jk +23aλu n -1j -1k +23aλu n -1jk +1-43aλu n -1jk +23aλu n -1jk -1+23aλu n j +1k -43aλu n jk +23aλu n j -1k +23aλu n jk +1-43aλu n jk +23aλu n jk -1令j =1:J -1,k =1:J -1U n =[u n 11,u n 21, ,u n (J -1)1,u n 12,u n22, ,u n (J -1)2, ,u n 1(J -1),u n 2(J -1), ,u n (J -1)(J -1)]T 取J -1阶方阵A ii ,B ii ,C ii ,F 1,F 2,F 3如下:A ii =1+83aλ-23aλ-23aλ1+83aλ-23aλ⋱⋱⋱-23aλ1+83aλ-23aλ-23aλ1+83aλæèçççççççççççöø÷÷÷÷÷÷÷÷÷÷÷B ii =-83aλ23aλ23aλ-83aλ23aλ⋱⋱⋱23aλ-83aλ23aλ23aλ-83aλæèçççççççççççöø÷÷÷÷÷÷÷÷÷÷÷C ii =1-83aλ23aλ23aλ1-83aλ23aλ⋱⋱⋱23aλ1-83aλ23aλ23aλ1-83aλæèçççççççççççöø÷÷÷÷÷÷÷÷÷÷÷F 1=diag -23aλ, ,-23aλ{}F 2=diag 23aλ, ,23aλ{}F 3=diag 23aλ, ,23aλ{}记A =A 11F 1F 1A 22F 1⋱⋱⋱F 1A J -2,J -2F 1F 1A J -1,J -1æèçççççççöø÷÷÷÷÷÷÷B =B 11F 2F 2B 22F 2⋱⋱⋱F 2B J -2,J -2F 2F 2B J -1,J -1æèçççççççöø÷÷÷÷÷÷÷641哈㊀尔㊀滨㊀理㊀工㊀大㊀学㊀学㊀报㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀第28卷㊀C =C 11F 3F 3C 22F 3⋱⋱⋱F 3C J -2,J -2F 3F 3C J -1,J -1æèçççççççöø÷÷÷÷÷÷÷显然A ,B ,C 均为(J -1)2阶方阵㊂用A ,B ,C 作为系数矩阵,于是有AU n +1+M =BU n +CU n -1+Q其中M =m 1m 2︙m J -2m J -1æèçççççççöø÷÷÷÷÷÷÷㊀㊀Q =q 1q 2︙q J -2q J -1æèçççççççöø÷÷÷÷÷÷÷m 1=-23aλu n +10,1-23aλu n +11,0-23aλu n +12,0-23aλu n +13,0︙-23aλu n +1J -2,0-23aλu n +1J ,1-23aλu n +1J -1,0æèççççççççççççççöø÷÷÷÷÷÷÷÷÷÷÷÷÷÷q 1=23aλu n 01+u n 10()+23aλu n -101+u n -110()23aλu n 20+23aλu n -12023aλu n 30+23aλu n -130︙23aλu n J -2,0+23aλu n -1J -2,023aλu n J ,1+u n J -1,0()+23aλu n -1J ,1+u n -1J -1,0()æèççççççççççççççöø÷÷÷÷÷÷÷÷÷÷÷÷÷÷m 2=-23aλu n +10,200︙0-23aλu n +1J ,2æèççççççççççöø÷÷÷÷÷÷÷÷÷÷㊀q 2=23aλu n 0,2+23aλu n -10,200︙023aλu n J ,2+23aλu n -1J ,2æèççççççççççöø÷÷÷÷÷÷÷÷÷÷m 3=-23aλu n +10,300︙0-23aλu n +1J ,3æèççççççççççöø÷÷÷÷÷÷÷÷÷÷㊀q 3=23aλu n 0,3+23aλu n -10,300︙023aλu n J ,3+23aλu n -1J ,3æèççççççççççöø÷÷÷÷÷÷÷÷÷÷m 4=-23aλu n +10,400︙0-23aλu n +1J ,4æèççççççççççöø÷÷÷÷÷÷÷÷÷÷㊀q 4=23aλu n 0,4+23aλu n -10,400︙023aλu n J ,4+23aλu n -1J ,4æèççççççççççöø÷÷÷÷÷÷÷÷÷÷一直到m J -2和q J -2㊂m J -1=-23aλu n +10,J -1+u n +11,J ()-23aλu n +12,J-23aλu n +13,J ︙-23aλu n +1J -2,J -23aλu n +1J ,J -1+u n +1J -1,J()æèçççççççççççççöø÷÷÷÷÷÷÷÷÷÷÷÷÷q J -1=23aλ(u n 0,J -1+u n 1,J )+23aλ(u n -10,J -1+u n -11,J )23aλu n 2,J +23aλu n -12,J ︙23aλu n J -2,J +23aλu n -1J -2,J 23aλ(u n J ,J -1+u n J -1,J )+23aλ(u n -1J ,J -1+u n -1J -1,J )æèçççççççççççöø÷÷÷÷÷÷÷÷÷÷÷显然M =0,Q =0㊂这样有U n +1=A -1BU n +A -1CU n -1(16)对于三层格式,U 0是已知的初始条件,U 1未知,这里选择二维Saulᶄev 差分格式[21]计算U 1:u ɵn +1jk=aλ1+(θ1+θ2)aλθ1u n +1j +1,k +[θ2u n +1j ,k +1+(1-θ1)u n j +1,k +(1-θ2)u n j ,k +1+u nj -1,k +u n j ,k -1-(4-θ1-θ2-1aλ)u n j ,k](17)741第5期刘㊀莹等:常系数扩散方程三层十五点差分格式的稳定性uɵn +1jk=aλ1+(θ1+θ2)aλθ1u n +1j -1,k +[θ2u n +1j ,k -1+(1-θ1)u n j -1,k +(1-θ2)u n j ,k -1+u nj +1,k +u nj ,k +1-(4-θ1-θ2-1aλ)u nj ,k](18)可以把式(17)㊁(18)结合来提高精确度㊂例如,u ɵn +1jk 和u ɵn +1jk 同时满足式(15)给定的初边值条件,则第一层可以由它们的算数平均数U 1=12(u ɵ1+u ɵ1)给出㊂Saulᶄev 格式的截断误差为Oτh+τ+h 2()㊂这里取θ1=θ2=12来计算U 1的值,容易证明矩阵A可逆,于是回到式(16)可得到n =2,3, ,N 各时间层的数值解㊂接下来给出数值解在不同时刻的图像㊂图1和图2分别给出了当λ=1时在T =1以及T =5时刻对应的数值解㊂表1给出了在不同网比下不同时刻的数值解与真解㊂图1㊀λ=1,T =1时数值解Fig.1㊀Numerical solution when λ=1,T =1图2㊀λ=1,T =5时数值解Fig.2㊀Numerical solution when λ=1,T =5表1㊀三层十五点差分格式下T =1和T =5时不同网格比下的数值解和误差Tab.1㊀Numerical solution anderror under different gridratios when T =1and T =5in three-levelfifteen-point difference scheme T =1T =5λ数值解误差λ数值解误差0.20.0126030.0000150.29.1046ˑ10-55.12ˑ10-70.60.0126030.0000150.69.1046ˑ10-55.12ˑ10-710.0126030.00001519.1047ˑ10-55.13ˑ10-730.0126030.00001539.1053ˑ10-55.19ˑ10-750.0126030.00001559.1064ˑ10-55.3ˑ10-75㊀结㊀论本文讨论了常系数扩散方程三层十五点差分格式的稳定性和误差估计问题㊂然后用Saulᶄev 格式求出第一层,再结合三层十五点差分格式求出数值解㊂数值结果表明,Saulᶄev 格式与三层十五点差分格式相结合误差小,精度高㊂λ的变化对误差的影响不大㊂算法的全部处理表明本文所讨论的三层十五点差分方法是无条件稳定㊁可行的㊂参考文献:[1]㊀WANG Tingchun,GUO Boling.Analysis of Some FiniteDifference Schemes for Two-Dimensional Ginzburg-Lan-dau Equation [J].Numer Methods Partial Differential E-quations,2011,27(5):1340.[2]㊀孙志忠.偏微分方程数值解法[M].北京:科学出版社,2005.[3]㊀武莉莉.不可压磁流体力学方程组的高精度紧致有限差分方法[D].宁夏:宁夏大学,2021.[4]㊀余德浩,汤华中.微分方程数值解法[M].北京:科学出版社,2003.[5]㊀张鲁明,常谦顺.复Schrödinger 场和实Klein-Gordon场相互作用下一类方程组守恒差分格式的收敛性和稳定性[J].高等学校计算数学学报,2000(4):362.ZHANG Luming,CHANG Qianshun.Convergence and Stability of a Conservative finite Difference Scheme for aClass of Equation system in Interaction of Complex Schrödinger Field and Real Klein-Gordon Field [J ].Journal of Computational Mathematics,2000(4):362.[6]㊀李小纲.流体力学中双曲守恒律方程的高精度差分方841哈㊀尔㊀滨㊀理㊀工㊀大㊀学㊀学㊀报㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀第28卷㊀法研究[D].西安:西安理工大学,2020.[7]㊀ZHANG Luming.Convergence of a Conservative Differ-ence Schemes for a Class of Klein-Gordon-Schrödinger E-quations in one Space Dimension[J].Applied Mathe-matics and Computation,2000,163(1):343. [8]㊀杨彩杰,孙同军.抛物型最优控制问题的Crank-Nicol-son差分方法[J].山东大学学报:理学版,2020,55(6):115.YANG Caijie,SUN Tongjun.Crank-Nicolson FiniteDifference Method for Parabolic Optimal Control Problem[J].Journal of Shandong University:Science Edition,2020,55(6):115.[9]㊀王廷春,张鲁明,陈芳启,等.求解Klein-Gordon-Schrödinger方程组的一个新型守恒差分算法的收敛性分析[J].高等应用数学学报,2008(1):41.WANG Tingchun,ZHANG Luming,CHEN Fangqi,etal.Convergence Analysis of a New Conservative Differ-ence Algorithm for Solving Klein-Gordon-Schrödinge E-quations[J].Applied Mathematics A Journal of Chi-nese,2008(1):41.[10]毕卉,陈莎莎.四阶线性方程局部间断Galerkin方法的误差估计[J].哈尔滨理工大学学报,2021,26(4):159.BI Hui,CHEN Shasha.Error Estimates for Local Discon-tinuous Galerkin Methods for Linear Fourth-order Equa-tions[J].Journal of Harbin University of Science andTechnology,2021,26(4):159.[11]毕卉,孟雄,孙阳.求解双曲守恒律方程的高阶TVD格式[J].哈尔滨理工大学学报,2010,15(3):54.BI Hui,MENG Xiong,SUN Yang.A Higher Order TVDScheme for Hyperbolic Conservation Laws[J].Journal ofHarbin University of Science and Technology,2010,15(3):54.[12]苏保金,姜子文.二维拟线性粘性波动方程的三层紧致差分格式[J].山东师范大学学报:自然科学版,2019,34(2):171.SU Baojin,JIANG Ziwen.A Three Level Compact Differ-ence Scheme For Solving A Two-Dumensional QuasilinearViscous Wave Equation[J].Journal of Shandong NormalUniversity:Natural Science Edition,2019,34(2):171.[13]李佳佳,张虹,王希,等.Rosenau-KdV-RLW方程的三层线性化差分格式[J].四川大学学报:自然科学版,2018,55(6):1137.LI Jiajia,ZHANG Hong,WANG Xi,et al.A Three-lev-el Linearized Difference Scheme for Rosenau-KdV-RLWEquation[J].Journal of Sichuan University:Natural Sci-ence Edition,2018,55(6):1137.[14]常红,丁丹平.Camassa-Holm方程的一种三层守恒有限差分格式[J].陕西科技大学学报,2017,35(3):186.CHANG Hong,DING Danping.A Three-level Conserva-tive Finite Difference Scheme for Camassa-Holm Equation[J].Journal of Shaanxi University of Science and Tech-nology,2017,35(3):186.[15]赵红伟,胡兵,郑茂波.General Improved KdV方程的三层加权平均线性差分格式[J].四川大学学报:自然科学版,2017,54(1):12.ZHAO Hongwei,HU Bing,ZHENG Maobo.Three-levelAverage Linear Difference Scheme for the General Im-proved KdV Equation[J].Journal of Sichuan University:Natural Science Edition,2017,54(1):12. [16]谢建强.一维粘性波动方程的三层紧致差分格式[J].南昌航空大学学报:自然科学版,2016,30(2):50.XIE Jianqiang.A Three level Compact Difference Schemefor Solving A One-dimensional Viscous Wave Equation[J].Journal of Nanchang Hangkong University:NaturalScience Edition,2016,30(2):50.[17]杜瑜,徐友才,胡兵.Rosenau-Burgers方程的三层差分格式(英文)[J].四川大学学报:自然科学版,2010,47(1):1.DU Yu,XU Youcai,HU Bing.The Three Level FiniteDifference Scheme for Rosenau-Burgers Equation[J].Journal of Sichuan University:Natural Science Edition,2010,47(1):1.[18]詹涌强,凌婷.求解一维热传导方程的一族三层隐格式[J].西南师范大学学报:自然科学版,2020,45(11):1.ZHAN Yongqiang,LING Ting.A Class of Three-level Im-plicit Difference Scheme for Solving One-dimensionalHeat Conduction Parabolic Equations[J].Journal ofSouthwest Normal University:Natural Science Edition,2020,45(11):1.[19]阿米娜㊃沙比尔,杨庆之.常系数扩散方程的三层九点差分格式的稳定性(英文)[J].高等学校计算数学学报,2020,42(2):148.AMINA Shabel,YANG Qingzhi.Stability of Three-levelNine-point Difference Scheme for Constant CoefficientDiffusion Equations[J].Numerical Mathematics A Journalof Chinese Universities,2020,42(2):148. [20]李荣华.偏微分方程数值解法[M].北京:高等教育出版社,2010.[21]孙洁.关于向前-向后热方程的数值方法[D].杭州:浙江大学,2008.(编辑:温泽宇)941第5期刘㊀莹等:常系数扩散方程三层十五点差分格式的稳定性。
二维热传导方程的差分格式
二维热传导方程的差分格式
二维热传导方程是描述热量在二维平面内传导的方程。
在数值计算中，我们通常采用差分法来求解二维热传导方程。
差分法的基本思想是将求解区域划分为若干个小区域，然后在每个小区域内近似求解热传导方程。
具体来说，我们可以采用有限差分法来离散化热传导方程，将偏导数转化为有限差分近似值，然后再用迭代方法求解离散化后的方程组。
常用的差分格式有显式差分格式、隐式差分格式、Crank-Nicolson差分格式等。
其中，显式差分格式计算简单，但是需要满足一定的时间步长条件才能保证稳定性；隐式差分格式较为稳定，但是计算量较大；Crank-Nicolson差分格式结合了显式和隐式的优点，是一种较为稳定且计算量较小的差分格式。
在实际应用中，我们可以根据具体问题的性质选择合适的差分格式来求解二维热传
导方程。
- 1 -。
【毕业设计(论文)】二维热传导方程有限差分法的MATLAB实现
第1章前言1.1问题背景在史策教授的《一维热传导方程有限差分法的MATLAB实现》和曹刚教授的《一维偏微分方程的基本解》中，对偏微分方程的解得MATLAB实现问题进行过研究，但只停留在一维中，而实际中二维和三维的应用更加广泛。
诸如粒子扩散或神经细胞的动作电位。
也可以作为某些金融现象的模型，诸如布莱克-斯科尔斯模型与Ornstein-uhlenbeck过程。
热方程及其非线性的推广形式也被应用与影响分析。
在科学和技术发展过程中，科学的理论和科学的实验一直是两种重要的科学方法和手段。
虽然这两种科学方法都有十分重要的作用，但是一些研究对象往往由于他们的特性（例如太大或太小，太快或太慢）不能精确的用理论描述或用实验手段来实现。
自从计算机出现和发展以来，模拟那些不容易观察到的现象，得到实际应用所需要的数值结果，解释各种现象的规律和基本性质。
科学计算在各门自然科学和技术科学与工程科学中其越来越大的作用，在很多重要领域中成为不可缺少的重要工具。
而科学与工程计算中最重要的内容就是求解科学研究和工程技术中出现的各种各样的偏微分方程或方程组。
解偏微分方程已经成为科学与工程计算的核心内容，包括一些大型的计算和很多已经成为常规的计算。
为什么它在当代能发挥这样大的作用呢？第一是计算机本身有了很大的发展；第二是数值求解方程的计算法有了很大的发展，这两者对人们计算能力的发展都是十分重要的。
1.2问题现状近三十年来，解偏微分方程的理论和方法有了很大的发展，而且在各个学科技术的领域中应用也愈来愈广泛，在我国，偏微分方程数值解法作为一门课程，不但在计算数学专业，而且也在其他理工科专业的研究生的大学生中开设。
同时，求解热传导方程的数值算法也取得巨大进展，特别是有限差分法方面，此算法的特点是在内边界处设计不同于整体的格式，将全局的隐式计算化为局部的分段隐式计算。
而且精度上更好。
目前，在欧美各国MATLAB的使用十分普及。
在大学的数学、工程和科学系科，MATLAB苏佳园：二维热传导方程有限差分法的MATLAB实现被用作许多课程的辅助教学手段，MATLAB也成为大学生们必不可少的计算工具，甚至是一项必须掌握的基本技能。
热传导方程有限差分法的MATLAB实现
△t
n
nn
关于
t
的二阶中心差商[10]：
坠2u 坠x2
≈
uj+1
-2uj +uj-1 (△x)2
，对方
程进行离散。 离散后的方程为：
n n-1
n
nn
uj -uj △t
=a2
uj+1
-2uj +uj-1 (△x)2
。
令
：r=
a2·△t (△x)2
，即
n
n
n
n-1
（1＋2r)uj -r·uj+1 -r·uj-1 =uj 。 可化为矩阵形式：
摘 要：对于有界热传导齐次方程的混合问题，用分离变量法求解往往很复杂。 为了更好地
理解热传导方程的解，使用 MATLAB 软件将方程的解用图像表示出来。 通过区域转换的思想，
利用 MATLAB 编程实现一定区域内热传导方程的有限差分方法，数值表明了方法的可行性和
稳定性。
关键词：热传导方程；有限差分；MATLAB
方法， 把控制方程中的导数用网格节点上的函数值
的差商代替进行离散，从而 建立以网格节点上的值
为未知数的代数方程组。
1 求解热传导方程的基本思想
基本思想是把连续的定解区域用有限个离散点
构成的网格来代替， 这些离散点称作网格的节点；
把连续定解区域上的连续变量的函数用在网格上定
义的离散变量函数来近似； 把原方程和定解条件中
x0(ii+1)=ii*ox; end u=sin(pi*x0/l); % t=0 时 u(x,t)的值 r=a^2*ot/(ox)^2; for ii=1:n
%数据的输入 B=zeros(M-1,1);%存放系数矩阵主对角线元素 A=zeros (M-2,1);%存放系数矩阵主对角线元素下 方次对角线的元素 C=zeros (M-2,1);%存放系数矩阵主对角线元素上 方次对角线的元素 S=zeros(M-1,1);%存放右端的常数项 for ii=1:M-2
求解热传导方程的高精度隐式差分格式毕业论文
新疆大学毕业论文(设计)题目:求解热传导方程的高精度隐式差分格式所属院系：数学与系统科学学院专业：信息与计算科学声明本人郑重声明该毕业论文（设计）是本人在开依沙尔老师指导下独立完成的，本人拥有自主知识产权，没有抄袭、剽窃他人成果，由此造成的知识产权纠纷由本人负责。
声明人（签名）：年月日亚库甫江.买买提同学在指导老师的指导下，按照任务书的内容，独立完成了该毕业论文（设计），指导教师已经详细审阅该毕业论文（设计）。
指导教师（签名）：年月日新疆大学毕业论文（设计）任务书班级：信计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】．求解该问题的数值方法主要有 差分法、有限元法、边界元法等，其中有限差分方法数值求解扩散方程的应用广泛的有效地方法之一。
2020年智慧树知道网课《传热学(华东交通大学)》课后章节测试满分答案
绪论单元测试1【判断题】(10分)导热，对流，辐射换热是热量传递的三种基本方式。
A.错B.对2【判断题】(10分)传热系数和导热系数单位不同。
A.对B.错3【多选题】(10分)下列哪几种传热过程不需要有物体的宏观运动？A.复合传热B.导热C.辐射D.对流4【多选题】(10分)热量传递的三种基本方式为（）。
A.热传导B.热辐射C.热对流D.传热5【单选题】(10分)太阳与地球间的热量传递属于下述哪种传热方式?A.热对流B.导热C.其他几种都不是D.热辐射6【单选题】(10分)温度对辐射换热的影响（）对对流换热的影响。
A.可能大于、小于B.小于C.大于D.等于7【单选题】(10分)物体不论()高低，都在相互辐射能量，只是辐射能量的大小不同。
A.导热B.温度C.热传导D.放热8【单选题】(10分)工程中常遇到热量从固体壁面一侧的高温流体，通过固体壁传递给另一侧低温流体的过程，称为()。
A.传热过程B.热对流C.热传导D.热辐射9【判断题】(10分)热辐射和流体对流及导热一样，需有温差才能发射辐射能。
A.对B.错10【单选题】(10分)传热学就是研究()引起的热量传递规律的学科。
A.焓差B.浓度差C.熵差D.温差第一章测试1【单选题】(10分)导热问题的第一类边界条件是已知（）。
A.温差B.热流密度C.壁温D.对流换热量2【单选题】(10分)下面材料中哪种材料的导热系数最小（）。
A.瓷砖B.铁C.铜D.硅藻土砖3【判断题】(10分)温度梯度表示温度场内的某一地点等温面法线方向的温度变化率。
A.对B.错4【单选题】(10分)表征材料导热能力的物理量是（）。
A.导温系数B.吸热系数C.导热系数D.传热系数5【单选题】(10分)按照导热机理，水的气、液、固三种状态中（）状态下的导热系数最小。
A.气态B.固态C.无法确定D.液态6【单选题】(10分)气体的导热系数随温度的升高而（）。
A.增加B.减小C.无法确定D.不变7【单选题】(10分)一般而言，金属比非金属（介电体）的导热系数值是（)。
热传导方程的差分解法
热传导方程的差分解法物理学中对热传导现象和扩散现象等物理过程的描述, 通常采用二阶偏微分方程, 统称为热传导方程.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。
数值计算方法解决二维热传导方程问题研究
数值计算方法解决二维热传导方程问题研究概述：热传导方程是描述物体中温度分布随时间演化的常见方程之一。
解决热传导方程的问题在工程、科学及实际应用中具有重要的意义。
然而，解析解往往难以得到，因此我们需要借助数值计算方法来求解这类问题。
本文将研究使用数值计算方法解决二维热传导方程问题，并介绍常用的数值方法及其应用。
引言：热传导方程是描述物体中温度分布的偏微分方程，通常形式为：∂u/∂t =α(∂^2u/∂x^2 + ∂^2u/∂y^2)，其中u(x, y, t)表示温度分布，α为热扩散系数。
本文将研究如何使用数值计算方法求解该方程的初始值问题。
数值方法介绍：1. 空间离散化在二维情况下，我们将区域划分为网格点，并对温度进行离散化。
常用的方法有有限差分法和有限元法。
有限差分法将二维空间离散化为矩形网格，根据差分近似导数并代入热传导方程，得到离散的方程组。
有限元法则通过将区域分解为多个小区域，利用试探函数对温度进行表示，在每个小区域内代入试探函数并求解线性方程组来得到温度分布。
2. 时间离散化对时间进行离散化也是求解二维热传导方程的重要步骤。
常用的方法有显式方法和隐式方法。
显式方法使用差分公式来逐步推进时间，从而求解温度在每个时间步长上的值。
隐式方法则利用迭代算法来求解线性方程组，通过反复迭代使得解逼近真实解。
数值方法应用与优缺点分析：1. 有限差分法有限差分法是最常用的数值方法之一，简单易于实现。
它将二维空间划分为网格点，并利用中心差分公式来近似偏导数。
在时间方向上，显式差分方法使用向前差分公式，而隐式差分方法则使用向后差分公式。
有限差分法的优点是计算效率高，在稳定性和精度上具有较好的表现，但对于非线性问题的处理稍显困难。
2. 有限元法有限元法是一种更为复杂的数值计算方法，对于复杂的边界条件和几何形状具有较好的适应性。
它将区域分解为小区域，并在每个小区域内引入试探函数。
通过求解线性方程组，可以得到温度的离散解。
基于变分原理的二维热传导方程差分格式
则 网格 上 建立 二 维热 传 导 方程 的差 分 格 式 就成 为 至
关重 要 的问 题 。 ． 尤其 对许 多应 用 性 问题 ， 激 光 如
惯性 约 束 聚变 研 究 中激 光 与靶 耦 合 过 程 ， 内爆 动 力 学压 缩 过 程 ， 都会 出 现 热 传 导 方 程 系数 和 温 度 函 数
流 通量 成 为独 立 的 未知 函 数 ， 通过 变 分原 理 ， 究 热 研 流泛 函 的极值 并 与 温 度 对 流 方 程 进 行 联 立 求 解 ； 与 此 同时 ， 出每个 网 格 内 的温 度 值 以及 通 过 网 格 四 解
周 边 界上 的热 流通量 ．
各 有 其 优缺 点 ： 分 途 径具 有 通 用 、 微 简便 等 优 点 ， 但
畸 变情 况下 数 值 求 解 结 果 的可 靠 性 ， 了对 上 述 插 除
值 公式 选 取具 有 一 定 明确 物 理含 义 和计 算 精度 的公 式外， 另一 种 有效 途径 是 本 文 将 要 讨 论 的基 于 变 分
原 理 的方 法 ． 先将 热 传导 方 程 写成 热 流形 式 ， 得 热 使
在空 间 和 时 间上 的强 烈 非 均 匀性 ， 为 问题 特 有 的 成
难点 ．
分格 式 克 服 了在非 正 交 网格 上 计算 热 流 通量 漏 失偏 大， 以及 温 度分 布 变 化 明显 依 赖 于 网 格 剖 分 形 状 等 缺点 ， 是 在计 算 公 式 中 因 涉及 到 网 格 结 点 上 的温 但 度值 ， 要 计算 每 个 网格 四周边 界 上 的热 传 导 系数 ， 需 而 这些 值 要借 助 插 值 公 式 求 得 ． 网 格 发 生 强 烈 当
［ 章 编 号 ］ １０ ．４ Ｘ（０２ ０ ．２９０ 文 ０ １ ６ ２０ ）４０ ９—６ ２
求解二维三温辐射热传导方程组的一种网格自适应方法-UCIMath
② 基于 Hessian 矩阵的网格自适应技术 ( 简记为 Hess 自适应) ,其中度量矩阵为
S = Hp .
③ 基于光子通量的自适应方法 ( 简记为 KGrad 自适应) ,其中度量矩阵 S 的表达式为
Kr S =
2 数值实验
我们通过数值实验来比较 3 种网格自适应方法在求解方程 ( 1 ) 时的效果 . 设方程的求解区域为 3 Ωxy ,0 ≤t ≤ { ( x , y , t) | ( x , y) ∈ 100} ,这里 Ωxy = ∪i = 1 Ωi ( 图 1) ,是由 3 个不同介质区域组成的一个连通但互
其中 H τ 是 H 在τ上的积分平均 . ② 存在两个正常数 β 0 和β 1 使得
∑d
i
2 τ, i
| τ|
2Πn
τ ∈ TN ; ≤β 0 , Π
其中| τ | 是 τ在新度量 Hp 下的体积 ,而 dτ, i 是τ的第 i 条边在该新度量下的长度 . 在上述关于单纯形网格的假设条件下 ,有如下定理成立 . 2 定理 1 设 u ∈C (Ω) ,网格 TN 满足假设条件 ① 和 ②, uI 是 TN 下关于 u 的线性有限元插值函数 ,则有 误差估计式
Hp = ( det H)
n
-
1 2 p+ n
记 TN 是 R 中的单纯形网格剖分 ,其中 N 是网格的单元数 . 设 TN 满足以下两个假设 : ① 存在两个正常数 α 0 和α 1 使得对任意单元 τ,
T T T α ξ ≤ξ ξ, ξ ∈ Rn , H ( x )ξ ≤α τ τ 0ξ H 1ξ H
2
2
,
Δ
Δ
Δ
由 H 2矩阵的定义知其实际上是 u 的 Hessian 矩阵 (
二维三温能量方程的差分格式及平面靶数值模拟
二维三温能量方程的差分格式及平面靶数值模拟
勇珩;段庆生;裴文兵;江松
【期刊名称】《计算物理》
【年(卷),期】2005(22)6
【摘要】通过统一二维直角坐标与柱坐标下三温能量方程的差分格式,实现了LARED-H(laser radiationelectron dynamic holehurm)程序既可以模拟正入射激光平面靶过程,也可模拟斜入射平面靶问题,从而扩展了LARED-H程序的功能.使用发展后的LARED-H程序进行了直角坐标系下线聚焦打靶耦合过程数值模拟,给出了不同入射角度下靶等离子体的膨胀过程.
【总页数】9页(P479-487)
【关键词】二维三温能量方程;平面坐标与柱坐标;平面靶数值模拟
【作者】勇珩;段庆生;裴文兵;江松
【作者单位】北京应用物理与计算数学研究所
【正文语种】中文
【中图分类】O532
【相关文献】
1.二维半线性伪抛物方程差分格式及数值模拟 [J], 涂慧;刘超;江成顺
2.三维椭圆方程基于微分矩阵的差分格式及数值模拟 [J], 李艳琴
3.二维三温能量方程的 RT 混合线性有限元及其数值模拟 [J], 陈铁军;李跃军;肖建新
4.五对角紧致差分格式优化及二维声波传播波动方程数值模拟 [J], 汪勇;穆鹏飞;蔡文杰;王鹏;桂志先
5.二维三温能量方程的九点差分格式及其迭代解法 [J], 符尚武;付汉清;沈隆钧;黄书科;陈光南
因版权原因，仅展示原文概要，查看原文内容请购买。
CP090-计算物理热传导方程的差分解法
t T
（9.6）
9.2 一维热传导方程的差分解法
边界条件： 3、 第三类边界条件：
u (, t ) x (t )u (, t ) g (t ) u (l , t ) (t )u (l , t ) g (t ) x
其中， (t )
9.1 热传导方程概述
四、三维非齐次热传导方程 由能量守恒定律，即
Q Q Q
可得：
或
u t dt V [c t ( Ku) F ]dV u c ( Ku ) F ( x, y, z, t ) t
t
（9.1）
式（9.1）称为各向同性介质有热源的热传导方程，也叫三维 非齐次热传导方程。
9.1 热传导方程概述
五、三维齐次热传导方程 当介质均匀（ c 、 和 K 为常数） 内无热源（ F ( x, y, z, t ) ）时： 、V
u c Ku, t
上式可表示为：
其中 x y z
u K u u u ( ) t c x y z
9.2 一维热传导方程的差分解法
例 9.1 求热传导方程混合问题：
u u t x u ( x,) x( x) u (, t ) , u (, t )
x , t x t
的数值解，取 N=10，h=0.1,计算到 k=36 为止。
9.2 一维热传导方程的差分解法
2、差分格式的稳定条件
h
3、具体步骤 （1）给定 , l , h, , T ； （2）计算 h / , N l / h, M T / ；
热传导方程的差分格式讲解
A(i,i-1)=-r;
end
end
u=zeros(N+1,M+1);
u(N+1,:)=u1;
fork=1:N
b=u(N+2-k,2:M);
u(N+1-k,2:M)=inv(A)*b';
end
uT=u(1,:);
x1=[0,x,1];
plot(x1,uT,'o')
hold
u_xt=exp(-pi*pi*T)*sin(pi.*x1);
u(N+1-k,2:M)=inv(A)*b;
end
uT=u(1,:);
x1=[0,x,1];
plot(x1,uT,'o')
hold
u_xt=exp(-pi*pi*T)*sin(pi.*x1);
plot(x1,u_xt,'r')
e=u_xt-uT;
plot(x1,u_xt,'r')
e=u_xt-uT;
（3）六点对称格式：
源程序：
function[e]=six(dx,dt,T)
M=1/dx;
N=T/dt;
u1=zeros(1,M+1);
x=[1:M-1]*dx;
u1([2:M])=sin(pi*x);
r=dt/dx/dx;r2=2+2*r;r3=2-2*r;
此格式为显格式。
其矩阵表达式如下：
（2）向后差分格式
向后差分格式，即
（2）
其中 （2）式可改写成
此种差分格式被称为隐格式。
其矩阵表达式如下：
（3）六点对称格式
六点差分格式：
二维热传导方程的三层显式差分格式
二维热传导方程的三层显式差分格式
刘继军
【期刊名称】《应用数学和力学》
【年(卷),期】2003(24)5
【摘要】对二维热传导方程构造了一个稳定的三层显式差分格式求其数值解,其背景源于高维热力学反问题迭代算法中对正问题小计算量算法的需求·首先建立一个含参数的一般差分格式去逼近微分方程,并得到了最优截断误差·然后导出了参数应满足的条件以保证差分格式的稳定性·最后给出了数值的例子并和其它算法进行比较。
【总页数】8页(P537-544)
【关键词】热传导方程;差分格式;误差估计;稳定性
【作者】刘继军
【作者单位】东南大学数学系
【正文语种】中文
【中图分类】O175.26
【相关文献】
1.解三维热传导方程的一族对称含参数高精度的显式差分格式 [J], 孙鸿烈
2.关于二维热传导方程的一族显式差分格式 [J], 金承日;刘家琦
3.二维热传导方程的一个新的两层显式格式 [J], 马小宁
4.解高维热传导方程的一族高精度显式差分格式 [J], 陈贞忠;马小霞
5.解高维热传导方程的一族高精度的显式差分格式 [J], 孙鸿烈
因版权原因，仅展示原文概要，查看原文内容请购买。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
n
28
T
ej eP*L
8 M 27 L / 4 < De Ω = Ωxy × [0, T ], Ωxy = (x, y) | 0 ≤ x ≤ 1; 0 ≤ y ≤ 1 , nBOpR)Y < ?^0M T $
Tα = 0, α = e, i, r, (x, y, t) ∈ ∂ Ωxy × [0, T ]. (1.5)
3
^0e
0 Tα (x, y, 0) = Tα (x, y ),
α = e, i, r, ωei
(1.4) ωer
∗
Y Te , Ti, Tr AiS L n ?Vdk] Zdk] dk] ?VdCZd VdC dM}g4 ~L ρ e@XMt] p 2001 6 I 4 #EK 2003 5 I 14 #EK uz )e P A (19871043 ) . 7 e P A (Q98A07115) c` y
1
_ '
,
1
K GKJ (ICF) MbfLUw;Y **(p < UVd 8G Z d 8G d 8Gj/Mq2n a j/h [1], YA@!_}gT-4 M (k 8Gj/hM <M@](MjV r0wP [2,3] X4j}B7;\bfV{ X"a{MHT'n{< n^:nL"%H [4−6], 2`)~<bf(k 8Gj /hMnL"@ 'n{< =\r/n}j'n{<MCb qX LJk =7MZ+ H 1 hLz' &<!qX ^:nL"je Q, bftY e) ftY < -B9 &R Xx&R q tYML w A;Ip m bnj/hM3 UtY
n+1 Teh,ij − Teh,ij n+ 1 2
n 2 Teh,ij − Teh,ij
n+ 1
(2.4a)
n+1 Cve,ij
∆t
=
1 n+1 n n δy (ke δy (Teh − Teh ))ij , ρ
1 ≤ j ≤ N − 1. (2.4b)
Zd 8Gj/ (1.2) MnL"4<'n{<e
Cve 1 ∂Te = div (K (ρ, Te ) grad Te ) + ωei (Ti − Te ) + ωer (Tr − Te ), ∂t ρ 1 ∂Ti Cvi = div (K (ρ, Ti ) grad Ti ) − ωei (Ti − Te ), ∂t ρ 1 ∂Tr = div (K (ρ, Tr ) grad Tr ) − ωer (Tr − Te ), Cvr ∂t ρ (1.1) (1.2) (1.3)
S 27 L S 1 1 H 2004
9<
N
ACTA MATHEMATICAE APPLICATAE SINICA
Vol. 27 No. 1 Jan., 2004
}
;
; 876 532 9 4 Z '@ F
r u
(
wy~v
z
∗
O
266071)
o
o 5* =cg)l 9Hk0iNoM#5=A (o|= ="~hkf 5(o|=Nl>8 , H iM %{("rY x )l 9Hk0i oM#(o|= Dc rY
2
4
2
n+1 n+1 n+1 n+1 n+1 n n n δx (kα δx Tαh )ij = h−2 kα,i + 1 ,j (Tαh,i+1,j − Tαh,ij ) − kα,i− 1 ,j (Tαh,ij − Tαh,i−1,j ) , n+1 n δy (kα δy Tαh )ij
(2.1) (2.2) (2.3)
=h
−2
n+1 n kα,i,j (Tαh,i,j 43;1 Tαh,ij )
−
n+1 n kα,i,j (Tαh,ij −1 2
−
n+1 Tαh,i,j −1 )
,
∇h ·
n+1 n (kα ∇h Tαh )ij
=
n+1 n δx (kα δx Tαh )ij
+
n+1 n δy (kα δy Tαh )ij .
n+1 Cvi,ij
∆t
=
1 n+1 n n δy (ki δy (Tih − Tih ))ij , ρ
1
cg)l 9Hk0iNoM#5=(o|= d 8Gj/ (1.3) MnL"4<'n{<e
n+1 Cvr,ij n 2 Trh,ij − Trh,ij n+ 1
J-Q
29
∆t 1 n+1 n n n n + δy (kr δy Trh )ij − ωer,ij (Trh − Teh )ij , ρ
Vd 8Gj/ (1.1) MnL"4<'n{<e
n+1 Cve,ij
1 1 n+ 1 n n n δx (ke δx Teh 2 )ij + δy (ke δy Teh )ij ∆t ρ ρ n+1 n+1 n n n n + ωei,ij (Tih − Teh )ij + ωer,ij (Trh − Teh )ij , 1 ≤ i ≤ N − 1, =
n%;
+4
0 < C∗ ≤ Cvα ≤ C ∗ , ∂K (ρ, Tα ) ≤ D∗ , ∂Tα 0 < K∗ ≤ K (ρ, Tα ) ≤ K ∗ , α = e, i, r. (1.6) (1.7)
7 C ∗ , C∗ , K ∗ , K∗ , D∗ MeP*L +4tYM<!j/M~LK@BEM 2 qt 1 n n n , xi = ih, yj = jh; ∆t = T mh= N L , t = n∆t, Wij = W (xi , yj , t ). : α e, i, r YM!.)~ : Tαh ? Tα M'nDQ (
n n n kα,i + 1 ,j = K (ρ, Tαh,i+1,j ) + K (ρ, Tαh,ij ) 2, n n n Kα,i + 1 ,j = K (ρ, Tα,i+1,j ) + K (ρ, Tα,ij ) 2,
2 2
?
(
n n kα,i,j , Kα,i,j +1 +1
2 2
YQX/
n+1 Cvi,ij n 2 Tih,ij − Tih,ij n+ 1
∆t 1 n+1 n n n n + δy (ki δy Tih )ij − ωei,ij (Tih − Teh )ij , ρ
n+1 2 Tih,ij − Tih,ij n+ 1
=
1 n+ 1 n δx (ki δx Tih 2 )ij ρ 1 ≤ i ≤ N − 1, 1 ≤ j ≤ N − 1. (2.5b) (2.5a)