研究报告有限差分格式稳定性的其他方法-报告

研究有限差分格式稳定性的其他方法

摘要

偏微分方程的求解一直是大家比较关心的一个问题，而有限差分格式则是求解偏微分方程时常用并且有效的一个方法。因此，研究有限差分格式的性质就显得尤为重要。在课上我们已经跟着老师学习了运用Fourier方法研究有限差分格式的稳定性，但是在很多研究有限差分格式稳定性的问题中仅仅会用Fourier方法是不够的，所以在本篇论文中，将会介绍其他三种常用的研究有限差分格式稳定性的方法，分别是：Hirt启示型方法、直接方法（或称矩阵方法）和能量不等式方法。

关键字：偏微分方程；有限差分格式；稳定性

Abstract

The solution of partial differential equations has been more concerned with a problem, and the finite difference scheme is a mon and effective method for solving partial differential equations. Therefore, it is very important to study the character of the finite difference scheme. We have followed the teacher to learn the use of Fourier method of finite difference scheme stability, but in a lot of research on the stability of finite difference scheme is only used Fourier method is not enough, so in this paper, will introduce the other three kinds of monly used in the study of finite difference scheme stability method, respectively is: Hirt enlightenment method, direct method (or matrix method) and energy inequality method.

Key words: partial differential equation; finite difference scheme; stability

1 前言

微分方程的定解问题就是在满足某些定解条件下求微分方程的解。在空间区域的边界上要满足的定解条件称为边值条件。如果问题与时间有关，在初始时刻所要满足的定解条件，称为初值条件。不含时间而只带边值条件的定解问题，称为边值问题。与时间有关而只带初值条件的定解问题，称为初值问题。同时带有两种定解条件的问题，称为初值边值混合问题。定解问题往往不具有解析解，或者其解析解不易计算。所以要采用可行的数值解法。有限差分方法就是一种数值解法，它的基本思想是先把问题的定义域进行网格剖分，然后在网格点上，按适当的数值微分公式把定解问题中的微商换成差商，从而把原问题离散化为差分格式，进而求出数值解。此外，还要研究差分格式的解的存在性和唯一性、解的求法、解法的数值稳定性、差分格式的解与原定解问题的真解的误差估计、差分格式的解当网格大小趋于零时是否趋于真解（即收敛性），等等。有限差分方法具有简单、灵活以及通用性强等特点，容易在计算机上实现。在课上我们已经跟着老师学习了运用Fourier 方法研究有限差分格式的稳定性，但是在很多研究有限差分格式稳定性的问题中仅仅会用Fourier 方法是不够的，所以在本篇论文中，将会介绍其他三种常用的研究有限差分格式稳定性的方法，分别是：Hirt 启示型方法、直接方法和能量不等式方法。

2 Hirt 启示性方法

2.1 方法概述

Hirt 启示性方法是一种近似分析方法。主要是把差分格式在某确定点上作泰勒级数近似展开，把高阶误差略去，只留下最低阶的误差项。如果差分格式是相容的，那么这样得到的新的微分方程（称之为第一微分近似或修正微分方程）与原来的微分方程相比只增加了一些含小参数的较高阶导数的附加项。Hirt 方法就是利用第一微分近似的适应性来研究差分格式的稳定性。Hirt 方法的判别准则是这样的：如果第一微分近似是适定的，那么原来微分方程的差分格式是稳定的，否则不稳定。其实所述的微分格式是原来微分方程问题的相容的差分格式，那么也可以看作第一微分近似问题的相容的差分格式。如果第一微分近似问题是不适

定的，那么它的差分格式将不稳定[1]

。

2.2 操作方法

先给出几个方程

0,,0,0>∈>=??+??t R x a x

u a t u （2.1） ,2,1,0,,2,1,0,011=±±==-+-++n j h

u u a

u u n

j

n j n j

n j τ

（2.2）

01

1=-+--+h

u u a

u u n j n j n j

