大气科学  2019, Vol. 43 Issue (1): 99-106 PDF
Runge-Kutta算法与Li差分法不同阶数配合对计算精度影响研究

1 中国科学院大气物理研究所季风系统研究中心, 北京 100190
2 中国科学院大气物理研究所大气科学和地球流体力学数值模拟国家重点实验室, 北京 100029
3 中国科学院大学, 北京 100049
4 中国科学技术大学地球和空间科学学院, 合肥 230026
5 中国科学院东亚区域气候-环境重点实验室, 北京 100029

A Study on the Precision of Runge-Kutta Method with Various Orders of Li Difference Scheme
WANG Pengfei1,2, CHU Pingxiang1,3, WANG Lizhi5, ZHOU Renjun4, HUANG Gang2
1 Center for Monsoon System Research, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100190
2 State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics(LASG), Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100029
3 University of Chinese Academy of Sciences, Beijing 100049
4 School of Earth and Space Sciences(SESS), University of Science and Technology of China, Hefei 230016
5 Key Laboratory of Regional Climate?Environment for Temperate East Asia(TEA), Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100029
Abstract: We implement the hybrid Runge-Kutta-Li (RKL) scheme for the purpose to take full advantage of Li's high order spatial differential method (Li, 2005). A set of numerical experiments has been conducted to analyze how the computation error is affected by the order of integration scheme. The results of the linear advection equation indicate that with the square-wave type initial values, the scheme can only obtain a third-order accuracy. However, for the Gaussian function type of initial values, the scheme can obtain a better result. The fifth (sixth) order Runge-Kutta (RK) integration scheme corresponds to 9th (10th) order Li's difference scheme in spatial direction and the total error can be controlled within 10-7 (10-8). The order of Li's scheme tends to increase while the RK order increases, and the total error gradually decreases. When we compute the nonlinear Burgers' equation, whether the RKL scheme can obtain good results is not only dependent on the form of the initial field, but also related to the target computation time. When the derivative is continuous (and infinite value does not appear) at the target observation time, 4th-6th order RKL scheme is effective. On the contrary, if the derivative is discontinuous, or the derivative tends to infinity, the RKL scheme cannot obtain high-precision numerical solution. In this case (Burgers' with smooth initial), the order of Li's scheme still increases while the RK order increases, but the relation between them shows a nonlinear tendency (which can be specified through some fitting methods). The results indicate that when the order of time integral is more than three, the corresponding optimal spatial difference order can be higher than six. This result confirms the finding of previous studies that the order of spatial difference above six makes no improvement to the results due to the lack of high-precision time integral scheme. Compared with Taylor-Li (Wang, 2017) scheme, the 5th-6th order RKL scheme is easier to program and can yield more precise results than the third-order scheme. To conclude, the high order RKL scheme can be applied to some complicated types of partial differential equations and is valuable for many other similar computation cases.
Keywords: Runge-Kutta-Li scheme    High-order    Burgers' equation
1 引言

Wang（2017）的研究表明：当算法的时间积分阶数达到5时，一维线性平流方程的算例中相匹配的空间差分阶数为10阶，误差为10−6量级。而无粘Burgers的算例，与5阶时间积分相匹配的空间差分阶数为16阶，误差为10−14量级。因而，对于一些需要高精度，但又不要求极致的高精度解时，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)

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)

 $\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)

$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)

 $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)

2.2 RKL算法的多精度实现

3 RKL方法的线性平流试验

 图 1 RKL（Runge-Kutta-Li）方法的线性平流试验：（a）Gauss波初始场；（b）计算误差随空间精度阶数的变化。（b）中横坐标为空间精度阶数，纵坐标为误差取对数（log10），粉、蓝、红、绿、黑色分别代表时间精度为2、3、4、5、6阶 Figure 1 The experiments of linear advection equation by RKL (Runge-Kutta-Li) method: (a) Gaussian-type initial condition; (b) error versus spatial difference order. In (b), the abscissa is the spatial difference order, the ordinate is the logarithm of the error (log10), and the purple, blue, red, green and black curves denote the second-order, third-order, fourth-order, fifth-order and sixth-order time-integration schemes, respectively

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

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

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

 $\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)

 图 2 RKL方法的无粘性Burgers方程试验：（a）t=0时的初始场；（b）t=0.8时u的理论解；（c）计算误差随空间精度阶数的变化。（c）中横坐标为空间精度阶数，纵坐标为误差取对数（log10），粉、蓝、红、绿、黑色分别代表时间精度为2、3、4、5、6阶 Figure 2 The experiments of the nonlinear Burgers' equation by RKL method: (a) Initial u at t=0; (b) analytical solution of u at t=0.8; (c) error versus spatial difference order. In (c), the abscissa is the spatial difference order and the ordinate is the logarithm of error (log10); purple, blue, red, green and black curves correspond to the second-order, third-order, fourth-order, fifth-order and sixth-order time-integration schemes, respectively

