热传导方程几种差分格式的MATLAB数值解法比较

来源：网络收集时间：2026-03-26

导读：第25卷第2期沈阳化工大学学报 2011．06文章编号： JOURNALOFSHENYANGUNIVERSITYOFCHEMICALTECHNOLOGY Vol．25No．2 Jun．2011 2095－2198（2011）02－0179－04 热传导方程几种差分格式的MATLAB数值解法比较冯立伟（沈阳化工大学数理系，辽宁沈阳110142）

第25卷第2期沈阳化工大学学报

2011．06文章编号：

JOURNALOFSHENYANGUNIVERSITYOFCHEMICALTECHNOLOGY

Vol．25No．2

Jun．2011

2095－2198（2011）02－0179－04

热传导方程几种差分格式的MATLAB数值解法比较

冯立伟

（沈阳化工大学数理系，辽宁沈阳110142）

摘

要：

对求解热传导方程的几种差分格式进行分析，并讨论使用MATLAB编程求解偏微分方

程的方法．编制几种差分格式的MATLAB程序，使用算例进行数值实验，在不同网格比情况下，比较几种算法的优劣．关键词：

热传导方程；

O241.82

MATLAB；

Crank-Nicolson离散

中图分类号：

文献标识码：

许多工程问题需要研究热量在物体内部的

传导情况或某种物质在液体中的扩散情况，因此研究热传导问题特别是非稳态热传导问题十分

［1－2］

．重要．目前热传导方程已有多种求解格式

MATLAB是目前最流行、应用最广泛的科学和工程计算软件．MATLAB基于矩阵运算，具有强大的数值运算能力和图形可视化能力，是方便实

［3］

用、功能强大的数学软件．用MATLAB求解常

［4－6］

．王飞等介绍了如微分方程已有大量的研究

何使用MATLAB实现有限差分法求解微分方程．高理平等给出了对两点边值问题有限元

［8］

方法的程序．李灿等对热传导问题的MAT-LAB数值计算进行了讨论［9］．本文讨论求解一维热传导方程几种不同差分格式的MATLAB编程方法，并使用算例进行检验和对结果进行分析．

uu

=a2+f（x，t）tx［7］

为空间步长与时间步长，用2族平行直线xj=jh，tk=kτ将矩形区域［0，T］×［0，L］分割成矩形网格．

显式格式：

Uk－Ukjj

=τ

Ukj+1－2Uj+Uj－1

+fkaJ．2

（2）

隐式格式：+1Uk－Ukjj

=τ

+1k+1+1Uk+Ukj+1－2Ujj－1a+fkJ．2

（3）

Crank-Nicolson格式：

+1+1k+1+1Uk－UkUk+Ukjjj+1－2Ujj－1

=a+τh2

