3差分格式

§3. 热传导方程

上一节曾指出，由于定解问题中的每一个偏导数都有多种差分近似，所以一个定解问题可以有多个不同的差分格式。下面以热传导方程为例，对此展开讨论。为简单起见，先不给出定解条件。

考虑热传导方程

2

抖uu

= 2 ¶t¶x

¶u¶2u

方程中出现了一阶时间导数 和二间空间导数 2 。

¶t¶x

对于一阶时间导数

¶u

，常用的差分近似就有三种，记 ¶t

n

+1nun-ujj

向前差分近似

¶u

禗t

Ü

j

t

n-1

un-ujj

向后差分近似

¶u

禗t

n

Ü

j

t

中心差分近似

¶u禗t

n

Ü

j

+1n-1

un-ujj

2t

这里，我们用“Ü”代表上一节推导差分近似的过程。

前面已经看到，二阶导数通常用中心差分近似。对这里的二阶空间导数的中心差分近似为

n

nn

un-2u+uj+1jj-1

¶u

禗x2

2

Ü

j

x

2

但是也可以考虑其他可能的方案。由泰勒展开，有

n+1

n

n

n

抖u

抖x2

2

=

j

u抖uu

+Dt+L=+O(Dt) 222

xj抖txj xju抖uu

-Dt+L=+O(Dt) 222

xj抖txj xj

n+1

2

n

3

n

2

n

232

抖u

抖x2

2

n-1

=

j

所以，如果用

¶u

¶x2

2

或

j

¶u¶x2

2

n-1

代替

j

¶u

，虽然会引入新的误差，2

¶xj

2

n

但这种误差与差分近似已有的误差为同一量级的，因而还是可以接受的。这样一来，二阶空间导数的差分近似又有了两种新的方案

n

n+1

+1n+1n+1un-2u+uj+1jj-1

抖u

抖x2

2

苘

j

u

x2

2

j

Dx

2

抖u抖x2

2

n

苘

j

u

x2

2

n-1

-1n-1n-1un-2u+uj+1jj-1

j

Dx

2

将前面讨论过的一阶时间导数的三种差分近似与这里给出的二阶空间导数的三种差分近似方案搭配起来，可以为热传导方程构造出九种不同的差分格式，如下表所示。

3

（注）二阶空间导数差分近似的三个方案可以叙述为：在用中心差分近似二阶空间导数时，时间自变量 t 的取值固定在 t=tn 、t=tn+1 或

t=tn-1 三个不同的时刻。

在上面的表中，除了给出差分格式，还给出了这个格式所用到的网格点的图示，称为格式的模板。另外，表中打 ´ 表示这个格式在网格点的使用上明显不合理，下面将不予考虑。在余下的格式中，我们挑出三个作为基本格式（打✓者，后面将会解释这三个名称的含义），它们是

+1n

un-ujj

nn

un-2u+uj+1jj-1

两层显式格式

Dt

+1n

un-ujj

=α

Dx

2

两层隐式格式

Dt

=α

+1n+1n+1

un-2u+uj+1jj-1

Dx

2

三层显式格式

+1n-1

un-ujj

2Dt

=α

nn

un-2u+uj+1jj-1

Dx

2

对表中标有 ❶、❷、❸、❹ 的这四个格式，有以下的分析。 对于 ❶ 和 ❷ ，令 n¢=n-1 ，则 n=n¢+1 ，于是它们可写成

+1n

un-ujj

ⅱ

nn

un-2u+uj+1jj-1

ⅱ

Dt

+1n

un-ujj

ⅱ

=α

Dx

ⅱ

2

Dt

=α

+1n+1n+1

un-2u+uj+1jj-1

Dx

2

显然它们与前两个基本格式是类似的。

而对于格式 ❸ 和 ❹ ，仍令 n¢=n-1 并取 Dt¢=2Dt ，则由于 tn+1=tn-1+2Dt=tn-1+Dt¢ ，所以如果用 Dt¢=2Dt 作为时间步长，tn+1 就是 tn-1 的下一个时刻，可记作 nⅱ+1 ，于是这两个格式也可写成

+1nnnn

un-uu-2u+ujjjj-1

=αj+1

2

¢DtDx+1nn+1n+1n+1

un-uu-2u+ujjjj-1

=αj+1

2

Dt¢Dxⅱ

ⅱ

ⅱ

ⅱ

ⅱ

ⅱ

ⅱ

可见它们也与前两个基本格式类似。

综合以上的分析，这四个格式与基本格式重复，也不必再考虑。 以下将集中分析上面选出的三个基本格式。前两个基本格式，只用到

