随着科技的进步和计算机技术的快速发展, 在自然科学和工程计算领域诞生了许多实际问题, 人们发现这些问题可以通过微分方程恰当的表示出来。双曲型方程作为一类非常重要的偏微分方程, 主要描述振动、波动等现象的运动过程。对双曲型方程的研究一直是非常重要的课题, 其应用范围包括流体力学、空气动力学、光学、电磁学, 以及广义相对论等^[1]。

现如今, 对于微分方程的解, 可以归结为两大类, 一类是精确解(解析解), 另一类是数值解。由于精确解只有在一些特殊的情况下才可以求出, 因此研究微分方程的数值解就变得非常重要。对于偏微分方程数值解的研究已经有很长的历史, 目前应用比较广泛的方法主要有有限差分法、有限元法、谱方法等^[2-4]。

本文主要考虑有限差分法。有限差分法是求解微分方程的一种成熟、高效的方法, 其主要思想为以差商近似导数, 具体步骤可分为: 对方程的求解区域进行网格剖分; 对微分算子进行离散; 建立以网格节点上的值为未知数的代数方程组等^[2]。有限差分法在实际生活中已有十分广阔的应用, 基于有限差分法进行数值模拟, 可以有效解决合金冶炼领域中铸件冷却不均匀的问题^[5]; 有限差分法在油藏模拟实验、地下水勘探等领域也发挥了不可替代的作用^[6]; 依托于计算机技术的快速发展, 有限差分法也被应用到新材料的研发当中^[7-8]。

将有限差分法应用于双曲型方程的求解中, 具有精度高、时效短等特点^[9-11]。但由于有限差分法对微分算子进行离散后的格式是多样的, 为了衡量差分格式的好坏, 还需要讨论差分格式的稳定性^[4]。差分格式的稳定性问题是决定差分格式在计算时能否求得理想数值解的关键问题, 只有稳定的差分格式才可能是有用的格式。本文基于分离变量法, 对二维二阶双曲型方程隐式差分格式的稳定性做了详细的讨论。

1 预备知识 1.1 差分格式的建立 1.1.1 网格剖分

为了构造微分方程的有限差分逼近, 首先要对求解区域进行网格覆盖, 以二维区域为例, 取沿x轴方向(空间方向)的步长为h, t轴方向(时间方向)步长为k, 分别作两族平行于x轴和t轴的直线

$ \begin{array}{l} x = {x_i} = ih, \quad i = 0, \pm 1, \pm 2, \cdots , \\ t = {t_n} = nk, \quad n = 0, 1, 2, \cdots 。\end{array} $

以上两簇直线的交点(x_i, t_n)称为网格的节点, 这里取的是与坐标轴平行的等距直线。

1.1.2 有限差分近似

利用Taylor展开式, 将微分方程在节点处进行离散化, 即以差商近似偏导数。以空间部分x为例, 由Taylor展开式得

$ u\left( {{x_{i + 1}}, {t_n}} \right) = u\left( {{x_i}, {t_n}} \right) + h\frac{{\partial u}}{{\partial x}}\left( {{x_i}, {t_n}} \right) + \frac{{{h^2}}}{2}\frac{{{\partial ^2}u}}{{\partial {x^2}}}\left( {{x_i}, {t_n}} \right) + \cdots , $

(1)

$ u\left( {{x_{i - 1}}, {t_n}} \right) = u\left( {{x_i}, {t_n}} \right) - h\frac{{\partial u}}{{\partial x}}\left( {{x_i}, {t_n}} \right) + \frac{{{h^2}}}{2}\frac{{{\partial ^2}u}}{{\partial {x^2}}}\left( {{x_i}, {t_n}} \right) + \cdots 。$

(2)

由(1)式, 忽略截断误差, 便得到了空间方向的一阶近似偏导

$ \frac{{\partial u}}{{\partial x}}\left( {{x_i}, {t_n}} \right) \approx \frac{{u\left( {{x_{i + 1}}, {t_n}} \right) - u\left( {{x_i}, {t_n}} \right)}}{h}, $

(3)

其截断误差为O(h)。

(1) 式和(2)式相加, 忽略截断误差, 便得到了空间方向的二阶近似偏导

