一维热传导方程的数值解
合集下载
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
Ａｂｓｔｒａｃｔ ： Ｕｓｅｉｎｇ ｔｈｅ ｄｉｆｆｅｒｅｎｃｅ ｔｏ ｇｏ ｏｎ ｎｕｍｅｒｉｃａｌ ｃａｌｃｕｌａｔｉｏｎ ｆｏｒ ｏｎｅ－ｄｉｍｅｎｓｉｏｎａｌ ｈｅａｔ ｃｏｎｄｕｃｔｉｏｎ ｅｑｕａｔｉｏｎ ，ａｎｄ ｇｉｖｅ ａｎ ｉｌｌｕｓｔｒａｔｉｏｎ ｏｆ ｔｈｅ ｓｏｌｕｔｉｏｎ ． Ｋｅｙ ｗｏｒｄｓ ： ｄｉｆｆｅｒｅｎｃｅ ｍｅｔｈｏｄ ；ｏｎｅ－ｄｉｍｅｎｓｉｏｎａｌ ｈｅａｔ ｃｏｎｄｕｃｔｉｏｎ ｅｑｕａｔｉｏｎ ；ｂｏｕｎｄａｒｙ ｃｏｎｄｉｔｉｏｎ
ｉ ，ｊ
Δｔ
（５）
以（４） 、（５） 代入（１） 式得
ｕ － ｕ ｉ ，ｊ＋ １
ｉ ，ｊ
Δｔ
＝
ａ２
ｕｉ －１ ，ｊ
－ ２ ｕｉ ，ｊ Δｘ２
＋
ｕｉ＋ １ ，ｊ
＋
ｆ （ ｉ ，ｊ）
（６）
解得
ｕｉ ，ｊ＋ １ ＝ ｃ （ ｕｉ －１ ，ｊ ＋ ｕｉ＋ １ ，ｊ ） ＋ （１ － ２ ｃ） ｕｉ ，ｊ ＋ Δｔｆ （ ｉ ，ｊ ）
式处理即可 ．
２） 如果是波动方程 ，那么初始条件中含有“初始速度” ，即 ｕｔ （ｘ ，ｔ） ｜ ｔ ＝ ０ ＝ 宝（ｘ） ，利用向后差分法得
抄 ｕ（ ｉ ，ｌ） 抄ｔ
＝
ｕｉ ，２ － ｕｉ ，０ ２ Δｔ
＝
宝（ ｘ ）
（３５）
则
ｕｉ ，０ ＝ ｕｉ ，２ － ２ Δｔ宝（ｘ ）
（３６）
利用（３６） 式便得以求解 ．
用中心差分近似代替对空间的偏微分即2u抄用向前差分近似代替对时间的偏微分即uij1uij2uxxfxt0xl0t123抄2xui1j2uijui1jx24抄u抄tt5以45代入1式得uij1uijta2ui1j2uijui1jx2fij6解得uij1cui1jui1j12cuijtfij7其中ctax228根据式7如果已知j不同i坐标每一个格点的温度值并且由11类边界条件可知两边界i1及in上的温度值那么就可以求出j1坐标上每一个格点上的温度值
次边界条件的原因 ．在图 ９ 中 ｘ ＝ ｌ 端的变化曲线开始时变化缓慢 ，而在中间 ｘ ＝ ０ ．５ｍ 的附近温度上升
较快 ，两边相对缓慢 ．比较可知 ，例 ３ 和例 ４ 导热细杆中的温度变化规律左右是对调的 ．
３ 结语
１） 如果数理方程为 ２２ 类边界条件 ，那么同时将 ｉ ＝ １ 和 ｉ ＝ Ｎ 两个边界点进行类似于（１８） 及（２４）
２１４
淮阴师范学院学报 （自然科学版 ）
第３卷
齐次边界条件的原因 ．图 ７ 中 ｘ ＝ ０ 的变化曲线开始时变化缓慢 ，而在 ｘ ＝ ０ ．５ｍ 的附近温度上升较快 ，
两边相对缓慢 ．
例 ４ 设一维热传导问题为
ｕｔ
＝
ａ２ ｕｘｘ
＋
ｇｓｉｎ
πｘ ｌ
（０
＜
ｘ＜
ｌ） （０ ＜
ｔ＜
∞）
（３２）
ｕ（０ ，ｔ） ＝ ０ ， ｕｘ （ ｌ ，ｔ） ＝ ０ （０ ≤ ｔ ≤ ∞ ）
图 １ 第一类边界条件下热传导方程的图解
例 １ 设定解问题为
ｕｔ
＝
ａ２ ｕｘｘ
＋
Ａｓｉｎ
πｘ ｌ
（０
＜
ｘ＜
ｌ） （０ ＜
ｔ＜
∞）
ｕ（０ ，ｔ） ＝ ０ ， ｕ（ ｌ ，ｔ） ＝ ０ （０ ≤ ｔ ≤ ∞ ）
ｕ（ｘ ，０） ＝ ０ （０ ≤ ｘ ≤ ｌ）
这是一个两端温度为 ０ ℃ ，并且具有热源的定解问题 ．本定解问题有解析解 ，其解为
１ 第一类边界条件的一维热传导混合问题的解法
设含第一类边界条件的一维热传导混合问题为 ：
ｕｔ ＝ ａ２ ｕｘｘ ＋ ｆ （ｘ ，ｔ） （０ ＜ ｘ ＜ ｌ） （０ ＜ ｔ ＜ ∞ ）
（１）
ｕ（０ ，ｔ） ＝ μ１ （ ｔ） ， ｕ（ ｌ ，ｔ） ＝ μ２ （ ｔ） （０ ≤ ｔ ≤ ∞ ）
ｕ（ｘ ，０） ＝ ｂｘ（ ｌ － ｘ）／ｌ２ （０ ≤ ｘ ≤ ｌ）
例 ２ 的数值计算结果如图 ３ 所示 ．
（１０） （１１） （１２）
（１３）
（１４） （１５） （１６）
图 ２ 例 １ 中相同 ｔ 不同 ｘ 的温度变化曲线 图 ３ 例 ２ 中相同 ｔ 不同 ｘ 的温度变化曲线
等份 ．令 ｉ 表示位置 ｘ 横轴 ，ｊ 表示时间 ｔ 纵轴 ．网格上每个格点对应一个温度值 ．用中心差分近似代替对 空间的偏微分 ，即
抄２ 抄２
ｕ ｘ
＝
ｕｉ －１ ，ｊ － ２ ｕｉ ，ｊ ＋ ｕｉ＋ １ ，ｊ Δｘ２
（４）
用向前差分近似代替对时间的偏微分 ，即
抄ｕ 抄ｔ
＝
ｕ － ｕ ｉ ，ｊ＋ １
及 ｉ ＝ Ｎ 上的温度值 ，那么就可以求出 ｊ ＋ １ 坐标上每一个格点上的温度值 ．因此 ，利用（７） 式从初始条
件 ｊ ＝ １ 开始 ，就可逐步算出每一个格点上的温度值 ，运算过程如图 １ 所示 ．
这里必须特别指出的是算法的稳定性问题 ，即解达到稳定的条件是
ｃ＝
Δ ｔａ２ Δｘ２
≤
１ ２
（９）
时刻的温度值 ．因 ｘ ＝ ｌ 处对应 ｉ ＝ Ｎ ＋ １ ，则由（１８） 式得
抄 ｕＮ＋ １ ，ｊ 抄ｘ
＝
－ ｕＮ ，ｊ ＋ ｕＮ＋ ２ ，ｊ ２ Δｘ
＝
μ（ｊ）
（２４）
所以
ｕＮ＋ ２ ，ｊ ＝ ｕ（Ｎ ，ｊ ） ＋ ２ Δｘμ（ ｊ）
（２５）
令（７） 式中 ｉ ＝ Ｎ ＋ １ 得
ｕ１ ，ｊ＋ １ ＝ ｃ（ ｕ０ ，ｊ ＋ ｕ２ ，ｊ ） ＋ （１ － ２ ｃ ） ｕ１ ，ｊ ＋ Δｔｆ （Ｎ ＋ １ ，ｊ ）
参考文献 ：
［１］ 佩卡姆 Ｈ Ｄ ．ＢＡＳＩＣ 解大学物理题 ［Ｍ］ ．北京 ：电子工业出版社 ，１９８５ ． ［２］ 沙瓦陀里 Ｍ Ｇ ．工程学中的数值方法 ［Ｍ］ ．北京 ：高等教育出版社 ，１９５９ ． ［３］ 梁昆淼 ．数学物理方法 ［Ｍ］ ．北京 ：高等教育出版社 ，１９９５ ． ［４］ 徐效海 ．数学物理方法引论 ［Ｍ］ ．南京 ：南京大学出版社 ，１９９９ ．
（７）
其中
ｃ＝
Δ ｔａ２ Δｘ２
（８）
收稿日期 ：２００４－０２－２１ 作者简介 ：徐建良 （１９６４－） ，男 ，江苏武进人 ，讲师 ，主要从事物理教学研究 ．
第３期
徐建良等 ：一维热传导方程的数值解
２１１
根据式（７） ，如果已知 ｊ（不同 ｉ） 坐标每一个格点的温度值 ，并且由 １１ 类边界条件可知两边界 ｉ ＝ １
（２６）
由（２５） 、（２６） 式解得
第３期
徐建良等 ：一维热传导方程的数值解
２１３
ｕＮ＋ １ ，ｊ＋ １ ＝ ２ ｃ（ ｕＮ ，ｊ ＋ Δｘμ（ ｊ）） ＋ （１ － ２ ｃ ） ｕＮ＋ １ ，ｊ ＋ Δｔｆ （Ｎ ＋ １ ，ｊ ）
（２７）
再令（７） 中 ｉ ＝ Ｎ 则得
ｕＮ ，ｊ＋ １ ＝ ｃ （ ｕＮ －１ ，ｊ ＋ ｕＮ＋ １ ，ｊ ） ＋ （１ － ２ ｃ ） ｕ２ ，ｊ ＋ Δｔｆ （Ｎ ，ｊ）
（２１）
再令（７） 中 ｉ ＝ ２ 则得
ｕ２ ，ｊ＋ １ ＝ α（ ｕ１ ，ｊ ＋ ｕ３ ，ｊ ） ＋ （１ － ２ α） ｕ２ ，ｊ ＋ Δｔｆ （２ ，ｊ ）
（２２）
这样就可先由（２２） 式算出 ｉ ＝ ２ 上的各点温度值 ，再由（２１） 式算出边界 ｘ ＝ ０ ，即 ｉ ＝ １ 上各点的温度
值 ．（２１） 和（２２） 式的运算过程如图 ４ ．
摘 要 ：利用差分法对数理方程的多个较复杂的一维热传导问题进行分析 ，并进行数值计算 ， 给出了直观的图象 ． 关键词 ：差分法 ；一维热传导方程 ；边界条件 中图分类号 ：Ｏ５５１ ．３ 文献标识码 ：Ａ 文章编号 ：１６７１－６８７６（２００４）０３－０２１０－０５
数理方程中混合问题的解法通常是比较复杂的 ，有时很难有解析解 ，即便解出 ，其解的形式也通常 是一个无穷级数形式［１ ４］ ．本文以一维热传导混合问题为例 ，探讨利用差分法对数理方程进行数值解的 方法 ．
（３３）
ｕ（ｘ ，０） ＝ １００ ℃ ｘ ／ｌ （０ ≤ ｔ ≤ ｌ）
（３４）
图 ８ 中 ｔ ＝ ０ 是初始状态 ，为一直线 ，右端为 ｘ ＝ ｌ ，每一条曲线左边的开始端平坦 ，这正是由于在 ｘ ＝ ｌ
处为第二类齐次边界条件的原因 ．而右端 ｘ ＝ ０ 的温度始终为 ０ ，这是由于该点的边界条件为第一类齐
图 ４ ２１ 类边界条件下的运行图 图 ５ １２ 类边界条件下的运行图
２ ．２ １２ 类边界条件的处理办法
设 ｘ ＝ ｌ 端满足第二类边界条件 ，即
ｕｘ （ ｘ ，ｔ） ｜ ｘ ＝ ｌ ＝ μ（ ｔ）
（２３）
那么也不能直接利用（７） 式 ，必须首先利用（２３） 式及初始条件（３） 逐步求出边界 ｘ ＝ ｌ 及其他各点处各
２１２
淮阴师范学院学报 （自然科学版 ）
第３卷
从结果可以知道热传导杆两端的温度始终保持为 ０ ，杆的中点温度总是高于其它点温度 ，各点温度 随着时间变化逐渐降低 ．
２ 含第二类边界条件的一维热传导混合问题的数值解法
２ ．１ ２１ 类边界条件的处理办法
设 ｘ ＝ ０ 端满足第二类边界条件 ，即
第 ３ 卷第 ３ 期 ２００４ 年 ８ 月
淮阴师范学院学报 （自然科学版 ）
ＪＯＵＲＮＡＬ ＯＦ ＨＵＡＩＹＩＮ ＴＥＡＣＨＥＲＳ ＣＯＬＬＥＧＥ （ＮＡＴＵＲＡＬ ＳＣＩＥＮＣＥ ＥＤＩＴＩＯＮ）
Ｖｏｌ畅３ Ｎｏ畅３ Ｍａｙ ．２００４
一维热传导方程的数值解
徐建良 ，汤炳书
（连云港高等师范专科学校 物理系 ，江苏 连云港 ２２２００６）
（２）
ｕ（ｘ ，０） ＝ φ（ｘ） （０ ≤ ｘ ≤ ｌ）
（３）
这是一个 １１ 类非齐次边界条件的一维热传导问题 ，通常这一类混合问题是很难解的 ，即便解出 ，其解也
通常是一个无穷级数的形式 ，对该解的物理意义不能直接讨论 ，不能给出直观的图象 ．
１ ．１ 计算方法
为求解方程（１） ，首先定义函数 ｕ（ｘ ，ｔ）的时间与空间的网格 ，将 ｘ 坐标分成 Ｎ 等份 ，将 ｔ 坐标分成 Ｍ
（２８）
这样就可先由（２８） 式算出 ｉ ＝ Ｎ 上的各点温度值 ，再由（２７） 式算出边界 ｘ ＝ ｌ ，即 ｉ ＝ Ｎ ＋ １ 上各点的
温度值 ．（２７） 和（２８） 式运算过程如图 ５ ．
图 ６ 例 ３ 中相同 ｔ 不同 ｘ 的温度变化曲线 图 ７ 例 ３ 中相同 ｘ 不同 ｔ 的温度变化曲线
＝
μ（ｊ）
（１８）
所以
ｕ０ ，ｊ ＝ ｕ２ ，ｊ － ２ Δｘμ（ ｊ ）
（１９）
令（７） 式中 ｉ ＝ １ 得
ｕ１ ，ｊ＋ １ ＝ α（ ｕ０ ，ｊ ＋ ｕ２ ，ｊ ） ＋ （１ － ２ α） ｕ１ ，ｊ ＋ Δｔｆ （１ ，ｊ ）
（２０）
由（１９） 、（２０） 式解得
ｕ１ ，ｊ＋ １ ＝ ２ α（ ｕ２ ，ｊ － Δｘμ（ ｊ ）） ＋ （１ － ２ α） ｕ１ ，ｊ ＋ Δｔｆ（１ ，ｊ）
Ｎｕｍｅｒｉｃａｌ Ｃａｌｃｕｌａｔｉｏｎ ｏｆ Ｏｎｅ－ｄｉｍｅｎｓｉｏｎａｌ Ｈｅａｔ Ｃｏｎｄｕｃｔｉｏｎ Ｅｑｕａｔｉｏｎ
ＸＵ Ｊｉａｎ－ｌｉａｎｇ ，ＴＡＮＧ Ｂｉｎｇ－ｓｈｕ
（Ｄｅｐａｒｔｍｅｎｔ ｏｆ Ｐｈｙｓｉｃｓ ，Ｌｉａｎｙｕｎｇａｎｇ Ｔｅａｃｔｈｅｒｓ Ｃｏｌｌｅｇｅ ，Ｌｉａｎｙｕｎｇａｎｇ ２２２００６ ，Ｃｈｉｎａ）
图 ８ 例 ４ 中相同 ｔ 不同 ｘ 的温度变化曲线 图 ９ 例 ４ 中相同 ｘ 不同 ｔ 的温度变化曲线
图 ６ 中 ｔ ＝ ０ 是初始状态 ，为一直线 ，左端为 ｘ ＝ ０ ，每一条曲线左边的开始端平坦 ，这正是由于在 ｘ ＝ ０ 处为第二类齐次边界条件的原因 ．而右端 ｘ ＝ ｌ 的温度始终为０ ，这是由于该点的边界条件为第一类
ｕ（ｘ ，ｔ） ＝
Ａ π２
ｌ２ ａ２
（１
－
ｅ － π２ ａ２ ｔ／ ｌ２ ）ｓｉｎ
πｘ ｌ
由数值计算得出的结果如图 ２ 所示 ．
例 ２ 设定解问题为
ｕｔ － ａ２ ｕｘｘ ＝ ０ （０ ＜ ｘ ＜ ｌ） （０ ＜ ｔ ＜ ∞ ）
ｕ（０ ，ｔ） ＝ ０ ， ｕ（ ｌ ，ｔ） ＝ ０ （０ ≤ ｔ ≤ ∞ ）
例 ３ 设一维热传导问题为
ｕｔ
＝
ａ２ ｕｘｘ
＋
ｇｓｉｎ
πｘ ｌ
（０
＜
ｘ＜
ｌ） （０ ＜
ｔ＜
Hale Waihona Puke ∞）（２９）ｕｘ （０ ，ｔ） ＝ ０ ， ｕ（ ｌ ，ｔ） ＝ ０ （０ ≤ ｔ ≤ ∞ ）
（３０）
ｕ（ｘ ，０） ＝ １００ ℃ （ ｌ － ｘ）／ｌ （０ ≤ ｘ ≤ ｌ）
（３１）
这是一个 ｘ ＝ ｌ 端为第二类齐次边界条件的并且具有热源的热传导混合问题 ，即 μ（ ｔ） ＝ ０ ．在编程时 ， 温度采用国际温标 ．设 ｌ ＝ １ｍ ，时间范围为 ０ ～ １ ｓ ，ｇ ＝ ２７３０ ｋ／ｓ ，将 ｌ 分割为 １００ 份 ，时间分为 ５００ 份 ． 时空网格有 １０１ × ５０１ 个格点 ．例 ３ 的数值计算结果如图 ６ 、图 ７ 所示 ．
ｕｘ （ ｘ ，ｔ） ｜ ｘ ＝ ０ ＝ μ（ ｔ）
（１７）
则此时就不能直接利用（７） 式求解 ，必须首先利用（１７） 式及初始条件（３） 逐步求出边界 ｘ ＝ ０ 及其
他各点处各时刻的温度值 ．因 ｘ ＝ ０ 处对应 ｉ ＝ １ ，则由（１７） 式得
抄 ｕ１ ，ｊ 抄ｘ
＝
－ ｕ０ ，ｊ ＋ ｕ２ ，ｊ ２ Δｘ
ｉ ，ｊ
Δｔ
（５）
以（４） 、（５） 代入（１） 式得
ｕ － ｕ ｉ ，ｊ＋ １
ｉ ，ｊ
Δｔ
＝
ａ２
ｕｉ －１ ，ｊ
－ ２ ｕｉ ，ｊ Δｘ２
＋
ｕｉ＋ １ ，ｊ
＋
ｆ （ ｉ ，ｊ）
（６）
解得
ｕｉ ，ｊ＋ １ ＝ ｃ （ ｕｉ －１ ，ｊ ＋ ｕｉ＋ １ ，ｊ ） ＋ （１ － ２ ｃ） ｕｉ ，ｊ ＋ Δｔｆ （ ｉ ，ｊ ）
式处理即可 ．
２） 如果是波动方程 ，那么初始条件中含有“初始速度” ，即 ｕｔ （ｘ ，ｔ） ｜ ｔ ＝ ０ ＝ 宝（ｘ） ，利用向后差分法得
抄 ｕ（ ｉ ，ｌ） 抄ｔ
＝
ｕｉ ，２ － ｕｉ ，０ ２ Δｔ
＝
宝（ ｘ ）
（３５）
则
ｕｉ ，０ ＝ ｕｉ ，２ － ２ Δｔ宝（ｘ ）
（３６）
利用（３６） 式便得以求解 ．
用中心差分近似代替对空间的偏微分即2u抄用向前差分近似代替对时间的偏微分即uij1uij2uxxfxt0xl0t123抄2xui1j2uijui1jx24抄u抄tt5以45代入1式得uij1uijta2ui1j2uijui1jx2fij6解得uij1cui1jui1j12cuijtfij7其中ctax228根据式7如果已知j不同i坐标每一个格点的温度值并且由11类边界条件可知两边界i1及in上的温度值那么就可以求出j1坐标上每一个格点上的温度值
次边界条件的原因 ．在图 ９ 中 ｘ ＝ ｌ 端的变化曲线开始时变化缓慢 ，而在中间 ｘ ＝ ０ ．５ｍ 的附近温度上升
较快 ，两边相对缓慢 ．比较可知 ，例 ３ 和例 ４ 导热细杆中的温度变化规律左右是对调的 ．
３ 结语
１） 如果数理方程为 ２２ 类边界条件 ，那么同时将 ｉ ＝ １ 和 ｉ ＝ Ｎ 两个边界点进行类似于（１８） 及（２４）
２１４
淮阴师范学院学报 （自然科学版 ）
第３卷
齐次边界条件的原因 ．图 ７ 中 ｘ ＝ ０ 的变化曲线开始时变化缓慢 ，而在 ｘ ＝ ０ ．５ｍ 的附近温度上升较快 ，
两边相对缓慢 ．
例 ４ 设一维热传导问题为
ｕｔ
＝
ａ２ ｕｘｘ
＋
ｇｓｉｎ
πｘ ｌ
（０
＜
ｘ＜
ｌ） （０ ＜
ｔ＜
∞）
（３２）
ｕ（０ ，ｔ） ＝ ０ ， ｕｘ （ ｌ ，ｔ） ＝ ０ （０ ≤ ｔ ≤ ∞ ）
图 １ 第一类边界条件下热传导方程的图解
例 １ 设定解问题为
ｕｔ
＝
ａ２ ｕｘｘ
＋
Ａｓｉｎ
πｘ ｌ
（０
＜
ｘ＜
ｌ） （０ ＜
ｔ＜
∞）
ｕ（０ ，ｔ） ＝ ０ ， ｕ（ ｌ ，ｔ） ＝ ０ （０ ≤ ｔ ≤ ∞ ）
ｕ（ｘ ，０） ＝ ０ （０ ≤ ｘ ≤ ｌ）
这是一个两端温度为 ０ ℃ ，并且具有热源的定解问题 ．本定解问题有解析解 ，其解为
１ 第一类边界条件的一维热传导混合问题的解法
设含第一类边界条件的一维热传导混合问题为 ：
ｕｔ ＝ ａ２ ｕｘｘ ＋ ｆ （ｘ ，ｔ） （０ ＜ ｘ ＜ ｌ） （０ ＜ ｔ ＜ ∞ ）
（１）
ｕ（０ ，ｔ） ＝ μ１ （ ｔ） ， ｕ（ ｌ ，ｔ） ＝ μ２ （ ｔ） （０ ≤ ｔ ≤ ∞ ）
ｕ（ｘ ，０） ＝ ｂｘ（ ｌ － ｘ）／ｌ２ （０ ≤ ｘ ≤ ｌ）
例 ２ 的数值计算结果如图 ３ 所示 ．
（１０） （１１） （１２）
（１３）
（１４） （１５） （１６）
图 ２ 例 １ 中相同 ｔ 不同 ｘ 的温度变化曲线 图 ３ 例 ２ 中相同 ｔ 不同 ｘ 的温度变化曲线
等份 ．令 ｉ 表示位置 ｘ 横轴 ，ｊ 表示时间 ｔ 纵轴 ．网格上每个格点对应一个温度值 ．用中心差分近似代替对 空间的偏微分 ，即
抄２ 抄２
ｕ ｘ
＝
ｕｉ －１ ，ｊ － ２ ｕｉ ，ｊ ＋ ｕｉ＋ １ ，ｊ Δｘ２
（４）
用向前差分近似代替对时间的偏微分 ，即
抄ｕ 抄ｔ
＝
ｕ － ｕ ｉ ，ｊ＋ １
及 ｉ ＝ Ｎ 上的温度值 ，那么就可以求出 ｊ ＋ １ 坐标上每一个格点上的温度值 ．因此 ，利用（７） 式从初始条
件 ｊ ＝ １ 开始 ，就可逐步算出每一个格点上的温度值 ，运算过程如图 １ 所示 ．
这里必须特别指出的是算法的稳定性问题 ，即解达到稳定的条件是
ｃ＝
Δ ｔａ２ Δｘ２
≤
１ ２
（９）
时刻的温度值 ．因 ｘ ＝ ｌ 处对应 ｉ ＝ Ｎ ＋ １ ，则由（１８） 式得
抄 ｕＮ＋ １ ，ｊ 抄ｘ
＝
－ ｕＮ ，ｊ ＋ ｕＮ＋ ２ ，ｊ ２ Δｘ
＝
μ（ｊ）
（２４）
所以
ｕＮ＋ ２ ，ｊ ＝ ｕ（Ｎ ，ｊ ） ＋ ２ Δｘμ（ ｊ）
（２５）
令（７） 式中 ｉ ＝ Ｎ ＋ １ 得
ｕ１ ，ｊ＋ １ ＝ ｃ（ ｕ０ ，ｊ ＋ ｕ２ ，ｊ ） ＋ （１ － ２ ｃ ） ｕ１ ，ｊ ＋ Δｔｆ （Ｎ ＋ １ ，ｊ ）
参考文献 ：
［１］ 佩卡姆 Ｈ Ｄ ．ＢＡＳＩＣ 解大学物理题 ［Ｍ］ ．北京 ：电子工业出版社 ，１９８５ ． ［２］ 沙瓦陀里 Ｍ Ｇ ．工程学中的数值方法 ［Ｍ］ ．北京 ：高等教育出版社 ，１９５９ ． ［３］ 梁昆淼 ．数学物理方法 ［Ｍ］ ．北京 ：高等教育出版社 ，１９９５ ． ［４］ 徐效海 ．数学物理方法引论 ［Ｍ］ ．南京 ：南京大学出版社 ，１９９９ ．
（７）
其中
ｃ＝
Δ ｔａ２ Δｘ２
（８）
收稿日期 ：２００４－０２－２１ 作者简介 ：徐建良 （１９６４－） ，男 ，江苏武进人 ，讲师 ，主要从事物理教学研究 ．
第３期
徐建良等 ：一维热传导方程的数值解
２１１
根据式（７） ，如果已知 ｊ（不同 ｉ） 坐标每一个格点的温度值 ，并且由 １１ 类边界条件可知两边界 ｉ ＝ １
（２６）
由（２５） 、（２６） 式解得
第３期
徐建良等 ：一维热传导方程的数值解
２１３
ｕＮ＋ １ ，ｊ＋ １ ＝ ２ ｃ（ ｕＮ ，ｊ ＋ Δｘμ（ ｊ）） ＋ （１ － ２ ｃ ） ｕＮ＋ １ ，ｊ ＋ Δｔｆ （Ｎ ＋ １ ，ｊ ）
（２７）
再令（７） 中 ｉ ＝ Ｎ 则得
ｕＮ ，ｊ＋ １ ＝ ｃ （ ｕＮ －１ ，ｊ ＋ ｕＮ＋ １ ，ｊ ） ＋ （１ － ２ ｃ ） ｕ２ ，ｊ ＋ Δｔｆ （Ｎ ，ｊ）
（２１）
再令（７） 中 ｉ ＝ ２ 则得
ｕ２ ，ｊ＋ １ ＝ α（ ｕ１ ，ｊ ＋ ｕ３ ，ｊ ） ＋ （１ － ２ α） ｕ２ ，ｊ ＋ Δｔｆ （２ ，ｊ ）
（２２）
这样就可先由（２２） 式算出 ｉ ＝ ２ 上的各点温度值 ，再由（２１） 式算出边界 ｘ ＝ ０ ，即 ｉ ＝ １ 上各点的温度
值 ．（２１） 和（２２） 式的运算过程如图 ４ ．
摘 要 ：利用差分法对数理方程的多个较复杂的一维热传导问题进行分析 ，并进行数值计算 ， 给出了直观的图象 ． 关键词 ：差分法 ；一维热传导方程 ；边界条件 中图分类号 ：Ｏ５５１ ．３ 文献标识码 ：Ａ 文章编号 ：１６７１－６８７６（２００４）０３－０２１０－０５
数理方程中混合问题的解法通常是比较复杂的 ，有时很难有解析解 ，即便解出 ，其解的形式也通常 是一个无穷级数形式［１ ４］ ．本文以一维热传导混合问题为例 ，探讨利用差分法对数理方程进行数值解的 方法 ．
（３３）
ｕ（ｘ ，０） ＝ １００ ℃ ｘ ／ｌ （０ ≤ ｔ ≤ ｌ）
（３４）
图 ８ 中 ｔ ＝ ０ 是初始状态 ，为一直线 ，右端为 ｘ ＝ ｌ ，每一条曲线左边的开始端平坦 ，这正是由于在 ｘ ＝ ｌ
处为第二类齐次边界条件的原因 ．而右端 ｘ ＝ ０ 的温度始终为 ０ ，这是由于该点的边界条件为第一类齐
图 ４ ２１ 类边界条件下的运行图 图 ５ １２ 类边界条件下的运行图
２ ．２ １２ 类边界条件的处理办法
设 ｘ ＝ ｌ 端满足第二类边界条件 ，即
ｕｘ （ ｘ ，ｔ） ｜ ｘ ＝ ｌ ＝ μ（ ｔ）
（２３）
那么也不能直接利用（７） 式 ，必须首先利用（２３） 式及初始条件（３） 逐步求出边界 ｘ ＝ ｌ 及其他各点处各
２１２
淮阴师范学院学报 （自然科学版 ）
第３卷
从结果可以知道热传导杆两端的温度始终保持为 ０ ，杆的中点温度总是高于其它点温度 ，各点温度 随着时间变化逐渐降低 ．
２ 含第二类边界条件的一维热传导混合问题的数值解法
２ ．１ ２１ 类边界条件的处理办法
设 ｘ ＝ ０ 端满足第二类边界条件 ，即
第 ３ 卷第 ３ 期 ２００４ 年 ８ 月
淮阴师范学院学报 （自然科学版 ）
ＪＯＵＲＮＡＬ ＯＦ ＨＵＡＩＹＩＮ ＴＥＡＣＨＥＲＳ ＣＯＬＬＥＧＥ （ＮＡＴＵＲＡＬ ＳＣＩＥＮＣＥ ＥＤＩＴＩＯＮ）
Ｖｏｌ畅３ Ｎｏ畅３ Ｍａｙ ．２００４
一维热传导方程的数值解
徐建良 ，汤炳书
（连云港高等师范专科学校 物理系 ，江苏 连云港 ２２２００６）
（２）
ｕ（ｘ ，０） ＝ φ（ｘ） （０ ≤ ｘ ≤ ｌ）
（３）
这是一个 １１ 类非齐次边界条件的一维热传导问题 ，通常这一类混合问题是很难解的 ，即便解出 ，其解也
通常是一个无穷级数的形式 ，对该解的物理意义不能直接讨论 ，不能给出直观的图象 ．
１ ．１ 计算方法
为求解方程（１） ，首先定义函数 ｕ（ｘ ，ｔ）的时间与空间的网格 ，将 ｘ 坐标分成 Ｎ 等份 ，将 ｔ 坐标分成 Ｍ
（２８）
这样就可先由（２８） 式算出 ｉ ＝ Ｎ 上的各点温度值 ，再由（２７） 式算出边界 ｘ ＝ ｌ ，即 ｉ ＝ Ｎ ＋ １ 上各点的
温度值 ．（２７） 和（２８） 式运算过程如图 ５ ．
图 ６ 例 ３ 中相同 ｔ 不同 ｘ 的温度变化曲线 图 ７ 例 ３ 中相同 ｘ 不同 ｔ 的温度变化曲线
＝
μ（ｊ）
（１８）
所以
ｕ０ ，ｊ ＝ ｕ２ ，ｊ － ２ Δｘμ（ ｊ ）
（１９）
令（７） 式中 ｉ ＝ １ 得
ｕ１ ，ｊ＋ １ ＝ α（ ｕ０ ，ｊ ＋ ｕ２ ，ｊ ） ＋ （１ － ２ α） ｕ１ ，ｊ ＋ Δｔｆ （１ ，ｊ ）
（２０）
由（１９） 、（２０） 式解得
ｕ１ ，ｊ＋ １ ＝ ２ α（ ｕ２ ，ｊ － Δｘμ（ ｊ ）） ＋ （１ － ２ α） ｕ１ ，ｊ ＋ Δｔｆ（１ ，ｊ）
Ｎｕｍｅｒｉｃａｌ Ｃａｌｃｕｌａｔｉｏｎ ｏｆ Ｏｎｅ－ｄｉｍｅｎｓｉｏｎａｌ Ｈｅａｔ Ｃｏｎｄｕｃｔｉｏｎ Ｅｑｕａｔｉｏｎ
ＸＵ Ｊｉａｎ－ｌｉａｎｇ ，ＴＡＮＧ Ｂｉｎｇ－ｓｈｕ
（Ｄｅｐａｒｔｍｅｎｔ ｏｆ Ｐｈｙｓｉｃｓ ，Ｌｉａｎｙｕｎｇａｎｇ Ｔｅａｃｔｈｅｒｓ Ｃｏｌｌｅｇｅ ，Ｌｉａｎｙｕｎｇａｎｇ ２２２００６ ，Ｃｈｉｎａ）
图 ８ 例 ４ 中相同 ｔ 不同 ｘ 的温度变化曲线 图 ９ 例 ４ 中相同 ｘ 不同 ｔ 的温度变化曲线
图 ６ 中 ｔ ＝ ０ 是初始状态 ，为一直线 ，左端为 ｘ ＝ ０ ，每一条曲线左边的开始端平坦 ，这正是由于在 ｘ ＝ ０ 处为第二类齐次边界条件的原因 ．而右端 ｘ ＝ ｌ 的温度始终为０ ，这是由于该点的边界条件为第一类
ｕ（ｘ ，ｔ） ＝
Ａ π２
ｌ２ ａ２
（１
－
ｅ － π２ ａ２ ｔ／ ｌ２ ）ｓｉｎ
πｘ ｌ
由数值计算得出的结果如图 ２ 所示 ．
例 ２ 设定解问题为
ｕｔ － ａ２ ｕｘｘ ＝ ０ （０ ＜ ｘ ＜ ｌ） （０ ＜ ｔ ＜ ∞ ）
ｕ（０ ，ｔ） ＝ ０ ， ｕ（ ｌ ，ｔ） ＝ ０ （０ ≤ ｔ ≤ ∞ ）
例 ３ 设一维热传导问题为
ｕｔ
＝
ａ２ ｕｘｘ
＋
ｇｓｉｎ
πｘ ｌ
（０
＜
ｘ＜
ｌ） （０ ＜
ｔ＜
Hale Waihona Puke ∞）（２９）ｕｘ （０ ，ｔ） ＝ ０ ， ｕ（ ｌ ，ｔ） ＝ ０ （０ ≤ ｔ ≤ ∞ ）
（３０）
ｕ（ｘ ，０） ＝ １００ ℃ （ ｌ － ｘ）／ｌ （０ ≤ ｘ ≤ ｌ）
（３１）
这是一个 ｘ ＝ ｌ 端为第二类齐次边界条件的并且具有热源的热传导混合问题 ，即 μ（ ｔ） ＝ ０ ．在编程时 ， 温度采用国际温标 ．设 ｌ ＝ １ｍ ，时间范围为 ０ ～ １ ｓ ，ｇ ＝ ２７３０ ｋ／ｓ ，将 ｌ 分割为 １００ 份 ，时间分为 ５００ 份 ． 时空网格有 １０１ × ５０１ 个格点 ．例 ３ 的数值计算结果如图 ６ 、图 ７ 所示 ．
ｕｘ （ ｘ ，ｔ） ｜ ｘ ＝ ０ ＝ μ（ ｔ）
（１７）
则此时就不能直接利用（７） 式求解 ，必须首先利用（１７） 式及初始条件（３） 逐步求出边界 ｘ ＝ ０ 及其
他各点处各时刻的温度值 ．因 ｘ ＝ ０ 处对应 ｉ ＝ １ ，则由（１７） 式得
抄 ｕ１ ，ｊ 抄ｘ
＝
－ ｕ０ ，ｊ ＋ ｕ２ ，ｊ ２ Δｘ