n 和 n+1 两个时间层的网格点，故称为两层格式。而第三个格式用到

了 n-1、 n、 n+1 三个时间层的网格点，所以称为三层格式。

下面为热传导方程补上定解条件，考虑热传导方程定解问题

2ìï抖uuï=α , a＃xb , t 0ï2ï¶t¶xïïïïïua,t=gt , ub,t=gt , t>0

í()()b()a()ïïïïïïu(x,0)=ϕ(x) , a

对这个定解问题，第一个基本格式为

n+1nnnnìïu-uu-2u+uïjjj+1jj-1ï=α2ïDtïDxïïï j=1,2,L,M-1 , n=1,2,L,N-1ïï íïïnnnnïu=gt , u=gt , n=1,2,L,N-10aMbïïïïïuj0=ϕ(xj) , j=1,2,L,M-1ïïî

()()

将它改写成便于计算的形式

0ìïu=ϕ(xj) , j=1,2,L,M-1ïjïïïïn+1nnnnïu=u+ru-2u+uïjj+1jj-1ïj

í

ïï j=1,2,L,M-1 , n=1,2,L,N-1ïïïïïnnnnïu=gt , u=gt , n=1,2,L,N-10aMbïïî

()

()()

式中 r=α

Dt

。这个格式的求解步骤如下： 2

Dx

t=0时刻（n=0） ：uj0=ϕxj

0000

t=t1时刻（n=1） ：u1=u+ru-2u+ujjj+1jj-1

()

()

，(j=1,2,L,M-1)

11

u0=gat1 ，uM=gbt1

()()

111

t=t2时刻（n=2） ：uj2=u1+ru-2u+ujj+1jj-1

()

，(j=1,2,L,M-1)

22

u0=gat2 ，uM=gbt2

()()

一般地，如果已经求出 t=tn 时刻的近似解，则

+1nnnn

t=tn+1时刻：un=u+ru-2u+ujjj+1jj-1

()

=1,2,L,M-1)