$ \frac{{{\partial ^2}u}}{{\partial {x^2}}}\left( {{x_i}, {t_n}} \right) \approx \frac{{u\left( {{x_{i + 1}}, {t_n}} \right) - 2u\left( {{x_i}, {t_n}} \right) + u\left( {{x_{i - 1}}, {t_n}} \right)}}{{{h^2}}}, $

(4)

其截断误差为O(h²)。与(3)式相比, (4)式具有更高的计算精度, 但提高精度的同时也带来了更高的计算复杂度。在实际问题中选用哪种差分格式, 需要具体情况具体分析。

在计算微分方程数值解的过程中, 不仅需要对求解区域进行离散, 还需要对边界条件和初始条件在网格点处进行离散。本文仅对差分格式的稳定性进行分析, 不涉及具体问题的求解, 这一部分内容不再展开介绍, 详细介绍可见参考文献^[2]。

1.2 分离变量法(Fourier方法)

在这一小节中, 我们将介绍分析差分格式稳定性的一种方法——分离变量法。分离变量法是分析差分方程稳定性常用的方法之一^[4], 首先给定初始时刻的误差, 假设在差分方程求解的过程中并未产生其它误差, 进而分析初始误差在沿时间轴传播的过程中是否带来明显的波动。

假设微分方程在节点(x_j, t_n)处的精确解为u(x_j, t_n), 利用差商近似导数的方法, 可以获得方程的近似解U(x_j, t_n), 此即差分方程的解。U(x_j, t_n)可表示为Fourier级数的形式U_jⁿ=ξ_ηⁿe^iβjh, 其中U_jⁿ为近似解U(x_j, t_n)的简写, j为格式所在的空间层数, n为格式所在的时间层数; i为虚数单位; β=2πη, η∈ $\mathbb{N} $; h为沿x轴方向的步长。

分离变量法的具体步骤如下:

(1) 将U_jⁿ=ξ_ηⁿe^iβjh代入具体的差分格式;

(2) 等式两边消去公因子e^iβjh;

(3) 若格式关于时间是两层的, 可以得到差分格式的增长因子或增长矩阵; 若格式关于时间是三层的, 可以得到差分格式的特征方程。

2 二维二阶线性双曲型方程隐式差分格式稳定性分析

考虑如下二维二阶线性双曲型方程

$ \frac{{{\partial ^2}u}}{{\partial {t^2}}} = \frac{{{\partial ^2}u}}{{\partial {x^2}}} + \frac{{{\partial ^2}u}}{{\partial {y^2}}}。$

(5)

取沿x轴和y轴方向的步长均为h, 沿t轴方向的步长为k。考虑如下的交替方向隐式格式

$ \frac{1}{{{k^2}}}\delta _t^2U_{i, j}^n = \frac{1}{{{h^2}}}\left( {\delta _x^2 + \delta _y^2} \right)\left( {\theta U_{i, j}^{n + 1} + (1 - 2\theta )U_{i, j}^n + \theta U_{i, j}^{n - 1}} \right), $

(6)

其中0≤θ≤1。当θ=0时, 格式(6)便转化为显式差分格式, 这里不再对其稳定性做详细的讨论^[2]。显然(6)式是一个三层格式, 其中δ_t², δ_x², δ_y²分别为关于t, x和y的二阶中心差分算子, 具体表达式如下

$ \begin{array}{l} \delta _t^2U_{i, j}^n = U_{i, j}^{n + 1} - 2U_{i, j}^n + U_{i, j}^{n - 1}, \\ \delta _x^2U_{i, j}^n = U_{i + 1, j}^n - 2U_{i, j}^n + U_{i - 1, j}^n, \\ \delta _y^2U_{i, j}^n = U_{i, j + 1}^n - 2U_{i, j}^n + U_{i, j - 10}^n。\end{array} $

引入网格比r的定义, 记

$ r = \frac{k}{h}, $

(7)

则(6)式可变形为