n j τ

（2.3）

考虑对流方程（2.1）的差分格式（2.3），在点）

，（n j t x 进行Taylor 技术展开，有 )(][2][)

,(),(2221h O x

u h x u h t u u t x u n

j n j n j n j +??-??=--

)(][2][)

,(),(2221ττ

O t

u h t u t u u t x u n

j n j n j n j +??-??=-+ 利用对流方程（2.1），有

22

222)(x

u a x u a t t u ??=??-??=?? 因此，在点）

，（n j t x 上，有差分方程（2.3）可以得到 )(2222222h O x

u a ah x u a t u ++??-=??+??ττ）（ 略去高阶误差项，得出第一微分方程近似

2

2222x u

a ah x u a t u ??-=??+??）（τ 要使上面的抛物型方程有意义，必须有

02

22>-τ

a ah 而上面的不等号改为等号，则就化为原来的对流方程。在这两种情况下，相应的问题是适定

的。即第一微分近似适定的条件是

02

22≥-τa ah 由此得出差分格式（2.3）的稳定性条件是1≤λa ，其中h

τ

λ=

。此结论与Fourier 方法分析

得到的结论是一致的。

下面我们再来分析逼近对流方程（2.1）（仍设0>a ）的差分格式（2.2）的稳定性。模仿上面的推导可以得到它的第一微分近似是

22222x

u a ah x u a t u ??+-=??+??）（τ 可以看出22x

u

??的系数小于0，因此第一微分近似是不适定的，从而推出差分格式（2.2）是

不稳定的。

3 直接方法

关于抛物型方程初值问题的差分格式的稳定性问题，可以用直接方法（或称矩阵方法）来研究。下面用具体例子来说明这个方法的基本思想及使用方法。

考虑常系数扩散方程的初值问题

???

???

?>==∈=>∈>??=??0,0),(),0(),0(),()0,(0

),,0(,0,022t t l u t u l x x u x u t l x a x u

a t

u （3.1） 采用显示差分格式来逼近，即

???

?

???>==-==-=>+-=--++0

,01

,2,1),(1,2,1,0,200002

1

11n u u J j x u u J j n h u u u a u u J n j j n j n j n j n j n j τ （3.2） 其中l Jh =。先把差分格式（3.2）写成

1,2,1

,)21(111-=+-+=-++J j u a u a u a u n j n j n j n j λλλ （3.3） 其中2

h

τ

λ=

。可以把（3.3）写成向量形式，即

]00