，(j

n+1n+1

u0=gatn+1 ，uM=gbtn+1

()()

对 n=0,1,2,L,N-1 反复进行这一步骤，直到求出 t=tN=T 时刻的近似解。

从这个计算过程可以看出，只要用初始条件给出第0层网格点上的解，就可以用这个格式一步一步地计算出各个时间层的近似解。这样的格式称为显式格式。

再来考虑第二个基本格式，它可以改写成

+1n+1n+1n

-run+1+2ru-ru=u()j-1jj+1j

+1n+1n+1

显然，第 n+1 层的三个近似解 un、、 是耦合在一起的，无uuj-1jj+1

法单独计算出来。所以，需要在第 n+1 层的每一个网格点上列出差分格式，联立求解，即

n+1n+1nn+1ìï1+2ru- ru =u+ru()ï1210

ïïïïn+1n+1n+1nï -ru+1+2ru-ru =u()1232ïïïïï Oïïïïïn+1n+1n+1nï í -ruj-1+(1+2r)uj-ruj+1 =ujïïïïï Oïïïïn+1n+1n+1nï -ru+1+2ru- ru=uï()M-3M-2M-1M-2ïïïïn+1n+1nn+1ï -ru+1+2ru=u+ru()ïM-2M-1M-1Mïî

n+1n+1方程组中的 u0=gatn+1 ，uM=gbtn+1 都是已知的（由边界条件

()()

事先确定的），所以可以移到方程右边。

如果将第 n+1 层的近似解和上述方程组的右端项写成向量

n+1nn+1骣骣uu+ru÷÷çç110÷÷çç÷÷çç÷÷çç÷÷çç÷n+1÷n

÷÷ççuu÷÷çç22÷÷çç÷÷çç÷÷çç÷÷çç÷÷÷÷ççMM÷÷çç÷÷çç÷÷çç÷÷çç÷÷n+1nn+1çç÷÷u÷u ， =ç÷b==çjj÷÷çç÷÷çç÷÷çç÷÷çç÷÷çç÷÷M÷M÷çç÷÷çç÷÷çç÷÷çç÷÷çç÷÷n+1nç÷ç÷uu÷÷ççM-2M-2÷÷çç÷÷çç÷÷çç÷÷çç÷n+1÷nn+1琪琪uu+ruçç桫M-1÷桫M-1M÷

un+1

nn+1骣u+rgtç1aççççnçuç2

çççççMçççççnçujçççççMçççççnçuMç-2ççççnn+1çu+rgtçb桫M-1

()

()

÷÷÷÷÷÷÷÷÷÷÷÷÷÷÷÷÷÷÷÷÷ ÷÷÷÷÷÷÷÷÷÷÷÷÷÷÷÷÷÷÷÷÷

再将方程组的系数矩阵记作

骣1+2r-rçççç-r1+2r-rçççççOOçççç

-rA=çççççççççççççççç桫

则方程组可写成

O1+2rO

-rO-r

÷÷÷÷÷÷÷÷÷÷÷÷÷÷÷÷ ÷÷÷÷÷÷O÷÷÷÷÷1+2r-r÷÷÷÷

÷-r1+2r÷

Aun+1=bn+1

这个格式的求解步骤如下：

t=0 时刻（n=0） ：由初始条件 uj0= xj 计算 u0 。

()

t=t1 时刻（n=1） ：利用 u0 和边界条件计算右端项 b1 ，

解方程组 Au1=b1 ，求出 u1 。

t=t2 时刻（n=2） ：利用 u1 和边界条件计算右端项 b2 ，

解方程组 Au2=b2 ，求出 u2 。

一般地，如果已经求出 t=tn 时刻的近似解 un ，则

t=tn+1 时刻 ：利用 un 和边界条件计算右端项 bn+1 ，

解方程组 Aun+1=bn+1 ，求出 un+1 。

对 n=1,2,L,N-1 反复求解方程组，直到求出 t=tN=T 时刻的近似解。

从上述计算过程可以看出，这个格式一个时间步的计算都需要求解一个方程组，才能得到这个时间层的近似解。这样的格式称为隐式格式。

上面对显式格式和隐式格式的分析表明，从计算量和编程工作量的角度来说，显式格式由于不需要求解方程组，计算量小，便于编程，所以优于隐式格式。

最后，简单分析一下第三个基本格式，即三层格式，它可以改写成

+1n-1nnn

un=u+2ru-2u+ujjj+1jj-1

()

可以看出，这是一个显式格式，无需解方程组，这是它的优点。另一方面，作为三层格式，计算第 n+1 层的近似解需要用到前两个时间层（第 n 层和第 n-1 层）的近似解，那么在开始计算时就只能从第2层开始（用第1层和第0层的近似解来计算）。但是与上一节的波动方程不同，这里的热传导方程定解问题只有一个初始条件 ux,0=ϕx ，它只能给出第0层的解。第1层的近似解必须用其他方法确定（有时这很困难）。这是三层格式的缺点（对于热传导方程而言，所有三层格式都有这个缺点）。

（注）热传导方程的差分格式，无论是两层格式还是三层格式，无论是显式格式还是隐式格式，都会出现一个无量纲的量，就是比值

()()

r=α

Dt

，下面给出的格式也是如此。由此可见，r 是一个重要的量，2

Dx

称为（热传导方程差分格式的）网格比。

除了本节已经讨论过的三个基本格式之外，热传导方程还有没有其它形式的差分格式？答案显然是肯定的。这里再给出一个比较常用的格式，作为这一节的结束。

在三个基本格式中有两个是两层格式，它们的区别在于：用中心差分近似二阶空间导数时，时间自变量 t 的取值固定在 t=tn 时刻还是t=tn+1 时刻。这两种选择，分别得到了显式格式和隐式格式。如果不想做这个选择，可以考虑用这两者的加权平均，于是得到差分格式

+1nun-ujjnnn轾un+1-2un+1+un+1u-2u+uj+1jj-1j+1jj-1 =α犏θ+1-θ()22DxDx犏臌Dt

式中的 θ 是权（0＃θ。这个格式称为六点格式，其模板是

1）

六点格式可改写成

+1n+1n+1 -θrun+1+2θru-θru()j-1jj+1

nn轾=(1-θ)run+1-21-θru+1-θru()()j-1jj+1犏臌

所以，只有 θ=0 时，六点格式才是显式格式，就是前面的两层显式格式。而当 θ¹0 时，六点格式总是隐式格式。特别是当 θ=1 时，就是前面给出的两层隐式格式。从现在起，为了与 θ 取其它非零值的两层隐式格式相区别，将它改称为两层完全隐式格式。

除了这两个已知的格式之外，还有一个常用的格式，就是取 θ=

称为两层平均隐式格式或Crank-Nicholson（C-N）格式，形式为

+1nun-ujjn+1n+1n+1nnn轾u-2u+uu-2u+u1犏j+1jj-1j+1jj-1 =α+222犏DxDx臌1 ，2Dt

对于三层格式也有类似的结果。例如，三层平均隐式格式（九点格式） +1n-1un-ujjn+1n+1n+1nnnn-1n-1n-1轾u-2u+uu-2u+uu-2u+u1j+1jj-1j+1jj-1j+1jj-1 =α

犏++2223犏DxDxDx臌2Dt