$ \begin{array}{l} U_{i, j}^{n + 1} - 2U_{i, j}^n + U_{i, j}^{n - 1}\\ = {r^2}\left[ {\theta \left( {U_{i + 1, j}^{n + 1} - 2U_{i, j}^{n + 1} + U_{i - 1, j}^{n + 1}} \right) + (1 - 2\theta )\left( {U_{i + 1, j}^n - 2U_{i, j}^n + U_{i - 1, j}^n} \right) + \theta \left( {U_{i + 1, j}^{n - 1} - 2U_{i, j}^{n - 1} + U_{i - 1, j}^{n - 1}} \right)} \right.\\ \left. { + \theta \left( {U_{i, j + 1}^{n + 1} - 2U_{i, j}^{n + 1} + U_{i, j - 1}^{n + 1}} \right) + (1 - 2\theta )\left( {U_{i, j + 1}^n - 2U_{i, j}^n + U_{i, j - 1}^n} \right) + \theta \left( {U_{i, j + 1}^{n - 1} - 2U_{i, j}^{n - 1} + U_{i, j - 1}^{n - 1}} \right)} \right]。\end{array} $

(8)

借助于1.2节介绍的分离变量法, 本文将利用两种不同的方法求差分格式(8)的特征方程; 然后借助于分类讨论的思想, 探究θ在不同情形下对应差分格式的稳定性。

2.1 计算差分格式的特征方程

本节将利用两种方法计算差分格式的特征方程。第一种方法直接使用分离变量法, 通过化简整理得到待求的特征方程; 第二种方法首先将给定的三层隐式格式改写为两层差分方程组, 再使用分离变量法将其转化为矩阵形式, 写出其增长矩阵, 增长矩阵的特征多项式即为该差分格式的特征方程。

2.1.1 方法一

将$U_{i, j}^n = \xi _{p, q}^n{e^{{\bf{i}}\left( {\alpha {x_i} + \beta {y_j}} \right)}}$代入(8)式, 其中i为格式所在的空间部分-x的层数, j为格式所在的空间部分-y的层数, n为格式所在的时间层数; p, q∈ $\mathbb{N} $; i为虚数单位; $\alpha = \frac{{2p\pi }}{l};\beta = \frac{{2q\pi }}{l}$。等式两边消去公因子e^{i(αx_i+βy_j)}, 得到

$ \begin{array}{l} \xi _{p, q}^{n + 1} - 2\xi _{p, q}^n + \xi _{p, q}^{n - 1}\\ = {r^2}\left[ {\theta \left( {\xi _{p, q}^{n + 1}{e^{{\bf{i}}\alpha h}} - 2\xi _{p, q}^{n + 1} + \xi _{p, q}^{n + 1}{e^{ - {\bf{i}}\alpha h}}} \right) + (1 - 2\theta )\left( {\xi _{p, q}^n{e^{{\bf{i}}\alpha h}} - 2\xi _{p, q}^n + \xi _{p, q}^n{e^{ - {\bf{i}}\alpha h}}} \right) + \theta \left( {\xi _{p, q}^{n - 1}{e^{{\bf{i}}\alpha h}} - 2\xi _{p, q}^{n - 1} + \xi _{p, q}^{n - 1}{e^{ - {\bf{i}}\alpha h}}} \right)} \right]\\ + {r^2}\left[ {\theta \left( {\xi _{p, q}^{n + 1}{e^{{\bf{i}}\beta h}} - 2\xi _{p, q}^{n + 1} + \xi _{p, q}^{n + 1}{e^{ - {\bf{i}}\beta h}}} \right) + (1 - 2\theta )\left( {\xi _{p, q}^n{e^{{\bf{i}}\beta h}} - 2\xi _{p, q}^n + \xi _{p, q}^n{e^{ - {\bf{i}}\beta h}}} \right) + \theta \left( {\xi _{p, q}^{n - 1}{e^{{\bf{i}}\beta h}} - 2\xi _{p, q}^{n - 1} + \xi _{p, q}^{n - 1}{e^{ - {\bf{i}}\beta h}}} \right)} \right]。\end{array} $

(9)

令ξ_{p, q}ⁿ⁺¹=λ², ξ_{p, q}ⁿ=λ, ξ_{p, q}^n-1=1, 结合欧拉公式e^ix=cos x+ isin x, 则(9)式可化简为

$ \begin{array}{l} {\lambda ^2}\left[ {1 - 2{r^2}\theta ((\cos \alpha h - 1) + (\cos \beta h - 1))} \right] - \lambda \left[ {2 + 2{r^2}(1 - 2\theta )((\cos \alpha h - 1) + (\cos \beta h - 1))} \right]\\ + 1 - 2{r^2}\theta ((\cos \alpha h - 1) + (\cos \beta h - 1)) = 0。\end{array} $

(10)

