对流扩散方程有限差分方法

对流扩散方程有限差分方法

求解对流扩散方程的差分格式有很多种，在本节中将介绍以下3种有限差分格式：中心差分格式、Samarskii格式、Crank-Nicolson型隐式差分格式。

3.1 中心差分格式

时间导数用向前差商、空间导数用中心差商来逼近，那么就得到了（1）式的中心差分格式[6]

1ununjj若令 aanunj1uj12hvnnunj12ujuj1h2 （3）

h，v1 unj，则（3）式可改写为

h21nnnnnun(uu)(u2uujj1j1j1jj1) （4）

21 从上式我们看到，在新的时间层n1上只包含了一个未知量unj，它可以nn由时间层n上的值unj1，uj，uj1直接计算出来。因此，中心差分格式是求解对

流扩散方程的显示格式。

1假定u(x,t)是定解问题的充分光滑的解，将unununj，j1，j1分别在(xj,tn)处

进行Taylor展开：

1unjuu(xj,tn1)u(xj,tn)O(2)

tjh22uu3u(xj1,tn)u(xj,tn)h2O(h)

xj2xjh22uu3u(xj1,tn)u(xj,tn)h2O(h)

xj2xj1ununjjnunj1uj1nnunj12ujuj1nnnunj1nnunj1代入(4)式，有 T(xj,tn)na2hnvh2

n2uuu2O()aaO(h)v2vO(h2) tjxjxj2uuuav2O()(av)O(h2) tjxjxj . 专业资料分享 .

nnn 完美WORD格式.整理

O(h2)

显然，当0，h0时，T(xj,tn)0，即中心差分格式与定解问题是相容的。由以上的讨论也可得知，对流扩散方程的中心差分格式的截断误差为

O(h2)。

对于我们上面构造的差分格式，是否可以直接用于实际计算呢？也就是说，如果初始值有误差，在计算过程中误差会不会扩大传播呢？这就是接下来我们要讨论的是差分方程的稳定性问题。下面用Fourier方法来分析中心差分格式的稳定性。 令ujnvneikjh，代入到（4）式

1vn1eikjhvneikjh(vneik(j1)hvneik(j1)h)(vneik(j1)h2vneikjhvneik(j1)h)2n1[12(1coskh)isinkh]vn 整理得 v所以该差分格式的增长因子为：

G(,k)12(1coskh)isinkh 其模的平方为

G(,k)[12(1coskh)]22(sinkh)2

kh)42(1coskh)22(sikhn)2 14(1cos 1(1coskh)[442(1coskh)2(1coskh)] 由于1coskh0，所以G(,k)1（即差分格式稳定）的充分条件为

kh)2(1coskh)0 442(1cos上式可以改写为

(2282)1cokhs4220 21注意到(1coskh)[0,1]，所以上面不等式满足的条件为

2 (2282)4220， 4220。

由此得到差分格式（3）的稳定性限制为

12v 2 ， v2。

2ha故有结论：对流扩散方程的中心差分格式是条件稳定的。根据Lax等价定理，我们可以知道，对流扩散方程的中心差分格式是条件收敛的。

. 专业资料分享 .

完美WORD格式.整理

3.2 Samarskii格式

设a>0，先对方程(1)作扰动，得到另一个对流扩散方程[7]

uu12uav （5） tx1Rx2其中 R1（5）式化为（1）式 ha，当h0时，

2v对于（5）式，构造迎风格式

1nunujjanunujj1h1v1Rnnun2uuj1jj1h2 （6）

差分格式（6）称为逼近对流扩散方程的Samarskii格式。

首先推导（6）的截断误差。设u(x,t)是对流扩散方程（1）式的充分光滑的解

Tjnu(xj,tn1)u(xj,tn)au(xj,tn)u(xj1,tn)hu(xj1,tn)2u(xj,tn)u(xj1,tn)1v1Rh2令

nju(xj,tn1)u(xj,tn)vu(xj1,tn)2u(xj,tn)u(xj1,tn)h2

用 Taylor级数展开有

un2un2 ()jv(2)jO(h)

txnj再令 anju(xj,tn)u(xj1,tn)hu(xj1,tn)2u(xj,tn)u(xj1,tn)1(1)v 1Rh2用 Taylor级数展开有

unah2un12unv(2)jO(h2) a()j(2)jt2x1Rxnjun2unR2uv(2)O(h2) a()jRv(2)jx1RxxunR22unv(2)jO(h2) a()jx1Rx . 专业资料分享 .

完美WORD格式.整理

2R由于

1Rh2a224v(1ha)2v22ha(4v2vh)a2O(h2)

所以 jna(njnjun)jO(h2) xnjunun2un T()ja()jv(2)jO(h2)

txx O(h2)

当0，h0时，Tjn0，所以Samarskii格式与定解问题是相容的，并且其截断误差为O(h2)。

现在看看Samarskii格式的稳定性。将（6）式两边同时加上

ann(uj12unjuj1)，把（6）式化为 2h 令 v1ununjjanunj1uj12hnnnvahuj12ujuj1() 1R2h2vah，则上式即为： 1R21ununjjanunj1uj12hvnnunj12ujuj1h2

根据中心显示格式稳定性的讨论，可以得到（6）式的稳定性条件为 1 v2 ，

2h2v 2

a即

vah1 v()2 ，

1R2h22vah 2()

a1R2稳定性的第二个条件等价于

2vha2(1R)a1

. 专业资料分享 .

完美WORD格式.整理

2vah(1R)2a2(1R)2()1R2v2vh而 2vhah(1R)ah(1R)112v2va2(1R)a利用不等式

ah(1R)2)ah(1R)2v1