5 结论

 Butcher J C. 2008. Numerical Methods for Ordinary Differential Equations[M]. 2nd ed. England: John Wiley & Sons: 463pp. 陈显尧, 宋振亚, 赵伟, 等. 2008. 气候模式系统模拟结果的不确定性分析[J]. 海洋科学进展, 26(2): 119-125. Chen Xianyao, Song Zhenya, Zhao Wei, et al. 2008. Uncertainity analysis of results simulated by climate model system[J]. Advances in Marine Science (in Chinese), 26(2): 119-125. DOI:10.3969/j.issn.1671-6647.2008.02.001 冯涛, 李建平. 2007. 高精度迎风偏斜格式的比较与分析[J]. 大气科学, 31(2): 245-253. Feng Tao, Li Jianping. 2007. A comparison and analysis of high order upwind-biased schemes[J]. Chinese Journal of Atmospheric Sciences (in Chinese), 31(2): 245-253. DOI:10.3878/j.issn.1006-9895.2007.02.06 Hairer E, Nørsett S P, Wanner G. 2000. Solving Ordinary Differential Equations I:Nonstiff Problems[M]. Berlin Heidelberg: Springer: 528pp. Hopf E. 1950. The partial differential equation ut + uux=μxx[J]. Commun. Pure Appl. Math., 3(3): 201-230. DOI:10.1002/cpa.3160030302 季仲贞, 王斌. 1994. 一类高时间差分精度的平方守恒格式的构造及其应用检验[J]. 自然科学进展, 4(2): 149-157. Ji Zhongzhen, Wang Bin. 1994. Construction and application test of a kind of high precision scheme with square-conservation[J]. Progress in Natural Science (in Chinese), 4(2): 149-157. DOI:10.3321/j.issn:1002-008X.1994.02.004 季仲贞, 李京, 王斌. 1999. 紧致平方守恒格式的构造和检验[J]. 大气科学, 23(3): 323-332. Ji Zhongzhen, Li Jing, Wang Bin. 1999. Construction and test of compact scheme with square-conservation[J]. Chinese Journal of Atmospheric Sciences (in Chinese), 23(3): 323-332. DOI:10.3878/j.issn.1006-9895.1999.03.08 Lele S K. 1992. Compact finite difference schemes with spectral-like resolution[J]. J. Comput. Phys., 103(1): 16-42. DOI:10.1016/0021-9991(92)90324-R Li J P. 2005. General explicit difference formulas for numerical differentiation[J]. J. Comput. Appl. Math., 183(1): 29-52. DOI:10.1016/j.cam.2004.12.026 Li J P, Zeng Q C, Chou J F. 2000. Computational uncertainty principle in nonlinear ordinary differential equations (I)——Numerical results[J]. Science in China (Series E), 43(5): 449-460. Liao S J. 2009. On the reliability of computed chaotic solutions of non-linear differential equations[J]. Tellus A, 61(4): 550-564. DOI:10.1111/j.1600-0870.2009.00402.x Ma Y W, Fu D X. 1996. Super compact finite difference method (SCFDM) with arbitrary high accuracy[J]. Comput. Fluid Dyn. J., 5(2): 259-276. Mastrandrea M D, Field C B, Stocker T F, et al. 2010. Guidance Note for Lead Authors of the IPCC Fifth Assessment Report on Consistent Treatment of Uncertainties[C]. Jasper Ridge, CA, USA: Intergovernmental Panel on Climate Change. Takacs L L. 1985. A two-step scheme for the advection equation with minimized dissipation and dispersion errors[J]. Mon. Wea. Rev., 113(6): 1050-1065. DOI:10.1175/1520-0493(1985)113<1050:ATSSFT>2.0.CO;2 Tal-Ezer H. 1986. Spectral methods in time for hyperbolic equations[J]. SIAM J. Numer. Anal., 23(1): 11-26. DOI:10.1137/0723002 Tal-Ezer H. 1989. Spectral methods in time for parabolic problems[J]. SIAM J. Numer. Anal., 26(1): 1-11. DOI:10.1137/0726001 Teixeira J, Reynolds C A, Judd K. 2007. Time step sensitivity of nonlinear atmospheric models:Numerical convergence, truncation error growth, and ensemble design[J]. J. Atmos. Sci., 64(1): 175-189. DOI:10.1175/JAS3824.1 Von Neumann J, Goldstine H H. 1947. Numerical inverting of matrices of high order[J]. Bull. Amer. Math. Soc., 53(11): 1021-1099. DOI:10.1090/S0002-9904-1947-08909-6 Wang P F. 2017. A high-order spatiotemporal precision-matching Taylor-Li scheme for time-dependent problems[J]. Adv. Atmos. Sci., 34(12): 1461-1471. DOI:10.1007/s00376-017-7018-1 王鹏飞, 王在志, 黄刚. 2007. 舍入误差对大气环流模式模拟结果的影响[J]. 大气科学, 31(5): 815-825. Wang Pengfei, Wang Zaizhi, Huang Gang. 2007. The influence of round-off error on the atmospheric general circulation model[J]. Chinese Journal of Atmospheric Sciences (in Chinese), 31(5): 815-825. DOI:10.3878/j.issn.1006-9895.2007.05.06 Wang P F, Li J P, Li Q. 2012. Computational uncertainty and the application of a high-performance multiple precision scheme to obtaining the correct reference solution of Lorenz equations[J]. Numerical Algorithms, 59(1): 147-159. DOI:10.1007/s11075-011-9481-6 Wang P F, Liu Y, Li J P. 2014. Clean numerical simulation for some chaotic systems using the parallel multiple-precision Taylor scheme[J]. Chinese Sci. Bull., 59(33): 4465-4472. DOI:10.1007/s11434-014-0412-5 吴声昌, 刘小清. 1996. KdV方程的时间谱离散方法[J]. 应用数学和力学, 17(4): 357-362. Wu Shengchang, Liu Xiaoqing. 1996. Spectral method in time for KdV equations[J]. Applied Mathematics and Mechanics (in Chinese), 17(4): 357-362.