[]][21212121[][0

012

2111

1

2

1

2

1

1J

n

n J n

J n n n J n J n n u u a u u u u a a a a a a a a a a u u u u

λλλλλλλλλλλ+----=--+-+-++ （3.4） 如果令

T

n J n n n u u u u ),,,(121-=

并考虑到00==n

J n u u ，则（3.4）式可以写成

n n Au u =+1 （3.5）

其中

]21212121[

λ

λλλ

λλ

λλλ

λa a a a a a a a a a A ----=

（3.6）

从显示格式出发，得到方程组（3.5）式，也可以理解为较为一般的形式，即对于逼近初值

问题（3.2）的其他二层格式也可以化为（3.5）式的形式。当然此时A 不是（3.6）式所表示的形式。如果差分格式是二层隐式格式。则A 为C B 1

-这种形式。因此（3.5）式这种形式可

理解为既包含二层显示格式又包含二层隐士格式的较为一般的形式。

引入误差向量~

n

n

n

u u z -=，其中n

u 是差分方程（3.5）的精确值（理论值），~n

u 是差

分方程（3.5）经数值求解得到的值（包括了舍入误差等）。显然，n

z 满足

n n Az z =+1 （3.7）

从而推出

0z A z n n = （3.8）

差分格式（3.5）的稳定性就要求

0,≥≤n K z n （3.9）

其中?为向量的2-数。由于

02

z A z n n ?≤

因此（3.9）式成立的充分必要条件为

M A n

≤2

（3.10）

上述采用2-数，当然也可以采用其他类型的数。对于稳定性条件（3.10），可以仿Fourier 方法中的推导，得到一些结论：

（1）谱半径条件

τρM A +≤1)( （3.11）

是差分格式稳定的一个必要条件，其中M 为常数。

（2）如果矩阵A 是一个正规矩阵，则（3.11）式也是格式稳定的一个充分条件。 下面讨论差分格式（3.5），（3.6）的稳定性。矩阵（3.6）是对称矩阵，所以只要使条件（3.11）成立即可。现在来计算A 的特征值。

令)1(-J 阶方阵

]0

1

101

1

110[

=S 则A 可以表示为

S a I a A λλ+-=)21(

其中I 为)1(-J 阶单位矩阵。由此可知，关键是求出S 的特征值和特征向量。

设γ和T J w w w w ),,,(121-= 分别为S 的特征值和特征向量，

w Sw γ=

写成分量的形式有

?

?

?==-==+-++02

,,1,0,0021J j j j w w J j w w w γ （3.12） 先求出j w ，再求出S 的特征值γ。由于S 为对称矩阵，所以其特征值γ为实数。由Gerschgorin 定理知，

∑-≠≤-1

J k

j kj kk s s γ

其中kj s 为矩阵S 的元素。由此得到2≤γ。（3.12）式的第一式为常系数线性差分方程。设其解具有如下形式：

0,≠=μμj j w

将它代入（3.12）式的第一式，便得到关于μ的一元二次方程

012=+-λμμ

此方程称为（3.12）式的第一式的特征方程。由于2≤γ，所以其解为

2)2

(12γ

γ

μ-±=

i 其中1-=

i 。可以看到

1)2

(1)2

(222

=-+=γγμ

取2)2

(1sin ,2cos λ

?γ

?-==

，则?μi e ±=。因此差分方程（3.12）的解可以表示为 J j e a e a w ij ij j ,,1,0,21 =+=-??

由00=w ，得到021=+a a 。再由0=J w ，得到021=+-?

?iJ iJ e

a e a ，从而有 02=--）（??

iJ iJ e e

a 由此可推0sin 2=?J a 。02≠a ，有1,2,1,-==J k k J π?。所以得到π?J

k

=，可以得到πγJ

k

k cos 2=。注意到J h 1=，则S 的特征值为πγkh k cos 2=。从而得到A 的特征

值为

1,2,1,2

sin 41cos 2212

-=-=+-=J k kh a kh a a k π

λπλλξ 当21≤λa 时，1)(≤A ρ。因此显示格式的稳定性条件为2

1

≤λa 。

下面讨论隐式格式

???

?

???===+-=-+-++++0

)

(20002

11

1111n J n j j n j n j n j n j n j u u x u u h u u u a u u τ 的稳定性。

可以把隐式格式写成向量形式

n n u B u 11-+=

其中T

n J n n n u u u u ),,,(121-= ，S a I a B λλ-+=)21(。

利用前面已经求得的S 的特征值，可以得到B 的特征值

1,,2,1),cos 1(21)(21)(-=-+=-+=J j kh a S a a B k k πλλμλμ

由此可知，1)(>B k μ，从而有1)(1≤-B k μ。注意的B 为对称矩阵，所以1

-B 也为对称矩阵，

利用直接方法结论（2）知，扩散方程隐式格式是无条件稳定的。 从上面的叙述看来，利用直接方法来分析抛物型方程的初值问题的差分格式并不困难。但在实际应用中却存在着一定的限制。上面讨论稳定性的两个例子中式依据了特殊矩阵S 才求出了)1(-J 阶矩阵A 、1

-B 的特征值。一般说来，计算高阶矩阵的特征值是相当困难的，因

此直接方法应用也就很困难了。

4 能量不等式方法

4.1 方法概述

在讨论线性常系数差分格式的稳定性问题时，建立了判别差分格式的稳定性准则，从而比较容易地判断一些差分格式的稳定性。但对于变系数问题和非线性问题，一般不能采用Fourier 方法和直觉法来讨论差分格式的稳定性。而对于上述这些问题，能量不等式方法是研究差分格式稳定性的有力工具。用能量不等式方法讨论差分格式稳定性是从稳定性的定义出发，通过一系列估计式来完成的。这个方法是偏微分方程中常用的能量方法的离散模拟，在此我们仅通过例子叙述其基本思想。

4.2 操作方法

考虑变系数对流方程的初值问题

??

???∈=≤<∈=??+??R x x g x u T

t R x x

u t x a t u

),()0,(0,,0),( （4.1） 假定0),(≥t x a ，建立差分格式

??

???==-+--+)(0

011j j n j n j n j n j n j x g u h

u u a u u τ （4.2） 其中),(n j n

j t x a a =。下面用能量不等式方法来讨论这个差分格式的稳定性。先把它改变形

式为

)(11n

j n j n j n j n j u u a u u -+--=λ

其中h

τ

λ=

为网格比。用1

+n j u 乘上式的两边，得

1

121)]([+-+--=n j n j n j n j n j n j u u u a u u λ）（

如果λ满足条件

1)(max ≤λn j j

a （4.3）

则有

])()[(2

])()[(2

12

1212

122

1+-++++

+-≤

n j n j n j n j n j n j n j u u a u u a u λλ）（

移项得

2

1222

1)(2

)()(n j n j n j n j

n j n j n j u a a u a u u -++

-≤λλλ）（

用h 乘上面不等式的两边，并对j 求和，令

∑∞

-∞

==

j n j

h

n h u

u

2

2)(

则有

∑∞

-∞

=++-+≤j n j n j n j h

n h

n h u a a u

u

2

122

1)(）（λ

如果

c x

a

t

x ≤??sup

, （4.4） 则有

2

21

)1(h

n h

n u c u τ+≤+

由此可得

T n u e

u c u h

cT

h

h

n ≤≤+≤ττ,)1(20202

由于问题是线性的，因此上述不等式就证明了差分格式（4.2）的稳定性。由此看出，条件

（4.4）是微分方程问题中给定的。而差分格式稳定性条件就是（4.3式）式。如果a t x a =),(即为常系数问题，那么（4.4）式满足，而条件（4.3）就化为1≤λa ，这与我们在课上所学的用Fourier 方法得到的结论一致。

5 结论

在本篇论文中，从微分方程的基本概念出发，先介绍了微分方程中比较基本的概念，然后又介绍了有限差分格式的性质。在介绍有限差分格式时从三种求解有限差分格式稳定性的方法出发，分别是：Hirt 启示性方法、直接方法（或矩阵方法）和能量不等式方法。在介绍这三种的方法时也是先从基本思想出发，然后分别阐述其方法原理、公式推导和实际应用等。但是求解有限差分格式稳定性的方法很多，作者也仅仅介绍了三种方法，希望能起到抛砖引玉的作用。

参考文献

[1] 陆金甫, 关治：《偏微分方程数值解法》，清华大学，，2003 [2] 康等编：《数值计算方法》，国防工业，，1978. [3] 胡祖炽编：《计算方法》，高等教育，，1959。 [4] 清华大学、大学《计算方法》编与组编：（计算方法），科学，，1980。 [5] 朱幼兰等著：《初边值问题差分法及绕流》，科学，，1980。 [6] R ．D ．里奇特迈尔著，何旭初等译：《初值问题差分方法》，科学，，1966

[7 ] R. D. Richtmyer ，Difference Methods for Initial-Value Problems ，Interscience Pub.，New York ，

1957.）

[8] R. D. Richtmyer ，K. W. Morton ，Difference Methods for Initial-Value Problems ，2nd ed.，

Interscience Pub.，New York ，1967.