1 引言

数值模式和数值模拟是定量研究天气和气候变化的主要工具，但是由于数学模型中各变量之间复杂的非线性作用关系，导致模拟结果中存在着不确定性。政府间气候变化专门委员会报告中明确指出正确分析模式的不确定性是IPCC报告的主要任务之一（Mastrandrea et al., 2010）。这些不确定性的来源既包含物理参数的选取、模式的框架和观测的误差，也包括了计算误差的影响（Von Neumann and Goldstine, 1947）。计算误差对数值模拟结果的影响可以从复杂的大气环流模式（王鹏飞等, 2007）、耦合模式（陈显尧等, 2008）运行结果看出，也可以从简单的混沌动力系统（Li et al., 2000; Liao, 2009; Wang et al., 2014）、准地转模式（Teixeira et al., 2007）的数值试验得到验证。因此，如何采取有效的方法来控制计算误差的增长，对长时间的数值计算和精确的数值模拟至关重要。

大气和海洋流体力学模式中有许多形如$\frac{{\partial F}}{{\partial t}} = \mathit{\boldsymbol{LF}}$的发展方程，其中L为包含F以及F关于空间变量导数的算子，使用数值模式求解这类方程组时，计算误差的大小主要取决于空间和时间方向的算法精度。空间导数的精确计算已经有过许多研究，早期人们认为2阶的空间差分格式已经可以满足使用上的要求，随着研究的深入，逐渐认识到直接模拟湍流等复杂流动，高阶的算法是必要的（Lele, 1992），已经有许多研究者使用高精度算法成功的模拟了复杂流体的运动（Lele, 1992; Ma and Fu, 1996）。但在时间积分方案的选取中多采用较低的阶数（3阶左右；季仲贞和王斌, 1994），这就导致了一些模式计算中空间、时间计算精度不匹配，进而影响了模式长时间计算的准确性。尤其是，Tal-Ezer（1986, 1989）、吴声昌和刘小清（1996）等研究表明，只有整体提高算法的阶数，才能使数值模拟中的总误差得到控制，因此，有必要研究适用于大气科学数值模式的时—空匹配的高阶算法，以减小天气和气候数值模拟中由于计算误差导致的不确定性。

使用数值方法求解依赖于时间的常微分方程组或偏微分方程组时，主要是进行空间差分和时间积分两类计算。减小差分方法空间计算误差的途径有加大算法阶数和减小网格尺度两种，但是对于三维的复杂流体力学系统，算法阶数的增加往往是困难的。而网格尺度减小一半时，往往意味着计算量以2⁴方式增加。在大气科学领域前人也进行了高阶算法的研究（季仲贞和王斌, 1994），高阶的差分格式往往存在一个问题，当计算精度提高了，其计算稳定性就降低了，季仲贞等（1999）使用守恒型的差分格式的办法，在计算精度和稳定性两方面取得了较好的平衡。高阶的空间差分算法有很多种实现，如紧致差分格式（Lele, 1992; Ma and Fu, 1996）已有一些应用。Li（2005）提出了一种空间导数的高阶计算方法，它的优点是采用显式计算，实现简单，而且省时，可以很方便的获得10阶甚至更高的阶数。冯涛和李建平（2007）利用此方法研究了一维平流方程、无粘Burgers方程（Hopf, 1950）、正压涡度方程的高精度差分格式，发现6阶空间差分算法的计算效果最好，而且增加的计算量并不多。但随着阶数的进一步增加，从7阶到10阶，计算结果改进的效果不明显，而且7阶以上的差分格式计算结果反而有变差的趋势。

减小时间积分计算误差也有两个途径：减小计算的时间步长和加大算法的阶数。缩小计算所用的时间步长是减小时间积分中误差的最简单方法，Li et al.（2000）和Teixeira et al.（2007）都对此做了研究，他们发现了计算结果对时间步长的敏感性。特别是获得混沌系统长时间数值解的方法不仅依赖于计算的步长，还与算法阶数和浮点计算精度有关。此外，还受制于实际消耗的计算时间，Wang et al.（2012）比较了4阶Runge-Kutta（RK）方法（简称RK4）和Taylor高阶算法在计算时间上的花费，结果表明：高阶算法能够将求解Lorenz方程花费的墙钟时间指数型的降低，高阶算法比使用最优步长计算的低阶算法更具效率。