由于r²>0, θ>0, β=2πη, η∈ $\mathbb{N} $, 步长h>0.又cos αh-1≤0; cos βh-1≤0, 所以

$ 1 - 2{r^2}\theta [(\cos \alpha h - 1) + (\cos \beta h - 1)] > 0。$

(11)

(10) 式两边消去公因子1-2r²θ[(cos αh-1)+(cos βh-1)], 可得特征方程为

$ {\lambda ^2} - \frac{{2 + 2{r^2}(1 - 2\theta )[(\cos \alpha h - 1) + (\cos \beta h - 1)]}}{{1 - 2{r^2}\theta [(\cos \alpha h - 1) + (\cos \beta h - 1)]}}\lambda + 1 = 0。$

(12)

2.1.2 方法二

由于差分格式(8)是一个三层格式, 在研究其稳定性之前, 考虑先将(8)式转化为如下的两层差分方程组

$ \left\{ {\begin{array}{*{20}{l}} {\left( {1 + 4{r^2}} \right)U_{i, j}^{n + 1} - {r^2}\theta \left( {U_{i + 1, j}^{n + 1} + U_{i - 1, j}^{n + 1}} \right) - {r^2}\theta \left( {U_{i, j + 1}^{n + 1} + U_{i, j - 1}^{n + 1}} \right)}\\ { = \left( {2 - 4{r^2}(1 - 2\theta )} \right)U_{i, j}^n + {r^2}(1 - 2\theta )\left( {U_{i + 1, j}^n + U_{i - 1, j}^n} \right) + {r^2}(1 - 2\theta )\left( {U_{i, j + 1}^n + U_{i, j - 1}^n} \right) + V_j^n}\\ {V_j^{n + 1} = - \left( {1 + 4{r^2}\theta } \right){r^2}U_{i, j}^n + {r^2}\theta \left( {U_{i + 1, j}^n + U_{i - 1, j}^n} \right) + {r^2}\theta \left( {U_{i, j + 1}^n + U_{i, j - 1}^n} \right)} \end{array}} \right.。$

(13)

由分离变量法, 将U_{i, j}ⁿ=ξ_{p, q}ⁿe^{i(αx_i+βy_j)}, V_jⁿ=γ_{p, q}ⁿe^{i(αx_i+βy_j)}代入上述差分方程组(13), 结合欧拉公式, 化简整理得

$ \left\{ {\begin{array}{*{20}{l}} {\xi _\eta ^{n + 1} = \frac{{2 - 4{r^2}(1 - 2\theta ) + 2{r^2}(1 - 2\theta )(\cos \alpha h + \cos \beta h)}}{{1 + 4{r^2}\theta - 2{r^2}\theta (\cos \alpha h + \cos \beta h)}}\xi _\eta ^n + \frac{1}{{1 + 4{r^2}\theta - 2{r^2}\theta (\cos \alpha h + \cos \beta h)}}\gamma _\eta ^n}\\ {\gamma _\eta ^{n + 1} = \left[ {2{r^2}\theta (\cos \alpha h + \cos \beta h) - 1 - 4{r^2}\theta } \right]\xi _\eta ^n} \end{array}} \right.。$

(14)

上式可写为如下的矩阵形式

$ \left( {\begin{array}{*{20}{c}} {\xi _\eta ^{n + 1}}\\ {\gamma _\eta ^{n + 1}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {\frac{{2 - 4{r^2}(1 - 2\theta ) + 2{r^2}(1 - 2\theta )(\cos \alpha h + \cos \beta h)}}{{1 + 4{r^2}\theta - 2{r^2}\theta (\cos \alpha h + \cos \beta h)}}}&{\frac{1}{{1 + 4{r^2}\theta - 2{r^2}\theta (\cos \alpha h + \cos \beta h)}}}\\ {2{r^2}\theta (\cos \alpha h + \cos \beta h) - 1 - 4{r^2}\theta }&0 \end{array}} \right)\left( {\begin{array}{*{20}{l}} {\xi _\eta ^n}\\ {\gamma _\eta ^n} \end{array}} \right), $

(15)

即差分方程的增长矩阵为

$ G(\beta , k) = \left( {\begin{array}{*{20}{c}} {\frac{{2 - 4{r^2}(1 - 2\theta ) + 2{r^2}(1 - 2\theta )(\cos \alpha h + \cos \beta h)}}{{1 + 4{r^2}\theta - 2{r^2}\theta (\cos \alpha h + \cos \beta h)}}}&{\frac{1}{{1 + 4{r^2}\theta - 2{r^2}\theta (\cos \alpha h + \cos \beta h)}}}\\ {2{r^2}\theta (\cos \alpha h + \cos \beta h) - 1 - 4{r^2}\theta }&0 \end{array}} \right)。$

(16)

该矩阵的特征多项式即为差分格式的特征方程(12)式。

上述两种方法相比较, 方法一在计算上更为直接, 方法二则突出了格式的构造特点。

2.2 差分格式的稳定性

选取不同的差商近似偏导数, 得到的差分格式是不同的, 为了衡量差分格式的好坏, 还需分析格式的稳定性^[12]。为了便于后面对差分格式稳定性进行讨论, 首先介绍引理1和引理2, 然后把θ所属的定义区间分为多段, 在每一段上分类讨论格式稳定时对应网格比r的取值范围。

引理1^[2]在求解区域Ω={(x, t)|0 < x < l, 0 < t≤T}内, 差分格式按谱范数稳定(或平方稳定)的充要条件是对于任意的0 < k < k₀, 0 < nk < T和所有的β=2πη, η∈ $\mathbb{N} $, 均有如下式子成立

$ |G(\beta , k)| \le 1 + O(k)。$

有了引理1, 在判别差分格式稳定性时, 只需先计算增长因子或增长矩阵, 并给出使引理1成立的条件。由于在2.1节中, 利用两种方法计算得到了差分格式的特征方程, 为了简化计算, 下面引入引理2。

引理2^[2]实系数二次方程λ²-bλ-c=0的根按模大于1的充要条件为|b|≤1-c≤2。

由引理2得, 方程(12)的根按模不大于1的充要条件为

$ \left| {\frac{{2 + 2{r^2}(1 - 2\theta )[(\cos \alpha h - 1) + (\cos \beta h - 1)]}}{{1 - 2{r^2}\theta [(\cos \alpha h - 1) + (\cos \beta h - 1)]}}} \right| \le 2。$

(17)

由(11)式知1-2r²θ[(cos αh-1)+(cos βh-1)]>0, 所以(17)式等价于

$ \left| {2 + 2{r^2}(1 - 2\theta )[(\cos \alpha h - 1) + (\cos \beta h - 1)]} \right| \le 2 - 4{r^2}\theta [(\cos \alpha h - 1) + (\cos \beta h - 1)]。$

(18)

至此, 差分格式(8)的稳定性问题转化为了求解不等式(18)。由于cos αh≤1;cos βh≤1, 所以(cos αh-1)+(cos βh-1)≤0恒成立, 即(18)式左端绝对值符号内的正负性由1-2θ(0≤θ≤1)所决定, 因此可以将θ= $\frac{1}{2}$作为分类讨论时的一个间断点。在下面的分析中, 我们发现θ= $\frac{1}{4}$也是一个特殊的点, 因此将θ= $\frac{1}{4}$单独考虑。综上, θ所属区间[0, 1]便被分为四段[0, $\frac{1}{4}$), $\frac{1}{4}$, ($\frac{1}{4}$, $\frac{1}{2}$), [$\frac{1}{2}$, 1], 如图 1所示。

下面将逐一对四种情况进行分析讨论。特别说明的是, 当θ=0时, 格式(8)便转化为显式差分格式, 显式差分格式稳定时所对应的网格比的取值范围为0 < r≤1, 这里不再对其稳定性做详细的讨论^[2]。

(1) 当0 < θ < $\frac{1}{4}$时格式的稳定性

当0 < θ < $\frac{1}{4}$时, (18)式左端项绝对值内的正负并不确定, 将其转化为

$ - 2 - 2{r^2}(1 - 2\theta )[(\cos \alpha h - 1) + (\cos \beta h - 1)] \le 12 + 2{r^2}(1 - 2\theta )[(\cos \alpha h - 1) + (\cos \beta h - 1)]|。$

结合(18)式可得

$ - 2 - 2{r^2}(1 - 2\theta )[(\cos \alpha h - 1) + (\cos \beta h - 1)] \le 2 - 4{r^2}\theta [(\cos \alpha h - 1) + (\cos \beta h - 1)], $

化简整理得

$ {r^2}((1 - \cos \alpha h) + (1 - \cos \beta h))\left( {\frac{{1 - 4\theta }}{2}} \right) \le 1 $

(19)

此时1-4θ>0, 即

$ {r^2}((1 - \cos \alpha h) + (1 - \cos \beta h)) \le \frac{2}{{1 - 4\theta }}。$

(20)

由于(1-cos αh)+(1-cos βh)最大值为4, 结合(20)式得

$ {r^2} \le \frac{1}{{2(1 - 4\theta )}}。$

又因为网格比r>0, 所以当0 < θ < $\frac{1}{4}$时, 差分格式(8)稳定的充要条件是r≤ $\frac{1}{{\sqrt {2(1 - 4\theta )} }}$。

(2) 当θ= $\frac{1}{4}$时格式的稳定性

将θ= $\frac{1}{4}$代入(12)式, 计算得到对应的特征方程为

$ {\lambda ^2} - \frac{{2 - 2{r^2}{{\sin }^2}\frac{{\alpha h}}{2} - 2{r^2}{{\sin }^2}\frac{{\beta h}}{2}}}{{1 + {r^2}{{\sin }^2}\frac{{\alpha h}}{2} + {r^2}{{\sin }^2}\frac{{\beta h}}{2}}}\lambda + 1 = 0, $

(21)

根据引理2, 方程(21)的根按模不大于1的充要条件为

$ \frac{{\left| {2 - 2{r^2}\left( {{{\sin }^2}\frac{{\alpha h}}{2} + {{\sin }^2}\frac{{\beta h}}{2}} \right)} \right|}}{{\left| {1 + {r^2}\left( {{{\sin }^2}\frac{{\alpha h}}{2} + {{\sin }^2}\frac{{\beta h}}{2}} \right)} \right|}} \le 2, $

上式显然成立。即θ= $\frac{1}{4}$时, 格式(8)是无条件稳定的, 与r的取值无关。

(3) 当$\frac{1}{4} < \theta < \frac{1}{2}$时格式的稳定性

当$\frac{1}{4} < \theta < \frac{1}{2}$时, (18)式左边绝对值内的正负仍不能确定, 同理可经化简整理得到(19)式。此时1-4θ < 0, 不等式(19)左边小于或等于0, 即不等式(19)恒成立。所以当$\frac{1}{4} < \theta < \frac{1}{2}$时, 对于任意r>0, 格式(8)是无条件稳定的。

(4) 当$\frac{1}{2}$ ≤θ≤1时格式的稳定性

当$\frac{1}{2}$≤θ≤1时, 由于2+2r²(1-2θ)[(cos αh-1)+(cos βh-1)]>0.则(18)式可化为

$ 2 + 2{r^2}(1 - 2\theta )[(\cos \alpha h - 1) + (\cos \beta h - 1)] \le 2 - 4{r^2}\theta [(\cos \alpha h - 1) + (\cos \beta h - 1)], $

(22)

上式两边同时减去2, 并除以2r²得

$ (1 - 2\theta )[(\cos \alpha h - 1) + (\cos \beta h - 1)] \le - 2\theta [(\cos \alpha h - 1) + (\cos \beta h - 1)], $

(23)

又

$ (\cos \alpha h - 1) + (\cos \beta h - 1) \le 0, $

当(cos αh-1)+(cos βh-1)=0时, (23)式恒成立。当(cos αh-1)+(cos βh-1)≠0时, 即

$ (1 - 2\theta ) \ge - 2\theta {\rm{, }} $

上式恒成立。对任意$\frac{1}{2}$ ≤θ≤1, (22)式恒成立, 即对任意r, 差分格式(8)是无条件稳定的。

综上, 便得到二维二阶线性双曲型方程隐式差分格式的稳定性条件, 见表 1。

3 结论

本文通过分离变量法讨论了二维二阶双曲型方程交替方向隐式差分格式的稳定性, 结果表明在不同的权重θ下, 差分格式的稳定性与网格比r的取值有关。当θ=0时, 格式便转化为了显式差分格式, 此时格式稳定的条件为0 < r≤1;当0 < θ < 1/4时, 差分格式稳定的充要条件为$0 < r \le 1/\sqrt {2(1 - 4\theta )} $; 当1/4≤θ < 1时, 差分格式是无条件稳定的。