ah(1R)2v12v所以

(

2vha2(1R)a2vah(1R)1ah2(1)2v()2 1Rh2v1R2vh利用稳定性的第一个条件，有

2vha2(1R)a211，从而可知稳定性条件的2第二个条件可由第一个条件推出，因此差分格式的稳定性条件为

1ah1 v()2，

1R2vh2v1ah即 。 2221ah/2vh 由Lax等价定理可知，Samarskii格式也是条件收敛的。

3.3 Crank-Nicolson型隐式差分格式

前面讨论了求解对流扩散方程的两种显示格式，它们都是条件稳定的，为

了放松稳定性条件，可以采用隐式格式进行求解。现在考虑Crank-Nicolson型隐式差分格式[8]

1nunujjnnn1n1auj1uj1uj1uj1() 22h2h

nnnn1n1n1vuj12ujuj1uj12ujuj1) （7） (222hh令av，2，则（7）式可化为 hh1n1n1(2)un4(1)u(2)uj1jj1

nn(2)un4(1)u(2)uj1jj1 （8）

. 专业资料分享 .

完美WORD格式.整理

把（8）式用矩阵的形式

24(1)(2)4(1)2 (2)4(1)2(2)4(1)4(1)224(1)= 22 4(1)224(1)n1u2n1u3 n1uJ2n1uJ1nu2nu3 nuJ2nuJ1(2)(u1n1u1n)0 （9） + 0n1n(2)(uuJJ)

24(1)(2)4(1)2 ， 设 A(2)4(1)2(2)4(1)4(1)224(1) B2 ， 4(1)224(1)2n(2)(u1n1u1n)u2n0u3 则有 Un ， FnuJ20nn1n(2)(uuuJJ)J1

AUn1BUnF

下面讨论Crank-Nicolson型格式的截断误差和精度。该格式涉及到时间层n和时间层n1上的xj1，xj，xj1处六个点。设u(x,t)是定解问题的充分光滑的解，

. 专业资料分享 .

完美WORD格式.整理

1n1n1nnn把(7)式中各un然后将unj的值用u(xj,tn)代替，j1，uj，uj1 ，uj1，uj，uj1分

别在点(xj,tn1/2)处进行Taylor展开[9]：

1ununjjun223un2()j()j t24t3un1h23un1()j(3)j x6x111111

1n1unj1uj12hun22un223un2h23un2)j()j()j ()j(x2xt8xt26x3nunj1uj12hun22un223un2h23un2()j()j()(3)j 2jx2xt8xt6x2un23un224un2h24un2(2)j(2)j(22)j()j

2xt8xt12x4x1111111111111n11ununj12ujj1h2nnunj12ujuj1h22un23un224un2h24un2(2)j(2)j(22)j(4)j

2xt8xt12xx这里出现的u的各阶偏导数假设都是存在而且连续的。于是（7）式的截断误差

Tjn23u()j24t3n12a23un2ah23un2v24un2vh24un2()j()j()j()j8xt26x38x2t212x41111 O(2h2)

显然，Crank-Nicolson型格式的精度是二阶的。再来看看该格式的稳定性情况，

我们还是用Fourier方法来分析。 令ujnvneikjh，代入到（8）式

(2)vn1eik(j1)h4(1)vn1eikjh(2)vn1eik(j1)h

(2)vneik(j1)h4(1)vneikjh(2)vneik(j1)h 整理得

[4(1)4coskhi2sinkh]vn1[4(1)4coskhi2sinkh]vn所以Crank-Nicolson型格式的增长因子是

(1coskh)i22sinkh G(,k)(1coskh)i

sinkh . 专业资料分享 .

完美WORD格式.整理

其模的平方

(1cokhs)2(sinkh)222 G(,k)

(1cokhs)2(sinkh)22改写上式

G(,k)124(1cokhs)(1cokhs)2(sinkh)22由于 1coskh0及上式的分母为正，故



G(,k)10

即 G(,k)1，从而得出Crank-Nicolson型格式是无条件稳定的。根据Lax 等价定理，Crank-Nicolson型格式也是无条件收敛的。

224、数值例子

给出如下对流扩散方程的初边值问题[10]：

uu2uav2 x[0,1] txx

u(0,t)e(av)t,u(1,t)e1(av)t t0

u(x,0)ex x[0,1]

所讨论的对流扩散方程的精确解为

x(av)tu(x,t)e x[0,1]

4.1 三种差分格式的比较

在各种对流扩散问题中，有许多对流相对于扩散来说在问题中起主导作用。

对流占有扩散问题的数值求解面临很多困难。因此，对流占有扩散问题的有效数值解法一直是计算数学中重要的研究内容[11]。

取a1，v0.001，h0.1，Tmax1，此时上面给出的就是一个对流占

优扩散问题。那么，本文讨论的三种差分格式对对流占有扩散问题的求解效果是怎样的呢？现在我们就来看看这个问题。

首先，根据差分格式的稳定性条件，确定的取值范围。

2v1（1）中心差分格式：根据稳定性条件 2，v2可知，要使中心差

2ah分格式稳定，的取值必须满足：0.002

. 专业资料分享 .

完美WORD格式.整理

v1ah（2）Samarskii格式：根据稳定性条件 可知，的取2221ah/2vh值必须满足：0.099

（3）Crank-Nicolson格式：该差分格式是无条件稳定的，所以可以取任意值。

要使三种差分格式都是稳定的，不妨取0.002。

首先，我们通过表格看看三种差分格式的数值解与准确解之间的相对误差。 表4.1 t0.2时三种差分格式结果的比较

x

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

x

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

. 中心差分格式 Samarskii格式 Crank-Nicolson格式

数值解 误差

数值解 误差 数值解 误差 1.2192 0

1.2192 0

1.2192

0

1.1009 -0.0023 1.1076 0.0044 1.1034 0.0002 0.9950 -0.0032 1.0051 0.0069 0.9985 0.0003 0.8999 -0.0033 0.9110 0.0078 0.9035 0.0003 0.8142 -0.0031 0.8250 0.0077 0.8175 0.0002 0.7367 -0.0028 0.7467 0.0072 0.7397 0.0002 0.6666 -0.0025 0.6757 0.0066 0.6693 0.0002 0.6031 -0.0023 0.6114 0.006

0.6056

0.0002 0.5459 -0.0019 0.5532 0.0054 0.5480 0.0002 0.4931 -0.0026 0.5006 0.0049 0.4959 0.0002 0.4485 0

0.4485 0

0.4485

0

表4.2 t0.4时三种差分格式结果的比较 中心差分格式 Samarskii格式 Crank-Nicolson格

式

数值解 误差 数值解 误差 数值解 误差 1.4894 0

1.4894 0

1.4894 0

1.3450 -0.0027 1.3537 0.006 1.3479

0.0002 1.2144 -0.0050 1.2301 0.0107 1.2199 0.0005 1.0969 -0.0065 1.1171 0.0137 1.1040 0.0006 0.9915 -0.0069 1.0137 0.0153 0.9990 0.0006 0.8966 -0.0068 0.9191 0.0157 0.9040 0.0006 0.8115 -0.0059 0.8327 0.0153 0.8179 0.0005 0.7333 -0.0063 0.7539 0.0143 0.7402 0.0006 0.6656 -0.0036 0.6824 0.0132 0.6696

0.0004 专业资料分享 .

准 确 解 1.2192

1.1032 0.9982 0.9032 0.8173 0.7395 0.6691 0.6054 0.5478 0.4957 0.4485

准

确 解 1.4894

1.3477 1.2194 1.1034 0.9984 0.9034 0.8174 0.7396 0.6692

完美WORD格式.整理

0.9 1.0

0.5982 -0.0074 0.6176 0.012 0.5479 0

0.5479 0

0.6062 0.5479

0.0006 0.6056 0

0.5479

表4.3 t0.6时三种差分格式的结果比较

x

数值解 误差

中心差分格式 Samarskii格式 数值解 误差 Crank-Nicolson格式 数值解 误差 准

确 解 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

. 1.8196 0

1.8196 1.6433 -0.0031 1.6538 1.4839 -0.0058 1.5032 1.3397 -0.0083 1.3660 1.2100 -0.0097 1.2410 1.0923 -0.0113 1.1269 0.9892 -0.0094 1.0226 0.8911 -0.0125 0.9274 0.8121 -0.0055 0.8404 0.7257 -0.0141 0.7612 0.6694 0

0.6694

专业资料分享 0

1.8196

0.0074 1.6467 0.0135 1.4902 0.018

1.3487

0.0213 1.2205 0.0233 1.1046 0.0240 0.9994 0.0238 0.9047 0.0228 0.8181 0.0214 0.7410 0

0.6694

0

0.0003 0.0005 0.0007 0.0008 0.001

0.0008 0.0011 0.0005 0.0012 0

.

1.8196

1.6464 1.4897 1.3480 1.2197 1.1036

0.9986 0.9036 0.8176 0.7398 0.6694

完美WORD格式.整理

x 0

表4.4 t0.8时三种差分格式结果的比较

中心差分格式 Samarskii格式 Crank-Nicolson格

式

数值解 误差 数值解 误差 数值解 误差 2.2229 0

2.2229 0

2.2229

0

准

确 解 2.2229

0.1 2.0073 -0.004 0.2 1.8132 -0.0067 0.3 1.6366 -0.0101 0.4 1.4794 -0.0106 0.5 1.3325 -0.0157 0.6 1.2087 -0.0112 0.7 1.0839 -0.0199 0.8 0.9915 -0.0073 0.9 0.8809 -0.0229 1.0 0.8177 0

2.0204 0.0091 2.0117 1.8364 0.0165 1.8205 1.6691 0.0224 1.6476 1.5169 0.0269 1.4909 1.3784 0.0302 1.3496 1.2521 0.0322 1.2209 1.1370 0.0332 1.1056 1.0318 0.033 0.9358 0.032 0.8177 0

0.9994 0.9057 0.8177

0.0004 2.0113 0.0006 1.8199 0.0009 1.6467 0.0009 1.4900 0.0014 1.3482 0.001

1.2199

0.0018 1.1038 0.0006 0.9988 0.0019 0.9038 0

0.8177

接下来，我们看看这三种差分格式在不同时间t的图形。

t=0.21.31.21.11中心格式Samarskii格式Crank-Nicolson格式准确解 t=0.41.51.41.31.21.1中心格式Samarskii格式Crank-Nicolson格式准确解 0.910.80.70.60.50.4 00.90.80.70.60.5 00.10.20.30.40.50.60.70.80.910.10.20.30.40.50.60.70.80.91

图4.1 t0.2时三种差分格式结果的比较 图4.2 t0.4时三种差分格式结果的比较

t=0.62中心格式Samarskii格式Crank-Nicolson格式准确解 2.62.42.221.41.81.61.411.20.81 00.8 0中心格式Samarskii格式Crank-Nicolson格式准确解t=0.8 1.81.61.20.10.20.30.40.50.60.70.80.910.10.20.30.40.50.60.70.80.91

. 专业资料分享 .

完美WORD格式.整理

图4.3 t0.6时三种差分格式结果的比较 图4.4 t0.8时三种差分格式结果的比较

4.2 结果分析

由表格中的数据以及图示可以看出，对于对流扩散方程的数值求解，三种差分格式的稳定性都比较好，其中以Crank-Nicolson格式的效果最好。

5.小结

对流扩散问题的数值求解一直是许多计算工作者比较重视的一类问题。本文分析了对流扩散方程的中心差分格式、Samarskii格式以及Crank-Nicolson格式。

中心差分格式和Samarskii格式是显式格式，所以很适合于并行计算，但由于稳定性条件的限制，必须采用非常小的时间步长来计算。Crank-Nicolson格式是隐式格式，它是无条件稳定的，但在每一时间层上要求解线性方程组，实现并行计算有一定困难。

中心差分格式的优点是简单易算，但由于截断误差为O(h2)，又仅当

2v1，时才稳定和收敛，所以想要算得略为精确一点，就要缩小。v222ah并且注意到最大为v加工作量。

h2，若h缩小一半，就要缩小为四分之一，这将大大增

Samarskii格式的截断误差也为O(h2)，稳定性条件相对于中心差分格式

v1ah有所放松，为，但计算量稍多于中心差分格式。 221ah/2v2hCrank-Nicolson格式的截断误差为O(2h2)，它关于时间以及空间均为二阶精度，但该格式形成的系数矩阵不一定满足对角占优条件，且计算量大。

鉴于本人的理论水平有限，论文中有些地方还需要加以完善。对于对流扩散方程的差分格式，在以后的理论发展中，还会提出一些具有更好性质的差分方法，使对流扩散方程的差分方法得到更加完善的补充。

. 专业资料分享 .