上述分析表明通过增加算法阶数来减小时间积分时的计算误差是非常必要的。虽然冯涛和李建平（2007）利用Li（2005）方法进行了大量数值试验研究，但是他们的研究中存在两个问题，一是他们只在空间方向上使用了较高阶的算法，而时间积分上只采用了3阶左右的RK法，有可能使计算结果中的误差主要受时间积分误差的影响，从而导致算法不能取得更好的计算结果；二是求解时选取的初值出现不连续空间导数或出现空间导数无穷大等使用差分方法无法适用的问题，因而导致高阶精度算法无法获得比低阶精度算法更好的计算结果，这并不表示高阶精度算法对那些空间导数处处连续且有限的问题也无法提高计算的精度。

为了解决这些问题，Wang（2017）给出一类发展型偏微分方程的高阶Taylor-Li计算方案（简称TLW），它的特点是可以调节时间积分方案的阶数，完成从3阶到30阶甚至更高阶数的算法。Wang（2017）发现时间积分方案的阶数大于3之后，对应的最优空间差分精度阶数可以比6阶提高很多，这就解释了冯涛和李建平（2007）研究中6阶以上空间差分格式对结果无改进的现象，是由于没有使用足够高精度的时间积分方案引起的。

Wang（2017）的研究表明：当算法的时间积分阶数达到5时，一维线性平流方程的算例中相匹配的空间差分阶数为10阶，误差为10⁻⁶量级。而无粘Burgers的算例，与5阶时间积分相匹配的空间差分阶数为16阶，误差为10⁻¹⁴量级。因而，对于一些需要高精度，但又不要求极致的高精度解时，5、6阶的时间积分方案就够用了。相比于Taylor-Li算法，6阶以内的RK方法编程实现要简单，因此，本文将冯涛和李建平（2007）所用的RK3算法，改进为2~6阶可调的RK方法，配合Li高阶微分公式形成Runge-Kutta-Li（简称RKL）算法。通过多组数值试验，研究RKL方法中时间积分阶数对计算误差的影响，评估总体的计算性能，以期找到一种能够对复杂方程适用的简易高阶算法方案，当使用Taylor-Li算法编程代价过高时，替代Taylor-Li算法进行偏微分方程的高精度计算。

2 RK方法与Li空间差分算法联合的RKL算法 2.1 2~6阶RK与Li联合算法

一维平流方程的形式为

$\frac{{\partial u}}{{\partial t}} + \frac{{\partial u}}{{\partial x}} = 0$,

(1)

其中，u为平流变量，可以用一个任意高阶精度的算法来求解。RK方法是一个经典的算法，常被用于求解常微分方程（ODEs），许多文献都对此方法如何微分方程系统进行了介绍（Hairer et al., 2000; Butcher, 2008）。

2阶RK方法如下：