(

（1）

u（x，t）表示在t时刻物体内部坐标为x处的温

a是热传导系数，f（x，t）为热源．度，

Ukkj+1－2Uj+Uj－1

+fj．2

（4）

DuFortFrankel格式：

+1－1Uk－Ukjj

2τ

k+1－1

Uk－Uk+Ukj+1－Ujjj－1a+fkj．2

1.1

热传导方程差分格式

差分格式的建立

t平面进行网格剖分．分别取h，首先对x-τ

（5）

由Taylor公式容易得出：它们都与一维热传导

O（τ+方程相容，其截断误差分别为O（τ+h），

收稿日期：作者简介：

2010－08－30

冯立伟（1980－），男，辽宁沈阳人，助教，硕士，主要从事偏微分方程数值解的研究．

h2），O（τ2+h2）和O（τ2+h2）［1］．1.2

初、边值条件的处理

可算出u在第一层各个节点处的近似值uj．重

可以逐层计算出所有的uj．复使用此式，

隐式格式：将（7）式与离散化的初边值条件联立，得差分方程组：

k+1k+1k+1kk

－rUj－1+（1+2r）Uj－rUj+1=Uj+τfj 2，…，N－1，j=0，1，…，M－1）（k=1，

（11）

U=（k=0，1，…，N）jj k

U0=g1k，Uk（j=0，1，…，M）L=g2k将上述方程组改写成矩阵形式

对定解条件进行离散化．由初始条件及第一

类边界条件，可直接得到：

U00）=j，（j=0，1，…，M）；j=u（xj，

U0=u（0，tk）=g1k，

Uktk）=g2k，（k=1，…，N）．L=u（L，1.3

稳定性分析令r=

aτ

为网格比，显式格式公式（2）变为：h2

+1kUk=rUkj+1+（1－2r）Uj+J

krUk（6）j－1+τfj

使用Fourier方法可知，当r≤1/2时显式格式稳定．

隐式格式变为：

+1k+1

－rUk－j－1+（1+2r）Uj

+1kkrUkj+1=Uj+τfj．

（7）

由Fourier方法可得隐式格式恒稳定．Crank-Nicolson格式变为：

－

rk+1+1Uj－1+（1+r）Uk－j2

k+1k

U1 U1+rg1j+1

k Uk+1 U22

︙ = ︙ k+1 UM－2 UkM－2 k+1 k UM－1 UM－1+rg2j+1

（12）

此方程组是三对角方程组，且系数矩阵严格对角

故解存在唯一．通过在每一时间层上求解占优，

一个这样的线性方程组得到在各个时刻各个网格点上的函数值．

（8）

Crank-Nicolson格式：将（8）式与离散化的

得差分方程组：初边值条件联立并整理，

r+1+1k+1

－rUk－Ukj－1+（1+r）Ujj+1=22

rkkk

2Uj+1+（1－r）Uj+2Uj－1+τfj

2，…，N－1，j=0，1，…，M－1）（k=1，

U0=（j=0，1，…，M）jj k U0=g1k，Uk（k=1，2，…，N）L=g2k

rk+1rk

U=U+2j+12j+1（1－r）Ukj+

Uj－1+τfkj．2

对于r＞0恒有增长因子|G|≤1，恒稳定．

DuFortFrankel格式变为：

+1－1

（1+2r）Uk=（1－2r）Uk+jj

2rU

j－1

+2rU

kj+1

+2τf．

（9）

（13）

由Fourier方法可得DuFortFrankel格式恒稳定．

2差分格式的求解

显式格式：将（6）式与离散化的初边值条件结合，得到求解此问题的差分方程组：

k+1kkkk Uj=rUj+1+（1－2r）Uj+rUj－1+τfj 2，…，N－1，j=0，1，…，M－1）（k=1， 0

Uj=j（j=0，1，…，M） k

U0=g1k，Uk（k=1，2，…，N）L=g2k

与隐式格式类似，用六点格式由第k时间层的值

k+1

Uk时，需求解三j计算第k+1时间层的值Uj对角方程组

：

（10）

由于初始时间层上的u值为已知，由（10）式即

end

DuFortFrankel格式：将（7）式与初始条件及第一类边界条件式联立并整理将格式改写为：

k+1k－1k

1+2r）U=（1－2r）U+2rU（j－1+jj

k 2rUkj+1+2τfj

（15）2，…，N－1，j=0，1，…，M－1）（k=1，

U0=（j=0，1，…，M－1）jj k

U0=g1k，Uk2，…，N－1）L=g2k（k=1，

k+1 U 2 ︙ k+1 UM－2 k+1 UM－1

k+11

fk+rUk+Uk+1）

12 k

+τ ︙

fk M－2

rk+1

fk（UkM－1+M+UM） 2

（14）

由于此格式是3层格式，需要事先知道前面2个时间层上的解，第1层上的解通过对初始条件离散可得，第2层上的解使用前面的任意1种2层差分格式得到，再用此格式求解其余时间层上的解．程序核心代码如下：

%GeneratesecondrowusingimplicitschemeA=（1+2*r）*eye（M－2）；temp=－r*linspace（1，1，M－3）；A=A+diag（temp，1）+diag（temp，－1）；b=U（1，2：M－1）'；

b（1）=b（1）+r*U（1，1）；b（end）=b（end）+r*U（1，M）；

U（2，2：M－1）=A/b；%solveequationAU=b%GenerateremainingrowsofUforkt=3：N

forjx=2：M－1U（kt，jx）=（2*r*U（kt－1，jx－1）+2*r*U（kt－1， …… 此处隐藏：2761字，全部文档内容请下载后查看。喜欢就下载吧 ……

热传导方程几种差分格式的MATLAB数值解法比较.doc 将本文的Word文档下载到电脑，方便复制、编辑、收藏和打印

下载这篇word文档

本文链接：https://www.jiaowen.net/wenku/97825.html（转载请注明文章来源）

上一篇：股权转让协议(适用于有限责任公司)2014年最新完整版
下一篇：热电公司年度安全生产目标责任落实情况自查汇报