$\left\{ \begin{array}{l} {Z_1} = {z^n},\\ {Z_2} = {z^n} + \frac{1}{2}\tau f({Z_1}),\\ {Z^{n + 1}} = {z^n} + \tau f({Z_2}), \end{array} \right.$

(2)

3阶RK方法（Heun方法）如下：

$\left\{ \begin{array}{l} {Z_1} = {z^n},\\ {Z_2} = {z^n} + \frac{1}{3}\tau f\left( {{Z_1}} \right),\\ {Z_3} = {z^n} + \frac{2}{3}\tau f\left( {{Z_2}} \right),\\ {Z^{n + 1}} = {z^n} + \frac{1}{4}\tau f\left( {{Z_1}} \right) + \frac{3}{4}\tau f\left( {{Z_3}} \right), \end{array} \right.$

(3)

4阶RK方法如下：

$\left\{ \begin{array}{l} {Z_1} = {z^n},\\ {Z_2} = {z^n} + \frac{1}{2}\tau f\left( {{Z_1}} \right),\\ {Z_3} = {z^n} + \frac{1}{2}\tau f\left( {{Z_2}} \right),\\ {Z_4} = {z^n} + \tau f\left( {{Z_3}} \right),\\ {Z^{n + 1}} = {z^n} + \frac{1}{6}\tau \left[ {f\left( {{Z_1}} \right) + 2f\left( {{Z_2}} \right) + 2f\left( {{Z_3}} \right) + f\left( {{Z_4}} \right)} \right]. \end{array} \right.$

(4)

当RK算法阶数提升时，计算步数增加，方程的个数也增加明显，为了简便书写，可将其写为系数表的形式，例如RK4写成系数表的形式为

$\begin{array}{*{20}{l}} {\begin{array}{*{20}{c}} 0\\ {\frac{1}{2}}\\ {\frac{1}{2}}\\ 1 \end{array}\left| {\underline {{\mkern 1mu} \begin{array}{*{20}{l}} {}&{}&{}&{}\\ {\frac{1}{2}}&{}&{}&{}\\ 0&{\frac{1}{2}}&{}&{}\\ 0&0&1&{} \end{array}{\mkern 1mu} } } \right.}\\ {\begin{array}{*{20}{c}} {\;\;\;\frac{1}{6}}&{\frac{1}{3}}&{\frac{1}{3}}&{\frac{1}{6}} \end{array}} \end{array}$ .

(5)

其中，第1行表示：${Z_1} = {z^n}$；第2行表示：${Z_2} = {z^n} + \frac{1}{2}\tau f\left({{Z_1}} \right)$，竖线右侧的第1列（$\frac{1}{2}$）为$\tau f\left({{Z_1}} \right)$的系数；第3行表示：${Z_3} = {z^n} + \frac{1}{2}\tau f\left({{Z_2}} \right)$，竖线右侧的第1和第2列（0，$\frac{1}{2}$）分别为$\tau f\left({{Z_1}} \right)$、$\tau f\left({{Z_2}} \right)$的系数；第4行表示：${Z_4} = {z^n} + \tau f\left({{Z_3}} \right)$，竖线右侧的第1、2、3列（0，0，1）分别为$\tau f\left({{Z_1}} \right)$、$\tau f\left({{Z_2}} \right)$、$\tau f\left({{Z_3}} \right)$的系数；依此类推。第5行（横线下）表示：${Z^{n + 1}} = {z^n} + \frac{1}{6}\tau [f\left({{Z_1}} \right) + 2f\left({{Z_2}} \right) + $ $2f\left({{Z_3}} \right) + f\left({{Z_4}} \right)]$，横线下为计算下一步数值${Z^{n + 1}}$，横线上的计算均为临时变量值。

5阶RK方法如下（Butcher, 2008；这里仅给出系数）：

$\begin{array}{r} \begin{array}{*{20}{c}} 0\\ {\frac{1}{4}}\\ {\begin{array}{*{20}{c}} {}\\ {} \end{array}\frac{1}{4}\begin{array}{*{20}{c}} {}\\ {} \end{array}}\\ {\begin{array}{*{20}{c}} {}\\ {} \end{array}\frac{1}{2}\begin{array}{*{20}{c}} {}\\ {} \end{array}}\\ {\begin{array}{*{20}{c}} {}\\ {} \end{array}\frac{3}{4}\begin{array}{*{20}{c}} {}\\ {} \end{array}}\\ {\begin{array}{*{20}{c}} {}\\ {} \end{array}1\begin{array}{*{20}{c}} {}\\ {} \end{array}} \end{array}\left| \!{\underline {{\mkern 1mu} {\begin{array}{*{20}{c}} {}&{}&{}&{}&{}&{}\\ {\frac{1}{4}}&{}&{}&{}&{}&{}\\ {\frac{1}{8}}&{\frac{1}{8}}&{}&{}&{}&{}\\ 0&0&{\begin{array}{*{20}{c}} {}\\ {} \end{array}\frac{1}{2}}&{}&{}&{}\\ {\begin{array}{*{20}{c}} {}\\ {} \end{array}\frac{3}{{16}}}&{ - \frac{3}{8}}&{\begin{array}{*{20}{c}} {}\\ {} \end{array}\frac{3}{8}}&{\begin{array}{*{20}{c}} {}\\ {} \end{array}\frac{9}{{16}}}&{}&{}\\ { - \frac{3}{7}}&{\frac{8}{7}}&{\frac{6}{7}}&{ - \frac{{12}}{7}}&{\begin{array}{*{20}{c}} {}\\ {} \end{array}\frac{8}{7}}&{\begin{array}{*{20}{c}} {} \end{array}} \end{array}\begin{array}{*{20}{c}} {} \end{array}} \,}} \right. \\ \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {}\\ {} \end{array}\frac{7}{{90}}}&{\begin{array}{*{20}{c}} {}\\ {} \end{array}0}&{\begin{array}{*{20}{c}} {}\\ {} \end{array}\frac{{32}}{{90}}}&{\begin{array}{*{20}{c}} {}\\ {} \end{array}\frac{{12}}{{90}}}&{\frac{{32}}{{90}}}&{\frac{7}{{90}}} \end{array} \end{array} $ .

(6)

6阶RK方法如下（Butcher, 2008；这里仅给出系数）：

$\begin{array}{r} \begin{array}{*{20}{c}} 0\\ {\frac{1}{3}}\\ {\frac{2}{3}}\\ {\frac{1}{3}}\\ {\frac{5}{6}}\\ {\frac{1}{6}}\\ {\begin{array}{*{20}{c}} {}\\ {} \end{array}1\begin{array}{*{20}{c}} {}\\ {} \end{array}} \end{array}\left| \!{\underline {{\mkern 1mu} {\begin{array}{*{20}{c}} {}&{}&{}&{}&{}&{}&{}\\ {\frac{1}{3}}&{}&{}&{}&{}&{}&{}\\ 0&{\frac{2}{3}}&{}&{}&{}&{}&{}\\ {\frac{1}{{12}}}&{\frac{1}{3}}&{ - \frac{1}{{12}}}&{}&{}&{}&{}\\ {\frac{{25}}{{48}}}&{ - \frac{{55}}{{44}}}&{\frac{{35}}{{48}}}&{\frac{{15}}{8}}&{}&{}&{}\\ {\frac{3}{{20}}}&{ - \frac{{11}}{{24}}}&{ - \frac{1}{8}}&{\frac{1}{2}}&{\frac{1}{{10}}}&{}&{}\\ { - \frac{{261}}{{260}}}&{\frac{{33}}{{13}}}&{\frac{{43}}{{156}}}&{ - \frac{{118}}{{39}}}&{\frac{{32}}{{195}}}&{\frac{{80}}{{39}}\begin{array}{*{20}{c}} {}\\ {} \end{array}}&{} \end{array}\begin{array}{*{20}{c}} {} \end{array}} \,}} \right. \\ \begin{array}{*{20}{c}} {\frac{{13}}{{200}}}&{\begin{array}{*{20}{c}} {}\\ {} \end{array}0\begin{array}{*{20}{c}} {}\\ {} \end{array}}&{\begin{array}{*{20}{c}} {}\\ {} \end{array}\frac{{11}}{{40}}\begin{array}{*{20}{c}} {}\\ {} \end{array}\begin{array}{*{20}{c}} {}\\ {} \end{array}}&{\frac{{11}}{{40}}}&{\begin{array}{*{20}{c}} {}\\ {} \end{array}\frac{4}{{25}}}&{\frac{4}{{25}}}&{\frac{{13}}{{200}}} \end{array} \end{array}$ .

(7)

需要注意的是，3阶以上的RK方法，系数形式并不唯一，这里选用了其中一种常用的形式。空间方向采用任意阶精度的Li算法（Li, 2005），计算公式为

${f_y}^{\left(m \right)}({y_i}) = \frac{1}{{{h^m}}}\sum\limits_{j = 0}^n {d_{n + 1, i, j}^{\left(m \right)}f({y_j})} $,

(8)

其中，$d_{n + 1, i, j}^{\left(m \right)}$的含义和计算方法参考Li（2005），本文选$m \equiv 1$，因此，n阶精度需要(n+1)个点即可；需要注意的是，此公式中的${y_i}$位置为其在从0到n个计算点之间的相对位置，而不是$u({x_i})$中实际网格的绝对位置；此公式的精度为$\left({n - m + 1} \right)$阶。

当$x$方向使用N个网格点时，$\frac{\partial }{{\partial x}}(u)$可以使用公式（8）计算，且能保证其达到n阶精度，使用n+1个格点（n+1为从N个网格点中选取出的格点数，n≤N－1）。例如计算x_i点时，如果n+1为奇数，这n+1个格点可选为$({x_{i - n/2}}, \cdots, {x_i}, \cdots, {x_{i + n/2}})$（${x_{i - n/2}}$对应${y_0}$，…，${x_i}$对应${y_{n/2}}$，…，${x_{i + n/2}}$对应${y_n}$），使x_i位于中间点；如果n+1为偶数，这n个格点可选为$({x_{i - \left({n + 1} \right)/2}}, \cdots, {x_i}, \cdots, {x_{i + \left({n + 1} \right)/2 - 1}})$，使x_i尽量靠近中间点；如果x_i靠近格点的左边界或右边界，而无法将其放于中心位置时，可以选取在左侧或右侧边界上能够覆盖x_i的连续n个格点即可。对于能使用周期边界条件而避免边界效应的方程，应尽量选取周期边界条件，这样可以保证x_i尽量靠近计算的中间点。

当$m \equiv 1$时，空间方向1阶导数的计算公式为

${f_y}^{(1)}({y_i}) = \frac{1}{h}\sum\limits_{j = 0}^n {d_{n + 1, i, j}^{(1)}f({y_j})} $.

(9)

计算时，空间导数的精度为n阶精度，为了保证空间导数能达到一定的精度，至少要求使用（n+1）个格点。例如取N＝64，若n＝5阶，则每次计算导数时需要6个格点，依此类推。具体的计算公式如下：

$ d_{2,0,1}^{\left( 1 \right)} = 1,\;\;\;d_{2,1,0}^{\left( 1 \right)} = - 1, $

(10)

$ d_{n + 1, i, j}^{(1)} = \frac{{{{(- 1)}^{(1 - j)}}a_{n - 1, i, j}^{(0)}}}{{j!(n - j)!}}, \;\;\;(当i \ne j) $

(11)

$ d_{n + 1, i, i}^{(1)} = - \sum\limits_{j = 0, j \ne i}^n {d_{n + 1, i, j}^{(1)}}, \;\;(当i \ne j)$

(12)

$ \begin{array}{l} a_{n - 1, i, j}^{\left(0 \right)} = {a_0}(- i, \cdots, k - i, \cdots, n - i), (k \ne i, k \ne j)\\ {\rm{ = }}(- i) \cdot (\cdots) \cdot (k - i) \cdot (\cdots) \cdot (n - i), (k \ne i, k \ne j) \end{array} $

(13)

存储空间的安排为$d[m][n][i][j]$，由于$m \equiv 1$，所以可以减少一维，变为${d^{\left(1 \right)}}[n][i][j]$，一般n与所需的时间精度M有关，一但n选定，可以将存储空间变为二维数组$d_{n + 1}^{(1)}[i][j]$，最大需要$(n + 1) \times (n + 1)$的内存数组空间。

求解时每步中逐个格点先计算空间差，后使用RK法得到下一时刻的数值，全部格点计算完毕进入下一步积分，直到目标时刻，此即RKL方法对一维线性平流方程的求解过程。

2.2 RKL算法的多精度实现

高阶算法在计算时还会碰到一个问题就是舍入误差，因为计算所用的时空精度阶数都可能超过10阶，Wang（2017）试验中，解的数值量级为$O(1)$，如果仅采用双精度计算，会发现高阶算法的绝对误差有时徘徊在10⁻¹⁵，这恰好是双精度计算时的相对误差极限，所以为了观察超高阶算法的计算效果，采用多精度（Multiple-Precision，简称MP）计算是必要的。本文使用了MP库，采用1024二进制位精度，相当于200位以上的十进制有效数字，足以区分绝对误差小到10⁻²⁰⁰的数值解。需注意MP版的程序比对应的使用双精度的程序慢约100倍左右，但仍可以较快的完成计算，相比于计算精度的极大提升，这些计算速度的减低是可以接受的，而且有并行计算的办法使之提升（Wang et al., 2014）。

3 RKL方法的线性平流试验

试验1：求解平流方程公式（1），使用Gauss波初始条件：$u(x)\left| {_{t = 0}} \right. = {e^{ - 400{{(x - 0.5)}^2}}}$。

仿照（Takacs, 1985）的做法，定义总误差为$E = \sqrt {\frac{1}{N}\sum\limits_{i = 1}^N {{{\left({{u_D} - {u_T}} \right)}^2}} } $，其中${u_D}$表示数值解，${u_T}$表示理论解，N为x方向的格点数，这里我们对平均误差方差取了根号，以保证总误差和原变量具有同样的量纲。

计算时x方向的空间步长为$h = 1/200$，时间步长为$\tau = 1/400$，计算区域为[0, 1]，计算400步，试验结果如图 1b，从图 1b可见，在RK算法阶数为2时，计算误差较大，而且即使空间差分阶数提升，误差也无明显变化；若M＝3阶，n达到6阶时，计算误差与n大于6阶时的计算误差相差不多；M＝4阶时，n达到7阶左右计算误差接近平缓变化；M＝5阶时，n达到9阶左右计算误差接近平缓变化；M＝6阶时，n达到10阶之后计算误差接近平缓变化。这些结果说明，本试验中总计算误差受时间积分精度的影响较大，当时间方向精度大于3阶时，空间方向的格式精度可以超过6阶，这与使用Taylor-Li方法（Wang, 2017）的结果一致。

对方波型初值也进行了试验[图略，与Wang（2017）的图 1相似]，高阶精度（时间方向）算法只能获得和3阶精度差不多的结果。而对于光滑的初值（图 1），高阶RK算法可以取得较好的计算效果，基本上随着RK算法阶数的提升，Li微分方案对应的有效阶数也在加大，总误差在逐渐减小。产生这种差别的原因是方波型初值不连续，导致$x$方向的各阶导数也不连续，使计算效果变差；而光滑型初值所对应的解各阶导数均存在且连续，因此可以获得较好的结果。

4 RKL方法的无粘Burgers方程试验

作为非线性问题的示例，分析一维无粘性Burgers方程：

$\frac{{\partial u}}{{\partial t}} + u\frac{{\partial u}}{{\partial x}} = 0$.

(14)

理论解的形式：如果初始场为$u(x, 0) = f(x)$，那么

$u(x, t) = f(x - u(x, t)t)$,

(15)

一定是方程（14）的解，这可以通过如下关系式得到：$\frac{{\partial u}}{{\partial t}} = \frac{{\partial f(x - u(x, t)t)}}{{\partial t}} = f'\left({ - u - t\frac{{\partial u}}{{\partial t}}} \right)$，即

$\frac{{\partial u}}{{\partial t}} = \frac{{ - uf'}}{{1 + tf'}}$,

(16)

$\frac{{\partial u}}{{\partial x}} = \frac{{\partial f\left( {x - u\left( {x,t} \right)t} \right)}}{{\partial x}} = f'\left( {1 - t\frac{{\partial u}}{{\partial x}}} \right) $, 即

$\frac{{\partial u}}{{\partial x}} = \frac{{f'}}{{1 + tf'}}$ .

(17)

将方程（16）和（17）代入方程（14），可见等式成立。方程（16）和（17）可以用于检验程序中一阶数值微分的精确性，而通过公式（15），使用隐函数求解的办法，可以得到任意时刻的理论解，以对数值解进行评估。

试验2：求解无粘性Burgers方程$\frac{{\partial u}}{{\partial t}} + u\frac{{\partial u}}{{\partial x}} = 0$，初始条件为：$u(x, 0) = - \sin (x)$。

试验时，计算区域选为$\left[ { - \pi, \pi } \right]$，网格数N=800，时间步长为$\tau = 0.001$，侧边界条件为周期边界条件。计算时进行检验计算误差的时刻选为：t＝0.8，共计算800步。

图 2给出了试验2的结果，由于计算时绝对误差会小于10⁻¹⁵，因此试验使用的是多精度的程序。从图 2c可以看出，在RK算法阶数为2时，计算误差较大，而且即使空间差分阶数提升，误差也无明显变化；当时间精度为3阶时，计算误差在空间精度达到6阶后就基本保持不变了；当时间精度为4、5、6阶时，有效的空间差分精度阶数分别可达12、25、40。这个试验说明了对非线性Burgers方程，只要时间积分精度阶数足够高，空间方向的差分精度阶数可以超过6阶。

对另外一组初值$u(x, 0) = - \arctan (x)$也进行了试验[图略，与Wang（2017）的图 4相似]，高阶精度（时间方向）算法只能获得和三阶精度差不多的结果。

综上所述，对于非线性无粘Burgers方程，高阶精度算法能否获得好的计算结果，除了受初值场分布形势的影响，还与计算的目标时刻有联系。当解的各阶导数连续且未出现$\infty $数值时，高阶精度的算法是有效的；若出现了导数间断、或导数为$\infty $，就会碰到冲击波这样的问题，此时高阶精度的算法也无法获得很精确的数值解。

值得一提的是，若数值解能够具有周期边界条件，求解精度会比非周期边界条件的高，这主要是周期边界条件满足时，总可以将目标格点放置于差分公式中的中点附近，计算时误差小。而非周期的边界条件，在边界点以及紧靠边界点处，即使能够用n阶精度的差分格式，目标格点也会被放置于偏离n个差分点中心的位置，实际计算表明这些位置的差分精度较差。

5 结论

本文利用Li提出的高阶显式微分公式，结合Runge-Kutta时间积分方案，实现了求解含时偏微分方程的Runge-Kutta-Li高阶算法格式。对一维平流方程，通过比较理论解和各阶计算格式的数值解，研究了计算误差随时间积分阶数的变化情况。结果表明：对于一些光滑、周期边界条件的初值，如果时间方向的积分精度为3阶，计算误差在空间差分阶数到达6阶后减小的不明显。而当时间积分的阶数大于3，例如选为4、5、6阶时，此时空间差分阶数超过6之后，计算误差仍有减小的趋势，相应的有效空间差分阶数可以超过6阶，说明了冯涛等所发现的计算误差在空间差分精度达到6阶出现饱和的现象是由于没有足够高的时间积分方案配合而引起的。对于方波这样的初值，本文计算结果与以前的研究一致。

对非线性的Burgers方程，一个光滑的且满足周期边界条件的算例表明，Runge-Kutta-Li格式能够取得较好的计算效果。相比于冯涛和李建平（2007）使用的3阶时间积分方案，本文实现了2~6阶的RK时间积分算法，而相应的空间差分精度也能从6阶拓展到30阶以上，因而计算精度也极大提升。

虽然可以通过减小时间步长来改进时间积分的精度，但是它的效率不如增加阶数高，这在以往的研究中已经被论证（Wang et al., 2012）。因此，相比于以往的仅增加空间差分精度的高阶算法，本文在时空两个方向提高算法的精度阶数，减少了总体计算误差。为了控制高阶算法计算过程中浮点舍入误差的影响，使用多精度计算的工具库，实现了可以求解平流方程和Burgers方程的RKL格式，此方法可以方便的拓展到其它类似方程的求解中。

本文实现的RKL算法程序可以灵活的指定时间积分阶数（2~6阶）、空间计算的阶数，还能够方便的调整空间格距、时间积分的步长，因此适用于研究算法阶数、时空步长和计算误差之间的复杂关系。本研究仅用到最高6阶RK方法，对于6阶以上的RK方法，实现较为复杂，特别是RK方法中的系数计算公式，若用到MP库计算则必须为完全准确的（即不能查算好的系数表，要在程序里自己算），因而会增加程序实现的困难，有兴趣的研究者可以仿照本文方法，进行更高阶数的试验。